summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--matlab/algorithms/DART/DARTalgorithm.m196
-rw-r--r--matlab/algorithms/DART/IterativeTomography.m455
-rw-r--r--matlab/algorithms/DART/IterativeTomography3D.m433
-rw-r--r--matlab/algorithms/DART/Kernels.m15
-rw-r--r--matlab/algorithms/DART/MaskingDefault.m17
-rw-r--r--matlab/algorithms/DART/MaskingGPU.m15
-rw-r--r--matlab/algorithms/DART/OutputDefault.m15
-rw-r--r--matlab/algorithms/DART/SegmentationDefault.m15
-rw-r--r--matlab/algorithms/DART/SmoothingDefault.m17
-rw-r--r--matlab/algorithms/DART/SmoothingGPU.m15
-rw-r--r--matlab/algorithms/DART/StatisticsDefault.m21
-rw-r--r--matlab/algorithms/DART/TomographyDefault.m73
-rw-r--r--matlab/algorithms/DART/TomographyDefault3D.m73
-rw-r--r--matlab/algorithms/DART/examples/example1.m88
-rw-r--r--matlab/algorithms/DART/examples/example2.m87
-rw-r--r--matlab/algorithms/DART/examples/example3.m78
-rw-r--r--matlab/algorithms/DART/tools/DARToptimizer.m118
-rw-r--r--matlab/algorithms/DART/tools/DARToptimizerBoneStudy.m118
-rw-r--r--matlab/algorithms/DART/tools/ProjDiffOptimFunc.m29
-rw-r--r--matlab/algorithms/DART/tools/dart_create_base_phantom.m17
-rw-r--r--matlab/algorithms/DART/tools/dart_scheduler.m53
-rw-r--r--matlab/algorithms/DART/tools/rNMPOptimFunc.m35
22 files changed, 680 insertions, 1303 deletions
diff --git a/matlab/algorithms/DART/DARTalgorithm.m b/matlab/algorithms/DART/DARTalgorithm.m
index b7e63b5..85a3ca0 100644
--- a/matlab/algorithms/DART/DARTalgorithm.m
+++ b/matlab/algorithms/DART/DARTalgorithm.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef DARTalgorithm < matlab.mixin.Copyable
@@ -20,7 +19,7 @@ classdef DARTalgorithm < matlab.mixin.Copyable
%----------------------------------------------------------------------
properties (GetAccess=public, SetAccess=public)
- tomography = TomographyDefault(); % POLICY: Tomography object.
+ tomography = IterativeTomography(); % POLICY: Tomography object.
segmentation = SegmentationDefault(); % POLICY: Segmentation object.
smoothing = SmoothingDefault(); % POLICY: Smoothing object.
masking = MaskingDefault(); % POLICY: Masking object.
@@ -29,10 +28,11 @@ classdef DARTalgorithm < matlab.mixin.Copyable
base = struct(); % DATA(set): base structure, should contain: 'sinogram', 'proj_geom', 'phantom' (optional).
- memory = 'yes'; % SETTING: reduce memory usage? (disables some features)
-
- testdata = struct();
-
+ memory = 'no'; % SETTING: reduce memory usage? (disables some features)
+
+ implementation = 'linear'; % SETTING: which type of projector is used ('linear', 'nonlinear')
+ t = 5; % SETTING: # ARMiterations, each DART iteration.
+ t0 = 100; % SETTING: # ARM iterations at DART initialization.
end
%----------------------------------------------------------------------
properties (GetAccess=public, SetAccess=private)
@@ -57,15 +57,25 @@ classdef DARTalgorithm < matlab.mixin.Copyable
methods
%------------------------------------------------------------------
- function this = DARTalgorithm(base)
+ function this = DARTalgorithm(varargin)
% Constructor
- % >> D = DARTalgorithm(base); base is a matlab struct (or the path towards one)
- % that should contain 'sinogram', 'proj_geom'.
-
- if ischar(base)
- this.base = load(base);
+ % >> D = DARTalgorithm(base); [base is a matlab struct that
+ % should contain 'sinogram' and
+ % 'proj_geom']
+ % >> D = DARTalgorithm('base_path'); [path to base struct file]
+ % >> D = DARTalgorithm(sinogram, proj_geom)
+ %
+ narginchk(1, 2)
+ if nargin == 1 && ischar(varargin{1})
+ this.base = load(varargin{1});
+ elseif nargin == 1 && isstruct(varargin{1})
+ this.base = varargin{1};
+ elseif nargin == 2
+ this.base = struct();
+ this.base.sinogram = varargin{1};
+ this.base.proj_geom = varargin{2};
else
- this.base = base;
+ error('invalid arguments')
end
end
@@ -74,7 +84,6 @@ classdef DARTalgorithm < matlab.mixin.Copyable
function D = deepcopy(this)
% Create a deep copy of this object.
% >> D2 = D.deepcopy();
-
D = copy(this);
props = properties(this);
for i = 1:length(props)
@@ -82,7 +91,6 @@ classdef DARTalgorithm < matlab.mixin.Copyable
D.(props{i}) = copy(this.(props{i}));
end
end
-
end
%------------------------------------------------------------------
@@ -91,14 +99,18 @@ classdef DARTalgorithm < matlab.mixin.Copyable
% >> D.initialize();
% Initialize tomography part
- this.tomography.initialize(this);
+ if ~this.tomography.initialized
+ this.tomography.sinogram = this.base.sinogram;
+ this.tomography.proj_geom = this.base.proj_geom;
+ this.tomography.initialize();
+ end
% Create an Initial Reconstruction
if isfield(this.base, 'V0')
this.V0 = this.base.V0;
else
this.output.pre_initial_iteration(this);
- this.V0 = this.tomography.createInitialReconstruction(this, this.base.sinogram);
+ this.V0 = this.tomography.reconstruct2(this.base.sinogram, [], this.t0);
this.output.post_initial_iteration(this);
end
this.V = this.V0;
@@ -113,85 +125,98 @@ classdef DARTalgorithm < matlab.mixin.Copyable
% iterate
function this = iterate(this, iters)
% Perform several iterations of the DART algorithm.
- % >> D.iterate(iterations);
-
+ % >> D.iterate(iterations);
+ if strcmp(this.implementation,'linear')
+ this.iterate_linear(iters);
+ elseif strcmp(this.implementation,'nonlinear')
+ this.iterate_nonlinear(iters);
+ end
+ end
+
+
+ %------------------------------------------------------------------
+ % iterate - linear projector implementation
+ function this = iterate_linear(this, iters)
+
this.start_tic = tic;
for iteration = 1:iters
+ this.iterationcount = this.iterationcount + 1;
+
+ % initial output
+ this.output.pre_iteration(this);
+
+ % update adaptive parameters
+ this.update_adaptiveparameter(this.iterationcount);
+
+ % segmentation
+ this.segmentation.estimate_grey_levels(this, this.V);
+ this.S = this.segmentation.apply(this, this.V);
+
+ % select update and fixed pixels
+ this.Mask = this.masking.apply(this, this.S);
+ this.V = (this.V .* this.Mask) + (this.S .* (1 - this.Mask));
+ F = this.V;
+ F(this.Mask == 1) = 0;
+
+ % compute residual projection difference
+ this.R = this.base.sinogram - this.tomography.project(F);
+ % ART update part
+ this.V = this.tomography.reconstruct2_mask(this.R, this.V, this.Mask, this.t);
+
+ % blur
+ this.V = this.smoothing.apply(this, this.V);
+
+ %calculate statistics
+ this.stats = this.statistics.apply(this);
+
+ % output
+ this.output.post_iteration(this);
+ end
+
+ end
+
+ %------------------------------------------------------------------
+ % iterate - nonlinear projector implementation
+ function this = iterate_nonlinear(this, iters)
+
+ this.start_tic = tic;
+
+ for iteration = 1:iters
this.iterationcount = this.iterationcount + 1;
- %----------------------------------------------------------
- % Initial Output
+ % Output
this.output.pre_iteration(this);
- %----------------------------------------------------------
- % Update Adaptive Parameters
- for i = 1:numel(this.adaptparam_name)
-
- for j = 1:numel(this.adaptparam_iters{i})
- if this.iterationcount == this.adaptparam_iters{i}(j)
- new_value = this.adaptparam_values{i}(j);
- eval(['this.' this.adaptparam_name{i} ' = ' num2str(new_value) ';']);
- disp(['value ' this.adaptparam_name{i} ' changed to ' num2str(new_value) ]);
- end
- end
-
- end
+ % update adaptive parameters
+ this.update_adaptiveparameter(this.iterationcount)
- %----------------------------------------------------------
% Segmentation
this.segmentation.estimate_grey_levels(this, this.V);
this.S = this.segmentation.apply(this, this.V);
- %----------------------------------------------------------
% Select Update and Fixed Pixels
this.Mask = this.masking.apply(this, this.S);
-
this.V = (this.V .* this.Mask) + (this.S .* (1 - this.Mask));
- %this.V(this.Mask == 0) = this.S(this.Mask == 0);
-
- F = this.V;
- F(this.Mask == 1) = 0;
- %----------------------------------------------------------
- % Create Residual Projection Difference
- %this.testdata.F{iteration} = F;
- this.R = this.base.sinogram - this.tomography.createForwardProjection(this, F);
- %this.testdata.R{iteration} = this.R;
+ % ART update part
+ this.V = this.tomography.reconstruct2_mask(this.base.sinogram, this.V, this.Mask, this.t);
- %----------------------------------------------------------
- % ART Loose Part
- %this.testdata.V1{iteration} = this.V;
- %this.testdata.Mask{iteration} = this.Mask;
-
- %X = zeros(size(this.V));
- %Y = this.tomography.createReconstruction(this, this.R, X, this.Mask);
- %this.V(this.Mask) = Y(this.Mask);
- this.V = this.tomography.createReconstruction(this, this.R, this.V, this.Mask);
-
- %this.testdata.V2{iteration} = this.V;
-
- %----------------------------------------------------------
- % Blur
+ % blur
this.V = this.smoothing.apply(this, this.V);
- %this.testdata.V3{iteration} = this.V;
- %----------------------------------------------------------
- % Calculate Statistics
+ % calculate statistics
this.stats = this.statistics.apply(this);
- %----------------------------------------------------------
- % Output
+ % output
this.output.post_iteration(this);
- end % end iteration loop
-
- %test = this.testdata;
- %save('testdata.mat','test');
+ end
+
+ end
+
- end
-
%------------------------------------------------------------------
% get data
function data = getdata(this, string)
@@ -208,9 +233,22 @@ classdef DARTalgorithm < matlab.mixin.Copyable
this.adaptparam_name{end+1} = name;
this.adaptparam_values{end+1} = values;
this.adaptparam_iters{end+1} = iterations;
- end
+ end
%------------------------------------------------------------------
+ % update adaptive parameter
+ function this = update_adaptiveparameter(this, iteration)
+ for i = 1:numel(this.adaptparam_name)
+ for j = 1:numel(this.adaptparam_iters{i})
+ if iteration == this.adaptparam_iters{i}(j)
+ new_value = this.adaptparam_values{i}(j);
+ eval(['this.' this.adaptparam_name{i} ' = ' num2str(new_value) ';']);
+ end
+ end
+ end
+ end
+
+ %------------------------------------------------------------------
function settings = getsettings(this)
% Returns a structure containing all settings of this object.
% >> settings = tomography.getsettings();
diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
deleted file mode 100644
index b30640e..0000000
--- a/matlab/algorithms/DART/IterativeTomography.m
+++ /dev/null
@@ -1,455 +0,0 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
-%
-% Copyright: iMinds-Vision Lab, University of Antwerp
-% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
-
-classdef IterativeTomography < matlab.mixin.Copyable
-
- % Algorithm class for 2D Iterative Tomography.
-
- %----------------------------------------------------------------------
- properties (SetAccess=public, GetAccess=public)
- superresolution = 1; % SETTING: Volume upsampling factor.
- proj_type = 'linear'; % SETTING: Projector type, only when gpu='no'.
- method = 'SIRT_CUDA'; % SETTING: Iterative method (see ASTRA toolbox documentation).
- gpu = 'yes'; % SETTING: Use gpu? {'yes', 'no'}
- gpu_core = 0; % SETTING: Which gpu core to use? Only when gpu='yes'.
- inner_circle = 'yes'; % SETTING: Do roi only? {'yes', 'no'}
- image_size = []; % SETTING: Overwrite default reconstruction size. Only if no vol_geom is specified.
- use_minc = 'no'; % SETTING: Use minimum constraint. {'no', 'yes'}
- end
- %----------------------------------------------------------------------
- properties (SetAccess=public, GetAccess=public)
- sinogram = []; % DATA: Projection data.
- proj_geom = []; % DATA: Projection geometry.
- V = []; % DATA: Volume data. Also used to set initial estimate (optional).
- vol_geom = []; % DATA: Volume geometry.
- end
- %----------------------------------------------------------------------
- properties (SetAccess=private, GetAccess=public)
- initialized = 0; % Is this object initialized?
- end
- %----------------------------------------------------------------------
- properties (SetAccess=protected, GetAccess=protected)
- proj_geom_sr = []; % PROTECTED: geometry of sinogram (with super-resolution)
- proj_id = []; % PROTECTED: astra id of projector (when gpu='no')
- proj_id_sr = []; % PROTECTED: astra id of super-resolution projector (when gpu='no')
- cfg_base = struct(); % PROTECTED: base configuration structure for the reconstruction algorithm.
- end
- %----------------------------------------------------------------------
-
- methods (Access=public)
-
- %------------------------------------------------------------------
- function this = IterativeTomography(varargin)
- % Constructor
- % >> tomography = IterativeTomography();
- % >> tomography = IterativeTomography(base); base struct should contain fields 'sinogram' and 'proj_geom'.
- % >> tomography = IterativeTomography(base_filename); mat-file should contain 'sinogram' and 'proj_geom'.
- % >> tomography = IterativeTomography(proj_geom, vol_geom);
- % >> tomography = IterativeTomography(sinogram, proj_geom);
-
- % Input: IterativeTomography(base)
- if nargin == 1
- if ischar(varargin{1})
- this.sinogram = load(varargin{1},'sinogram');
- this.proj_geom = load(varargin{1},'proj_geom');
- else
- this.sinogram = varargin{1}.sinogram;
- this.proj_geom = varargin{1}.proj_geom;
- end
- end
- % Input: IterativeTomography(proj_geom, vol_geom)
- if nargin == 2 && isstruct(varargin{1}) && isstruct(varargin{2})
- this.proj_geom = varargin{1};
- this.vol_geom = varargin{2};
- % Input: IterativeTomography(sinogram, proj_geom)
- elseif nargin == 2
- this.sinogram = varargin{1};
- this.proj_geom = varargin{2};
- end
-
- end
-
- %------------------------------------------------------------------
- function delete(this)
- % Destructor
- % >> clear tomography;
- if strcmp(this.gpu,'no') && numel(this.proj_id) > 0
- astra_mex_projector('delete', this.proj_id, this.proj_id_sr);
- end
- end
-
- %------------------------------------------------------------------
- function settings = getsettings(this)
- % Returns a structure containing all settings of this object.
- % >> settings = tomography.getsettings();
- settings.superresolution = this.superresolution;
- settings.proj_type = this.proj_type;
- settings.method = this.method;
- settings.gpu = this.gpu;
- settings.gpu_core = this.gpu_core;
- settings.inner_circle = this.inner_circle;
- settings.image_size = this.image_size;
- settings.use_minc = this.use_minc;
- end
-
- %------------------------------------------------------------------
- function ok = initialize(this)
- % Initialize this object. Returns 1 if succesful.
- % >> tomography.initialize();
-
- % create projection geometry with super-resolution
- this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);
-
- % if no volume geometry is specified by the user: create volume geometry
- if numel(this.vol_geom) == 0
- if numel(this.image_size) < 2
- this.image_size(1) = this.proj_geom.DetectorCount;
- this.image_size(2) = this.proj_geom.DetectorCount;
- end
- this.vol_geom = astra_create_vol_geom(this.image_size(1) * this.superresolution, this.image_size(2) * this.superresolution, ...
- -this.image_size(1)/2, this.image_size(1)/2, -this.image_size(2)/2, this.image_size(2)/2);
- else
- this.image_size(1) = this.vol_geom.GridRowCount;
- this.image_size(2) = this.vol_geom.GridColCount;
- end
-
- % create projector
- if strcmp(this.gpu, 'no')
- this.proj_id = astra_create_projector(this.proj_type, this.proj_geom, this.vol_geom);
- this.proj_id_sr = astra_create_projector(this.proj_type, this.proj_geom_sr, this.vol_geom);
- end
-
- % create reconstruction configuration
- this.cfg_base = astra_struct(upper(this.method));
- if strcmp(this.gpu,'no')
- this.cfg_base.ProjectorId = this.proj_id;
- this.cfg_base.ProjectionGeometry = this.proj_geom;
- this.cfg_base.ReconstructionGeometry = this.vol_geom;
- this.cfg_base.option.ProjectionOrder = 'random';
- end
- this.cfg_base.option.DetectorSuperSampling = this.superresolution;
- %this.cfg_base.option.DetectorSuperSampling = 8;
- if strcmp(this.gpu,'yes')
- this.cfg_base.option.GPUindex = this.gpu_core;
- end
- this.cfg_base.option.UseMinConstraint = this.use_minc;
-
- this.initialized = 1;
- ok = this.initialized;
- end
-
- %------------------------------------------------------------------
- function P = project(this, volume)
- % Compute forward projection.
- % Stores sinogram in tomography.sinogram if it is still empty.
- % >> P = tomography.project(); projects 'tomography.V'.
- % >> P = tomography.project(volume); projects 'volume'.
-
- if ~this.initialized
- this.initialize();
- end
-
- if nargin == 1 % tomography.project();
- P = this.project_c(this.V);
-
- elseif nargin == 2 % tomography.project(volume);
- P = this.project_c(volume);
- end
-
- % store
- if numel(this.sinogram) == 0
- this.sinogram = P;
- end
- end
-
- %------------------------------------------------------------------
- function V = reconstruct(this, iterations)
- % Compute reconstruction.
- % Uses tomography.sinogram
- % Initial solution (if available) should be stored in tomography.V
- % >> V = tomography.reconstruct(iterations);
-
- if ~this.initialized
- this.initialize();
- end
-
- this.V = this.reconstruct_c(this.sinogram, this.V, [], iterations);
- if strcmp(this.inner_circle,'yes')
- this.V = this.selectROI(this.V);
- end
- V = this.V;
- end
-
- %------------------------------------------------------------------
- function I = reconstructMask(this, mask, iterations)
- % Compute reconstruction with mask.
- % Uses tomography.sinogram
- % Initial solution (if available) should be stored in tomography.V
- % >> V = tomography.reconstructMask(mask, iterations);
-
- if ~this.initialized
- this.initialize();
- end
-
- if strcmp(this.inner_circle,'yes')
- mask = this.selectROI(mask);
- end
- I = this.reconstruct_c(this.sinogram, this.V, mask, iterations);
- end
- %------------------------------------------------------------------
-
- end
-
- %----------------------------------------------------------------------
- methods (Access = protected)
-
- %------------------------------------------------------------------
- % Protected function: create FP
- function sinogram = project_c(this, volume)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- % data is stored in astra memory
- if numel(volume) == 1
-
- if strcmp(this.gpu, 'yes')
- sinogram_tmp = astra_create_sino_cuda(volume, this.proj_geom_sr, this.vol_geom, this.gpu_core);
- else
- sinogram_tmp = astra_create_sino(volume, this.proj_id);
- end
-
- % sinogram downsampling
- if this.superresolution > 1
- sinogram_data = astra_mex_data2d('get', sinogram_tmp);
- astra_mex_data2d('delete', sinogram_tmp);
- sinogram_data = downsample_sinogram(sinogram_data, this.superresolution);
- %if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- % sinogram = sinogram / this.superresolution;
- %end
- sinogram = astra_mex_data2d('create','sino', this.proj_geom, sinogram_data);
- else
- sinogram = sinogram_tmp;
- end
-
- % data is stored in matlab memory
- else
-
- % 2D and 3D slice by slice
- sinogram_tmp = zeros([astra_geom_size(this.proj_geom_sr), size(volume,3)]);
- sinogram_tmp2 = zeros([astra_geom_size(this.proj_geom), size(volume,3)]);
- for slice = 1:size(volume,3)
-
- if strcmp(this.gpu, 'yes')
- [tmp_id, sinogram_tmp2(:,:,slice)] = astra_create_sino_sampling(volume(:,:,slice), this.proj_geom, this.vol_geom, this.gpu_core, this.superresolution);
- astra_mex_data2d('delete', tmp_id);
- else
- [tmp_id, tmp] = astra_create_sino(volume(:,:,slice), this.proj_id_sr);
- sinogram_tmp2(:,:,slice) = downsample_sinogram(tmp, this.superresolution) * (this.superresolution^2);
- astra_mex_data2d('delete', tmp_id);
- end
-
- end
-
- % sinogram downsampling
- if strcmp(this.gpu, 'yes')
- %sinogram = downsample_sinogram(sinogram_tmp, this.superresolution);
- sinogram2 = sinogram_tmp2;
- if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- sinogram2 = sinogram2 / this.superresolution;
- elseif strcmp(this.proj_geom.type,'parallel')
- sinogram2 = sinogram2 / (this.superresolution * this.superresolution);
- end
- sinogram = sinogram2;
- else
- sinogram = sinogram_tmp2;
- end
-
- end
-
- end
-
- %------------------------------------------------------------------
- % Protected function: reconstruct
- function V = reconstruct_c(this, sinogram, V0, mask, iterations)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- % data is stored in astra memory
- if numel(sinogram) == 1
- V = this.reconstruct_c_astra(sinogram, V0, mask, iterations);
-
- % data is stored in matlab memory
- else
- V = this.reconstruct_c_matlab(sinogram, V0, mask, iterations);
- end
-
- end
-
- %------------------------------------------------------------------
- % Protected function: reconstruct (data in matlab)
- function V = reconstruct_c_matlab(this, sinogram, V0, mask, iterations)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- % parse method
- method2 = upper(this.method);
- if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
- iterations = iterations * size(sinogram,1);
- elseif strcmp(method2, 'ART')
- iterations = iterations * numel(sinogram);
- end
-
- % create data objects
- V = zeros(this.vol_geom.GridRowCount, this.vol_geom.GridColCount, size(sinogram,3));
- reconstruction_id = astra_mex_data2d('create', '-vol', this.vol_geom);
- sinogram_id = astra_mex_data2d('create', '-sino', this.proj_geom);
- if numel(mask) > 0
- mask_id = astra_mex_data2d('create', '-vol', this.vol_geom);
- end
-
- % algorithm configuration
- cfg = this.cfg_base;
- cfg.ProjectionDataId = sinogram_id;
- cfg.ReconstructionDataId = reconstruction_id;
- if numel(mask) > 0
- cfg.option.ReconstructionMaskId = mask_id;
- end
- alg_id = astra_mex_algorithm('create', cfg);
-
- % loop slices
- for slice = 1:size(sinogram,3)
-
- % fetch slice of initial reconstruction
- if numel(V0) > 0
- astra_mex_data2d('store', reconstruction_id, V0(:,:,slice));
- else
- astra_mex_data2d('store', reconstruction_id, 0);
- end
-
- % fetch slice of sinogram
- astra_mex_data2d('store', sinogram_id, sinogram(:,:,slice));
-
- % fecth slice of mask
- if numel(mask) > 0
- astra_mex_data2d('store', mask_id, mask(:,:,slice));
- end
-
- % iterate
- astra_mex_algorithm('iterate', alg_id, iterations);
-
- % fetch data
- V(:,:,slice) = astra_mex_data2d('get', reconstruction_id);
-
- end
-
- % correct attenuation factors for super-resolution
- if this.superresolution > 1 && strcmp(this.gpu,'yes')
- if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- if numel(mask) > 0
- V(mask > 0) = V(mask > 0) ./ this.superresolution;
- else
- V = V ./ this.superresolution;
- end
- end
- end
-
- % garbage collection
- astra_mex_algorithm('delete', alg_id);
- astra_mex_data2d('delete', sinogram_id, reconstruction_id);
- if numel(mask) > 0
- astra_mex_data2d('delete', mask_id);
- end
-
- end
-
- %------------------------------------------------------------------
- % Protected function: reconstruct (data in astra)
- function V = reconstruct_c_astra(this, sinogram, V0, mask, iterations)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- if numel(V0) > 1 || numel(mask) > 1 || numel(sinogram) > 1
- error('Not all required data is stored in the astra memory');
- end
-
- if numel(V0) == 0
- V0 = astra_mex_data2d('create', '-vol', this.vol_geom, 0);
- end
-
- % parse method
- method2 = upper(this.method);
- if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
- iterations = iterations * astra_geom_size(this.proj_geom, 1);
- elseif strcmp(method2, 'ART')
- s = astra_geom_size(this.proj_geom);
- iterations = iterations * s(1) * s(2);
- end
-
- % algorithm configuration
- cfg = this.cfg_base;
- cfg.ProjectionDataId = sinogram;
- cfg.ReconstructionDataId = V0;
- if numel(mask) > 0
- cfg.option.ReconstructionMaskId = mask;
- end
- alg_id = astra_mex_algorithm('create', cfg);
-
- % iterate
- astra_mex_algorithm('iterate', alg_id, iterations);
-
- % fetch data
- V = V0;
-
- % correct attenuation factors for super-resolution
- if this.superresolution > 1
- if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- if numel(mask) > 0
- astra_data_op_masked('$1./s1', [V V], [this.superresolution this.superresolution], mask, this.gpu_core);
- else
- astra_data_op('$1./s1', [V V], [this.superresolution this.superresolution], this.gpu_core);
- end
- end
- end
-
- % garbage collection
- astra_mex_algorithm('delete', alg_id);
-
- end
-
- %------------------------------------------------------------------
- function V_out = selectROI(~, V_in)
-
- if numel(V_in) == 1
- cfg = astra_struct('RoiSelect_CUDA');
- cfg.DataId = V_in;
- alg_id = astra_mex_algorithm('create',cfg);
- astra_mex_algorithm('run', alg_id);
- astra_mex_algorithm('delete', alg_id);
- V_out = V_in;
- else
- V_out = ROIselectfull(V_in, min([size(V_in,1), size(V_in,2)]));
- end
-
- end
- %------------------------------------------------------------------
-
- end
-
-end
-
diff --git a/matlab/algorithms/DART/IterativeTomography3D.m b/matlab/algorithms/DART/IterativeTomography3D.m
deleted file mode 100644
index 3f89457..0000000
--- a/matlab/algorithms/DART/IterativeTomography3D.m
+++ /dev/null
@@ -1,433 +0,0 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
-%
-% Copyright: iMinds-Vision Lab, University of Antwerp
-% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
-
-classdef IterativeTomography3D < matlab.mixin.Copyable
-
- % Algorithm class for 3D Iterative Tomography.
-
- %----------------------------------------------------------------------
- properties (SetAccess=public, GetAccess=public)
- superresolution = 1; % SETTING: Volume upsampling factor.
- proj_type = 'linear'; % SETTING: Projector type, only when gpu='no'.
- method = 'SIRT3D_CUDA'; % SETTING: Iterative method (see ASTRA toolbox documentation).
- gpu = 'yes'; % SETTING: Use gpu? {'yes', 'no'}
- gpu_core = 0; % SETTING: Which gpu core to use? Only when gpu='yes'.
- inner_circle = 'yes'; % SETTING: Do roi only? {'yes', 'no'}
- image_size = []; % SETTING: Overwrite default reconstruction size. Only if no vol_geom is specified.
- use_minc = 'no'; % SETTING: Use minimum constraint. {'no', 'yes'}
- maxc = +Inf; % SETTING: Maximum constraint. +Inf means off.
- end
- %----------------------------------------------------------------------
- properties (SetAccess=public, GetAccess=public)
- sinogram = []; % DATA: Projection data.
- proj_geom = []; % DATA: Projection geometry.
- V = []; % DATA: Volume data. Also used to set initial estimate (optional).
- vol_geom = []; % DATA: Volume geometry.
- end
- %----------------------------------------------------------------------
- properties (SetAccess=private, GetAccess=public)
- initialized = 0; % Is this object initialized?
- end
- %----------------------------------------------------------------------
- properties (SetAccess=protected, GetAccess=protected)
- proj_geom_sr = []; % PROTECTED: geometry of sinogram (with super-resolution)
- proj_id = []; % PROTECTED: astra id of projector (when gpu='no')
- proj_id_sr = []; % PROTECTED: astra id of super-resolution projector (when gpu='no')
- cfg_base = struct(); % PROTECTED: base configuration structure for the reconstruction algorithm.
- end
- %----------------------------------------------------------------------
-
- methods (Access=public)
-
- %------------------------------------------------------------------
- function this = IterativeTomography(varargin)
- % Constructor
- % >> tomography = IterativeTomography();
- % >> tomography = IterativeTomography(base); base struct should contain fields 'sinogram' and 'proj_geom'.
- % >> tomography = IterativeTomography(base_filename); mat-file should contain 'sinogram' and 'proj_geom'.
- % >> tomography = IterativeTomography(proj_geom, vol_geom);
- % >> tomography = IterativeTomography(sinogram, proj_geom);
-
- % Input: IterativeTomography(base)
- if nargin == 1
- if ischar(varargin{1})
- this.sinogram = load(varargin{1},'sinogram');
- this.proj_geom = load(varargin{1},'proj_geom');
- else
- this.sinogram = varargin{1}.sinogram;
- this.proj_geom = varargin{1}.proj_geom;
- end
- end
- % Input: IterativeTomography(proj_geom, vol_geom)
- if nargin == 2 && isstruct(varargin{1}) && isstruct(varargin{2})
- this.proj_geom = varargin{1};
- this.vol_geom = varargin{2};
- % Input: IterativeTomography(sinogram, proj_geom)
- elseif nargin == 2
- this.sinogram = varargin{1};
- this.proj_geom = varargin{2};
- end
-
- end
-
- %------------------------------------------------------------------
- function delete(this)
- % Destructor
- % >> clear tomography;
- if strcmp(this.gpu,'no') && numel(this.proj_id) > 0
- astra_mex_projector('delete', this.proj_id, this.proj_id_sr);
- end
- end
-
- %------------------------------------------------------------------
- function settings = getsettings(this)
- % Returns a structure containing all settings of this object.
- % >> settings = tomography.getsettings();
- settings.superresolution = this.superresolution;
- settings.proj_type = this.proj_type;
- settings.method = this.method;
- settings.gpu = this.gpu;
- settings.gpu_core = this.gpu_core;
- settings.inner_circle = this.inner_circle;
- settings.image_size = this.image_size;
- settings.use_minc = this.use_minc;
- settings.maxc = this.maxc;
- end
-
- %------------------------------------------------------------------
- function ok = initialize(this)
- % Initialize this object. Returns 1 if succesful.
- % >> tomography.initialize();
-
-% % create projection geometry with super-resolution
-% this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);
-
- % if no volume geometry is specified by the user: create volume geometry
- if numel(this.vol_geom) == 0
- if numel(this.image_size) < 2
- this.image_size(1) = this.proj_geom.DetectorRowCount;
- this.image_size(2) = this.proj_geom.DetectorColCount;
- end
- this.vol_geom = astra_create_vol_geom(this.proj_geom.DetectorColCount, this.proj_geom.DetectorColCount, this.proj_geom.DetectorRowCount);
- else
- this.image_size(1) = this.vol_geom.GridRowCount;
- this.image_size(2) = this.vol_geom.GridColCount;
- end
-
- % create projector
- if strcmp(this.gpu, 'no')
- this.proj_id = astra_create_projector(this.proj_type, this.proj_geom, this.vol_geom);
- this.proj_id_sr = astra_create_projector(this.proj_type, this.proj_geom_sr, this.vol_geom);
- end
-
- % create reconstruction configuration
- this.cfg_base = astra_struct(upper(this.method));
- if strcmp(this.gpu,'no')
- this.cfg_base.ProjectorId = this.proj_id;
- this.cfg_base.ProjectionGeometry = this.proj_geom;
- this.cfg_base.ReconstructionGeometry = this.vol_geom;
- this.cfg_base.option.ProjectionOrder = 'random';
- end
- this.cfg_base.option.DetectorSuperSampling = this.superresolution;
- %this.cfg_base.option.DetectorSuperSampling = 8;
- if strcmp(this.gpu,'yes')
- this.cfg_base.option.GPUindex = this.gpu_core;
- end
- this.cfg_base.option.UseMinConstraint = this.use_minc;
- if ~isinf(this.maxc)
- this.cfg_base.option.UseMaxConstraint = 'yes';
- this.cfg_base.option.MaxConstraintValue = this.maxc;
- end
-
- this.initialized = 1;
- ok = this.initialized;
- end
-
- %------------------------------------------------------------------
- function P = project(this, volume)
- % Compute forward projection.
- % Stores sinogram in tomography.sinogram if it is still empty.
- % >> P = tomography.project(); projects 'tomography.V'.
- % >> P = tomography.project(volume); projects 'volume'.
-
- if ~this.initialized
- this.initialize();
- end
-
- if nargin == 1 % tomography.project();
- P = this.project_c(this.V);
-
- elseif nargin == 2 % tomography.project(volume);
- P = this.project_c(volume);
- end
-
- % store
- if numel(this.sinogram) == 0
- this.sinogram = P;
- end
- end
-
- %------------------------------------------------------------------
- function V = reconstruct(this, iterations)
- % Compute reconstruction.
- % Uses tomography.sinogram
- % Initial solution (if available) should be stored in tomography.V
- % >> V = tomography.reconstruct(iterations);
-
- if ~this.initialized
- this.initialize();
- end
-
- this.V = this.reconstruct_c(this.sinogram, this.V, [], iterations);
- if strcmp(this.inner_circle,'yes')
- this.V = this.selectROI(this.V);
- end
- V = this.V;
- end
-
- %------------------------------------------------------------------
- function I = reconstructMask(this, mask, iterations)
- % Compute reconstruction with mask.
- % Uses tomography.sinogram
- % Initial solution (if available) should be stored in tomography.V
- % >> V = tomography.reconstructMask(mask, iterations);
-
- if ~this.initialized
- this.initialize();
- end
-
- if strcmp(this.inner_circle,'yes')
- mask = this.selectROI(mask);
- end
- I = this.reconstruct_c(this.sinogram, this.V, mask, iterations);
- end
- %------------------------------------------------------------------
-
- end
-
- %----------------------------------------------------------------------
- methods (Access = protected)
-
- %------------------------------------------------------------------
- % Protected function: create FP
- function sinogram = project_c(this, volume)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- % data is stored in astra memory
- if numel(volume) == 1
-
- if strcmp(this.gpu, 'yes')
- sinogram_tmp = astra_create_sino_cuda(volume, this.proj_geom_sr, this.vol_geom, this.gpu_core);
- else
- sinogram_tmp = astra_create_sino(volume, this.proj_id);
- end
-
- % sinogram downsampling
- if this.superresolution > 1
- sinogram_data = astra_mex_data2d('get', sinogram_tmp);
- astra_mex_data2d('delete', sinogram_tmp);
- sinogram_data = downsample_sinogram(sinogram_data, this.superresolution);
- %if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- % sinogram = sinogram / this.superresolution;
- %end
- sinogram = astra_mex_data2d('create','sino', this.proj_geom, sinogram_data);
- else
- sinogram = sinogram_tmp;
- end
-
- % data is stored in matlab memory
- else
-
- [tmp_id, sinogram] = astra_create_sino3d_cuda(volume, this.proj_geom, this.vol_geom);
- astra_mex_data3d('delete', tmp_id);
-
- end
-
- end
-
- %------------------------------------------------------------------
- % Protected function: reconstruct
- function V = reconstruct_c(this, sinogram, V0, mask, iterations)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- % data is stored in astra memory
- if numel(sinogram) == 1
- V = this.reconstruct_c_astra(sinogram, V0, mask, iterations);
-
- % data is stored in matlab memory
- else
- V = this.reconstruct_c_matlab(sinogram, V0, mask, iterations);
- end
-
- end
-
- %------------------------------------------------------------------
- % Protected function: reconstruct (data in matlab)
- function V = reconstruct_c_matlab(this, sinogram, V0, mask, iterations)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- % parse method
- method2 = upper(this.method);
- if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
- iterations = iterations * size(sinogram,1);
- elseif strcmp(method2, 'ART')
- iterations = iterations * numel(sinogram);
- end
-
- % create data objects
-% V = zeros(this.vol_geom.GridRowCount, this.vol_geom.GridColCount, size(sinogram,3));
- reconstruction_id = astra_mex_data3d('create', '-vol', this.vol_geom);
- sinogram_id = astra_mex_data3d('create', '-proj3d', this.proj_geom);
- if numel(mask) > 0
- mask_id = astra_mex_data3d('create', '-vol', this.vol_geom);
- end
-
- % algorithm configuration
- cfg = this.cfg_base;
- cfg.ProjectionDataId = sinogram_id;
- cfg.ReconstructionDataId = reconstruction_id;
- if numel(mask) > 0
- cfg.option.ReconstructionMaskId = mask_id;
- end
- alg_id = astra_mex_algorithm('create', cfg);
-
-% % loop slices
-% for slice = 1:size(sinogram,3)
-
- % fetch slice of initial reconstruction
- if numel(V0) > 0
- astra_mex_data3d('store', reconstruction_id, V0);
- else
- astra_mex_data3d('store', reconstruction_id, 0);
- end
-
- % fetch slice of sinogram
- astra_mex_data3d('store', sinogram_id, sinogram);
-
- % fecth slice of mask
- if numel(mask) > 0
- astra_mex_data3d('store', mask_id, mask);
- end
-
- % iterate
- astra_mex_algorithm('iterate', alg_id, iterations);
-
- % fetch data
- V = astra_mex_data3d('get', reconstruction_id);
-
-% end
-
- % correct attenuation factors for super-resolution
- if this.superresolution > 1 && strcmp(this.gpu,'yes')
- if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- if numel(mask) > 0
- V(mask > 0) = V(mask > 0) ./ this.superresolution;
- else
- V = V ./ this.superresolution;
- end
- end
- end
-
- % garbage collection
- astra_mex_algorithm('delete', alg_id);
- astra_mex_data3d('delete', sinogram_id, reconstruction_id);
- if numel(mask) > 0
- astra_mex_data3d('delete', mask_id);
- end
-
- end
-
- %------------------------------------------------------------------
- % Protected function: reconstruct (data in astra)
- function V = reconstruct_c_astra(this, sinogram, V0, mask, iterations)
-
- if this.initialized == 0
- error('IterativeTomography not initialized');
- end
-
- if numel(V0) > 1 || numel(mask) > 1 || numel(sinogram) > 1
- error('Not all required data is stored in the astra memory');
- end
-
- if numel(V0) == 0
- V0 = astra_mex_data2d('create', '-vol', this.vol_geom, 0);
- end
-
- % parse method
- method2 = upper(this.method);
- if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
- iterations = iterations * astra_geom_size(this.proj_geom, 1);
- elseif strcmp(method2, 'ART')
- s = astra_geom_size(this.proj_geom);
- iterations = iterations * s(1) * s(2);
- end
-
- % algorithm configuration
- cfg = this.cfg_base;
- cfg.ProjectionDataId = sinogram;
- cfg.ReconstructionDataId = V0;
- if numel(mask) > 0
- cfg.option.ReconstructionMaskId = mask;
- end
- alg_id = astra_mex_algorithm('create', cfg);
-
- % iterate
- astra_mex_algorithm('iterate', alg_id, iterations);
-
- % fetch data
- V = V0;
-
- % correct attenuation factors for super-resolution
- if this.superresolution > 1
- if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
- if numel(mask) > 0
- astra_data_op_masked('$1./s1', [V V], [this.superresolution this.superresolution], mask, this.gpu_core);
- else
- astra_data_op('$1./s1', [V V], [this.superresolution this.superresolution], this.gpu_core);
- end
- end
- end
-
- % garbage collection
- astra_mex_algorithm('delete', alg_id);
-
- end
-
- %------------------------------------------------------------------
- function V_out = selectROI(~, V_in)
-
- if numel(V_in) == 1
- cfg = astra_struct('RoiSelect_CUDA');
- cfg.DataId = V_in;
- alg_id = astra_mex_algorithm('create',cfg);
- astra_mex_algorithm('run', alg_id);
- astra_mex_algorithm('delete', alg_id);
- V_out = V_in;
- else
- V_out = ROIselectfull(V_in, min([size(V_in,1), size(V_in,2)]));
- end
-
- end
- %------------------------------------------------------------------
-
- end
-
-end
-
diff --git a/matlab/algorithms/DART/Kernels.m b/matlab/algorithms/DART/Kernels.m
index b5e3134..558f611 100644
--- a/matlab/algorithms/DART/Kernels.m
+++ b/matlab/algorithms/DART/Kernels.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef Kernels
%KERNELS Summary of this class goes here
diff --git a/matlab/algorithms/DART/MaskingDefault.m b/matlab/algorithms/DART/MaskingDefault.m
index 6bd25a5..f91a6f5 100644
--- a/matlab/algorithms/DART/MaskingDefault.m
+++ b/matlab/algorithms/DART/MaskingDefault.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef MaskingDefault < matlab.mixin.Copyable
@@ -179,7 +178,7 @@ classdef MaskingDefault < matlab.mixin.Copyable
%------------------------------------------------------------------
function Mask = apply_3D_gpu(this, S)
- vol_geom = astra_create_vol_geom(size(S));
+ vol_geom = astra_create_vol_geom(size(S,2), size(S,1), size(S,3));
data_id = astra_mex_data3d('create', '-vol', vol_geom, S);
cfg = astra_struct('DARTMASK3D_CUDA');
diff --git a/matlab/algorithms/DART/MaskingGPU.m b/matlab/algorithms/DART/MaskingGPU.m
index c4ef2b7..7d31aa8 100644
--- a/matlab/algorithms/DART/MaskingGPU.m
+++ b/matlab/algorithms/DART/MaskingGPU.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef MaskingGPU < matlab.mixin.Copyable
diff --git a/matlab/algorithms/DART/OutputDefault.m b/matlab/algorithms/DART/OutputDefault.m
index a34a430..33e908b 100644
--- a/matlab/algorithms/DART/OutputDefault.m
+++ b/matlab/algorithms/DART/OutputDefault.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef OutputDefault < matlab.mixin.Copyable
diff --git a/matlab/algorithms/DART/SegmentationDefault.m b/matlab/algorithms/DART/SegmentationDefault.m
index c1d7d99..98265c0 100644
--- a/matlab/algorithms/DART/SegmentationDefault.m
+++ b/matlab/algorithms/DART/SegmentationDefault.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef SegmentationDefault < matlab.mixin.Copyable
diff --git a/matlab/algorithms/DART/SmoothingDefault.m b/matlab/algorithms/DART/SmoothingDefault.m
index 58a8baa..1974fa1 100644
--- a/matlab/algorithms/DART/SmoothingDefault.m
+++ b/matlab/algorithms/DART/SmoothingDefault.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef SmoothingDefault < matlab.mixin.Copyable
@@ -153,7 +152,7 @@ classdef SmoothingDefault < matlab.mixin.Copyable
%------------------------------------------------------------------
function V_out = apply_3D_gpu(this, V_in)
- vol_geom = astra_create_vol_geom(size(V_in));
+ vol_geom = astra_create_vol_geom(size(V_in,2),size(V_in,1),size(V_in,3));
data_id = astra_mex_data3d('create', '-vol', vol_geom, V_in);
cfg = astra_struct('DARTSMOOTHING3D_CUDA');
diff --git a/matlab/algorithms/DART/SmoothingGPU.m b/matlab/algorithms/DART/SmoothingGPU.m
index 857da37..1249e19 100644
--- a/matlab/algorithms/DART/SmoothingGPU.m
+++ b/matlab/algorithms/DART/SmoothingGPU.m
@@ -1,13 +1,12 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef SmoothingGPU < matlab.mixin.Copyable
diff --git a/matlab/algorithms/DART/StatisticsDefault.m b/matlab/algorithms/DART/StatisticsDefault.m
index 7822c5f..9d33256 100644
--- a/matlab/algorithms/DART/StatisticsDefault.m
+++ b/matlab/algorithms/DART/StatisticsDefault.m
@@ -1,22 +1,21 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
%
-% Copyright: iMinds-Vision Lab, University of Antwerp
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
classdef StatisticsDefault < matlab.mixin.Copyable
% Default policy class for statistics for DART.
properties (Access=public)
- pixel_error = 'yes'; % SETTING: Store pixel error? {'yes','no'}
- proj_diff = 'yes'; % SETTING: Store projection difference? {'yes','no'}
- timing = 'yes'; % SETTING: Store timings? {'yes','no'}
+ pixel_error = 'no'; % SETTING: Store pixel error? {'yes','no'}
+ proj_diff = 'no'; % SETTING: Store projection difference? {'yes','no'}
+ timing = 'no'; % SETTING: Store timings? {'yes','no'}
end
diff --git a/matlab/algorithms/DART/TomographyDefault.m b/matlab/algorithms/DART/TomographyDefault.m
deleted file mode 100644
index 4db3905..0000000
--- a/matlab/algorithms/DART/TomographyDefault.m
+++ /dev/null
@@ -1,73 +0,0 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
-%
-% Copyright: iMinds-Vision Lab, University of Antwerp
-% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
-
-classdef TomographyDefault < IterativeTomography
-
- % Policy class for tomography for DART.
-
- %----------------------------------------------------------------------
- properties (Access=public)
- t = 5; % SETTING: # ARMiterations, each DART iteration.
- t0 = 100; % SETTING: # ARM iterations at DART initialization.
- end
- %----------------------------------------------------------------------
-
- methods
-
- %------------------------------------------------------------------
- function settings = getsettings(this)
- % Returns a structure containing all settings of this object.
- % >> settings = DART.tomography.getsettings();
-% settings = getsettings@IterativeTomography();
- settings.t = this.t;
- settings.t0 = this.t0;
- end
-
- %------------------------------------------------------------------
- function initialize(this, DART)
- % Initializes this object.
- % >> DART.tomography.initialize();
- this.proj_geom = DART.base.proj_geom;
- this.initialize@IterativeTomography();
- end
-
- %------------------------------------------------------------------
- function P = createForwardProjection(this, ~, volume)
- % Compute forward projection.
- % >> DART.tomography.createForwardProjection(DART, volume);
- P = this.project_c(volume);
- end
-
- %------------------------------------------------------------------
- function I = createReconstruction(this, ~, sinogram, V0, mask)
- % Compute reconstruction (with mask).
- % >> DART.tomography.createReconstruction(DART, sinogram, V0, mask);
- if strcmp(this.inner_circle,'yes')
- mask = ROIselectfull(mask, size(mask,1));
- end
- I = this.reconstruct_c(sinogram, V0, mask, this.t);
- end
-
- %------------------------------------------------------------------
- function I = createInitialReconstruction(this, ~, sinogram)
- % Compute reconstruction (initial).
- % >> DART.tomography.createInitialReconstruction(DART, sinogram);
- I = this.reconstruct_c(sinogram, [], [], this.t0);
- if strcmp(this.inner_circle,'yes')
- I = ROIselectfull(I, size(I,1));
- end
- end
- %------------------------------------------------------------------
-
- end
-
-end
-
diff --git a/matlab/algorithms/DART/TomographyDefault3D.m b/matlab/algorithms/DART/TomographyDefault3D.m
deleted file mode 100644
index 2be1b17..0000000
--- a/matlab/algorithms/DART/TomographyDefault3D.m
+++ /dev/null
@@ -1,73 +0,0 @@
-% This file is part of the
-% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
-%
-% Copyright: iMinds-Vision Lab, University of Antwerp
-% License: Open Source under GPLv3
-% Contact: mailto:astra@ua.ac.be
-% Website: http://astra.ua.ac.be
-%
-% Author of this DART Algorithm: Wim van Aarle
-
-
-classdef TomographyDefault3D < IterativeTomography3D
-
- % Policy class for 3D tomography for DART.
-
- %----------------------------------------------------------------------
- properties (Access=public)
- t = 5; % SETTING: # ARMiterations, each DART iteration.
- t0 = 100; % SETTING: # ARM iterations at DART initialization.
- end
- %----------------------------------------------------------------------
-
- methods
-
- %------------------------------------------------------------------
- function settings = getsettings(this)
- % Returns a structure containing all settings of this object.
- % >> settings = DART.tomography.getsettings();
-% settings = getsettings@IterativeTomography();
- settings.t = this.t;
- settings.t0 = this.t0;
- end
-
- %------------------------------------------------------------------
- function initialize(this, DART)
- % Initializes this object.
- % >> DART.tomography.initialize();
- this.proj_geom = DART.base.proj_geom;
- this.initialize@IterativeTomography3D();
- end
-
- %------------------------------------------------------------------
- function P = createForwardProjection(this, ~, volume)
- % Compute forward projection.
- % >> DART.tomography.createForwardProjection(DART, volume);
- P = this.project_c(volume);
- end
-
- %------------------------------------------------------------------
- function I = createReconstruction(this, ~, sinogram, V0, mask)
- % Compute reconstruction (with mask).
- % >> DART.tomography.createReconstruction(DART, sinogram, V0, mask);
- if strcmp(this.inner_circle,'yes')
- mask = ROIselectfull(mask, size(mask,1));
- end
- I = this.reconstruct_c(sinogram, V0, mask, this.t);
- end
-
- %------------------------------------------------------------------
- function I = createInitialReconstruction(this, ~, sinogram)
- % Compute reconstruction (initial).
- % >> DART.tomography.createInitialReconstruction(DART, sinogram);
- I = this.reconstruct_c(sinogram, [], [], this.t0);
- if strcmp(this.inner_circle,'yes')
- I = ROIselectfull(I, size(I,1));
- end
- end
- %------------------------------------------------------------------
-
- end
-
-end
-
diff --git a/matlab/algorithms/DART/examples/example1.m b/matlab/algorithms/DART/examples/example1.m
index daa3ce8..9a836f8 100644
--- a/matlab/algorithms/DART/examples/example1.m
+++ b/matlab/algorithms/DART/examples/example1.m
@@ -1,62 +1,67 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
clear all;
-addpath('..');
+addpath('../');
-%
-% Example 1: parallel beam, three slices.
+%%
+% Example 1: 2D parallel beam, cuda
%
% Configuration
-proj_count = 20;
-slice_count = 3;
+proj_count = 20;
dart_iterations = 20;
-filename = 'cylinders.png';
-outdir = './';
-prefix = 'example1';
-rho = [0, 1];
-tau = 0.5;
-gpu_core = 0;
+filename = 'cylinders.png';
+outdir = './';
+prefix = 'example1';
+rho = [0, 255];
+tau = 128;
+gpu_core = 0;
-% Load phantom.
-I = double(imread(filename)) / 255;
+% Load phantom
+I = imreadgs(filename);
-% Create projection and volume geometries.
+% Create projection and volume geometries
det_count = size(I, 1);
-angles = linspace(0, pi - pi / proj_count, proj_count);
-proj_geom = astra_create_proj_geom('parallel3d', 1, 1, slice_count, det_count, angles);
+angles = linspace2(0, pi, proj_count);
+proj_geom = astra_create_proj_geom('parallel', 1, det_count, angles);
vol_geom = astra_create_vol_geom(det_count, det_count, 1);
% Create sinogram.
-[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
-astra_mex_data3d('delete', sinogram_id);
+[sinogram_id, sinogram] = astra_create_sino_cuda(I, proj_geom, vol_geom);
+astra_mex_data2d('delete', sinogram_id);
-%
+%%
% DART
%
-base.sinogram = sinogram;
-base.proj_geom = proj_geom;
+%base.sinogram = sinogram;
+%base.proj_geom = proj_geom;
-D = DARTalgorithm(base);
+D = DARTalgorithm(sinogram, proj_geom);
+D.t0 = 100;
+D.t = 10;
-D.tomography = TomographyDefault3D();
-D.tomography.t0 = 100;
-D.tomography.t = 10;
-D.tomography.method = 'SIRT3D_CUDA';
-D.tomography.gpu_core = gpu_core;
-D.tomography.use_minc = 'yes';
-% D.tomography.maxc = 0.003; % Not a sensible value, just for demonstration.
+D.tomography.method = 'SIRT_CUDA';
+D.tomography.gpu_core = gpu_core;
+D.tomography.use_minc = 'yes';
-D.segmentation.rho = rho;
-D.segmentation.tau = tau;
+D.segmentation.rho = rho;
+D.segmentation.tau = tau;
-D.smoothing.b = 0.1;
-D.smoothing.full3d = 'yes';
-D.smoothing.gpu_core = gpu_core;
+D.smoothing.b = 0.1;
+D.smoothing.gpu_core = gpu_core;
-D.masking.random = 0.1;
-D.masking.conn = 6;
-D.masking.gpu_core = gpu_core;
+D.masking.random = 0.1;
+D.masking.gpu_core = gpu_core;
D.output.directory = outdir;
D.output.pre = [prefix '_'];
@@ -69,11 +74,10 @@ D.statistics.proj_diff = 'no';
D.initialize();
-disp([D.output.directory D.output.pre]);
-
D.iterate(dart_iterations);
+%%
% Convert middle slice of final iteration to png.
-load([outdir '/' prefix '_results_' num2str(dart_iterations) '.mat']);
-imwritesc(D.S(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_S.png']);
-imwritesc(D.V(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_V.png']);
+%
+imwritesc(D.S, [outdir '/' prefix '_S.png']);
+imwritesc(D.V, [outdir '/' prefix '_V.png']);
diff --git a/matlab/algorithms/DART/examples/example2.m b/matlab/algorithms/DART/examples/example2.m
index 8ee8cba..7fe6988 100644
--- a/matlab/algorithms/DART/examples/example2.m
+++ b/matlab/algorithms/DART/examples/example2.m
@@ -1,52 +1,56 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
clear all;
-addpath('..');
+addpath('../');
-%
-% Example 2: cone beam, full cube.
+%%
+% Example Z: 3D parallel beam, cuda
%
% Configuration
-det_count = 128;
-proj_count = 45;
-slice_count = det_count;
+proj_count = 20;
dart_iterations = 20;
-outdir = './';
-prefix = 'example2';
-rho = [0 0.5 1];
-tau = [0.25 0.75];
-gpu_core = 0;
-
-% Create phantom.
-% I = phantom3d([1 0.9 0.9 0.9 0 0 0 0 0 0; -0.5 0.8 0.8 0.8 0 0 0 0 0 0; -0.5 0.3 0.3 0.3 0 0 0 0 0 0], det_count);
-% save('phantom3d', 'I');
-load('phantom3d'); % Loads I.
-
-% Create projection and volume geometries.
-angles = linspace(0, pi - pi / proj_count, proj_count);
-proj_geom = astra_create_proj_geom('cone', 1, 1, slice_count, det_count, angles, 500, 0);
-vol_geom = astra_create_vol_geom(det_count, det_count, slice_count);
-
-% Create sinogram.
+outdir = './';
+prefix = 'example2';
+rho = [0, 0.5, 1];
+tau = [0.25, 0.75];
+gpu_core = 0;
+
+% Load phantom
+load('phantom3d');
+
+% Create projection and volume geometries
+det_count = size(I, 1);
+slice_count = size(I,3);
+angles = linspace2(0, pi, proj_count);
+proj_geom = astra_create_proj_geom('parallel3d', 1, 1, slice_count, det_count, angles);
+vol_geom = astra_create_vol_geom(size(I));
+
+% Create sinogram
[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
astra_mex_data3d('delete', sinogram_id);
-%
+%%
% DART
%
-base.sinogram = sinogram;
-base.proj_geom = proj_geom;
+D = DARTalgorithm(sinogram, proj_geom);
+D.t0 = 100;
+D.t = 10;
-D = DARTalgorithm(base);
-
-D.tomography = TomographyDefault3D();
-D.tomography.t0 = 100;
-D.tomography.t = 10;
-D.tomography.method = 'SIRT3D_CUDA';
-D.tomography.gpu_core = gpu_core;
-D.tomography.use_minc = 'yes';
-% D.tomography.maxc = 0.003; % Not a sensible value, just for demonstration.
+D.tomography = IterativeTomography3D();
+D.tomography.method = 'SIRT3D_CUDA';
+D.tomography.gpu_core = gpu_core;
+D.tomography.use_minc = 'yes';
D.segmentation.rho = rho;
D.segmentation.tau = tau;
@@ -56,7 +60,7 @@ D.smoothing.full3d = 'yes';
D.smoothing.gpu_core = gpu_core;
D.masking.random = 0.1;
-D.masking.conn = 6;
+D.masking.conn = 4;
D.masking.gpu_core = gpu_core;
D.output.directory = outdir;
@@ -70,11 +74,10 @@ D.statistics.proj_diff = 'no';
D.initialize();
-disp([D.output.directory D.output.pre]);
-
D.iterate(dart_iterations);
-% Convert middle slice of final iteration to png.
-load([outdir '/' prefix '_results_' num2str(dart_iterations) '.mat']);
-imwritesc(D.S(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_S.png']);
-imwritesc(D.V(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_V.png']);
+%%
+% Convert output of final iteration to png.
+%
+imwritesc(D.S(:,:,64), [outdir '/' prefix '_S_slice_64.png']);
+imwritesc(D.V(:,:,64), [outdir '/' prefix '_V_slice_64.png']);
diff --git a/matlab/algorithms/DART/examples/example3.m b/matlab/algorithms/DART/examples/example3.m
index 1c04f86..895630b 100644
--- a/matlab/algorithms/DART/examples/example3.m
+++ b/matlab/algorithms/DART/examples/example3.m
@@ -1,51 +1,56 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
clear all;
-addpath('..');
+addpath('../');
-%
-% Example 3: parallel beam, 2D
+%%
+% Example 3: 3D cone beam, cuda
%
% Configuration
-proj_count = 20;
+proj_count = 20;
dart_iterations = 20;
-filename = 'cylinders.png';
-outdir = './';
-prefix = 'example3';
-rho = [0, 1];
-tau = 0.5;
-gpu_core = 0;
+outdir = './';
+prefix = 'example3';
+rho = [0, 0.5, 1];
+tau = [0.25, 0.75];
+gpu_core = 0;
-% Load phantom.
-I = double(imread(filename)) / 255;
+% Load phantom
+load('phantom3d');
-% Create projection and volume geometries.
+% Create projection and volume geometries
det_count = size(I, 1);
-angles = linspace(0, pi - pi / proj_count, proj_count);
-proj_geom = astra_create_proj_geom('parallel', 1, det_count, angles);
-vol_geom = astra_create_vol_geom(det_count, det_count);
+slice_count = size(I,3);
+angles = linspace2(0, pi, proj_count);
+proj_geom = astra_create_proj_geom('cone', 1, 1, slice_count, det_count, angles, 500, 0);
+vol_geom = astra_create_vol_geom(size(I));
-% Create sinogram.
-[sinogram_id, sinogram] = astra_create_sino_cuda(I, proj_geom, vol_geom);
-astra_mex_data2d('delete', sinogram_id);
+% Create sinogram
+[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
+astra_mex_data3d('delete', sinogram_id);
-%
+%%
% DART
%
-base.sinogram = sinogram;
-base.proj_geom = proj_geom;
+D = DARTalgorithm(sinogram, proj_geom);
+D.t0 = 100;
+D.t = 10;
-D = DARTalgorithm(base);
-
-%D.tomography = TomographyDefault3D();
-D.tomography.t0 = 100;
-D.tomography.t = 10;
-D.tomography.method = 'SIRT_CUDA';
-D.tomography.proj_type = 'strip';
-D.tomography.gpu_core = gpu_core;
-D.tomography.use_minc = 'yes';
-% D.tomography.maxc = 0.003; % Not a sensible value, just for demonstration.
+D.tomography = IterativeTomography3D();
+D.tomography.method = 'SIRT3D_CUDA';
+D.tomography.gpu_core = gpu_core;
+D.tomography.use_minc = 'yes';
D.segmentation.rho = rho;
D.segmentation.tau = tau;
@@ -69,11 +74,10 @@ D.statistics.proj_diff = 'no';
D.initialize();
-disp([D.output.directory D.output.pre]);
-
D.iterate(dart_iterations);
+%%
% Convert output of final iteration to png.
-load([outdir '/' prefix '_results_' num2str(dart_iterations) '.mat']);
-imwritesc(D.S, [outdir '/' prefix '_S.png']);
-imwritesc(D.V, [outdir '/' prefix '_V.png']);
+%
+imwritesc(D.S(:,:,64), [outdir '/' prefix '_S_slice_64.png']);
+imwritesc(D.V(:,:,64), [outdir '/' prefix '_V_slice_64.png']);
diff --git a/matlab/algorithms/DART/tools/DARToptimizer.m b/matlab/algorithms/DART/tools/DARToptimizer.m
new file mode 100644
index 0000000..7aeabbe
--- /dev/null
+++ b/matlab/algorithms/DART/tools/DARToptimizer.m
@@ -0,0 +1,118 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
+classdef DARToptimizer < handle
+
+ %----------------------------------------------------------------------
+ properties (SetAccess=public,GetAccess=public)
+
+ % optimization options
+ max_evals = 100; % SETTING: Maximum number of function evaluations during optimization.
+ tolerance = 0.1; % SETTING: Minimum tolerance to achieve.
+ display = 'off'; % SETTING: Optimization output. {'off','on','iter'}
+ verbose = 'yes'; % SETTING: verbose? {'yes','no}
+
+ metric = ProjDiffOptimFunc(); % SETTING: Optimization object. Default: ProjDiffOptimFunc.
+
+ % DART options
+ DART_iterations = 20; % SETTING: number of DART iterations in each evaluation.
+
+ D_base = [];
+
+ end
+
+ %----------------------------------------------------------------------
+ properties (SetAccess=private,GetAccess=public)
+
+ stats = Statistics();
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ % Constructor
+ function this = DARToptimizer(D_base)
+
+ this.D_base = D_base;
+
+ % statistics
+ this.stats = Statistics();
+ this.stats.register('params');
+ this.stats.register('values');
+ this.stats.register('score');
+ end
+
+ %------------------------------------------------------------------
+ function opt_values = run(this, params, initial_values)
+
+ if nargin < 3
+ for i = 1:numel(params)
+ initial_values(i) = eval(['this.D_base.' params{i} ';']);
+ end
+ end
+
+ % fminsearch
+ options = optimset('display', this.display, 'MaxFunEvals', this.max_evals, 'TolX', this.tolerance);
+ opt_values = fminsearch(@this.optim_func, initial_values, options, params);
+
+ % save to D_base
+ for i = 1:numel(params)
+ eval(sprintf('this.D_base.%s = %d;',params{i}, opt_values(i)));
+ end
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=protected)
+
+ %------------------------------------------------------------------
+ function score = optim_func(this, values, params)
+
+ % copy DART
+ D = this.D_base.deepcopy();
+
+ % set parameters
+ for i = 1:numel(params)
+ eval(sprintf('D.%s = %d;',params{i}, values(i)));
+ D.output.pre = [D.output.pre num2str(values(i)) '_'];
+ end
+
+ % evaluate
+ if D.initialized == 0
+ D.initialize();
+ end
+ rng('default');
+ D.iterate(this.DART_iterations);
+
+ % compute score
+ score = this.metric.calculate(D, this);
+
+ % statistics
+ this.stats.add('params',params);
+ this.stats.add('values',values);
+ this.stats.add('score',score);
+
+ % output
+ if strcmp(this.verbose,'yes')
+ disp([num2str(values) ': ' num2str(score)]);
+ end
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/tools/DARToptimizerBoneStudy.m b/matlab/algorithms/DART/tools/DARToptimizerBoneStudy.m
new file mode 100644
index 0000000..5bcf6e9
--- /dev/null
+++ b/matlab/algorithms/DART/tools/DARToptimizerBoneStudy.m
@@ -0,0 +1,118 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
+classdef DARToptimizerBoneStudy < handle
+
+ %----------------------------------------------------------------------
+ properties (SetAccess=public,GetAccess=public)
+
+ % optimization options
+ max_evals = 100;
+ tolerance = 0.1;
+ display = 'off';
+
+ % DART options
+ DART_iterations = 50;
+
+ D_base = [];
+
+ end
+
+ %----------------------------------------------------------------------
+ properties (SetAccess=private,GetAccess=public)
+
+ stats = struct();
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ % Constructor
+ function this = DARToptimizerBoneStudy(D_base)
+
+ this.D_base = D_base;
+
+ this.stats.params = {};
+ this.stats.values = [];
+ this.stats.rmse = [];
+ this.stats.f_250 = [];
+ this.stats.f_100 = [];
+ this.stats.w_250 = [];
+ this.stats.w_125 = [];
+
+ end
+
+ %------------------------------------------------------------------
+ function opt_values = run(this, params, initial_values)
+
+ if nargin < 3
+ for i = 1:numel(params)
+ initial_values(i) = eval(['this.D_base.' params{i} ';']);
+ end
+ end
+
+ % fminsearch
+ options = optimset('display', this.display, 'MaxFunEvals', this.max_evals, 'TolX', this.tolerance);
+ opt_values = fminsearch(@optim_func, initial_values, options, this.D_base, params, this);
+
+ % save to D_base
+ for i = 1:numel(params)
+ eval(sprintf('this.D_base.%s = %d;',params{i}, opt_values(i)));
+ end
+
+ end
+ %------------------------------------------------------------------
+ end
+
+end
+
+%--------------------------------------------------------------------------
+function rmse = optim_func(values, D_base, params, Optim)
+
+ % copy DART
+ D = D_base.deepcopy();
+
+ % set parameters
+ for i = 1:numel(params)
+ eval(sprintf('D.%s = %d;',params{i}, values(i)));
+ D.output.pre = [D.output.pre num2str(values(i)) '_'];
+ end
+
+ % evaluate
+ if D.initialized == 0
+ D.initialize();
+ end
+ rng('default');
+ D.iterate(Optim.DART_iterations);
+
+ % compute rmse
+ ROI = load('roi.mat');
+ [rmse, f_250, f_100, w_250, w_125] = compute_rmse(D.S, ROI);
+ %projection = D.tomography.createForwardProjection(D, D.S);
+ %proj_diff = sum((projection(:) - D.base.sinogram(:)).^2);
+
+ % save
+ Optim.stats.params{end+1} = params;
+ Optim.stats.values(end+1,:) = values;
+ Optim.stats.rmse(end+1) = rmse;
+ Optim.stats.f_250(end+1) = f_250;
+ Optim.stats.f_100(end+1) = f_100;
+ Optim.stats.w_250(end+1) = w_250;
+ Optim.stats.w_125(end+1) = w_125;
+
+ disp([num2str(values) ': ' num2str(rmse)]);
+
+end
+
+
+
+
diff --git a/matlab/algorithms/DART/tools/ProjDiffOptimFunc.m b/matlab/algorithms/DART/tools/ProjDiffOptimFunc.m
new file mode 100644
index 0000000..64f3091
--- /dev/null
+++ b/matlab/algorithms/DART/tools/ProjDiffOptimFunc.m
@@ -0,0 +1,29 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
+classdef ProjDiffOptimFunc < handle
+
+ %----------------------------------------------------------------------
+ properties (SetAccess=private, GetAccess=public)
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ function proj_diff = calculate(~, D, ~)
+ projection = D.tomography.createForwardProjection(D, D.S);
+ proj_diff = sum((projection(:) - D.base.sinogram(:)).^2);
+ end
+
+ end
+ %----------------------------------------------------------------------
+
+end
diff --git a/matlab/algorithms/DART/tools/dart_create_base_phantom.m b/matlab/algorithms/DART/tools/dart_create_base_phantom.m
new file mode 100644
index 0000000..4c7061e
--- /dev/null
+++ b/matlab/algorithms/DART/tools/dart_create_base_phantom.m
@@ -0,0 +1,17 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
+function base = dart_create_base_phantom(Im, angle_count, angle_range, gpu_core)
+
+base.phantom = imreadgs(Im);
+base.proj_geom = astra_create_proj_geom('parallel', 1, size(base.phantom,1), linspace2(angle_range(1),angle_range(2),angle_count));
+base.vol_geom = astra_create_vol_geom(size(base.phantom,1),size(base.phantom,2));
+[sino_id, base.sinogram] = astra_create_sino_cuda(base.phantom, base.proj_geom, base.vol_geom, gpu_core);
+astra_mex_data2d('delete',sino_id);
diff --git a/matlab/algorithms/DART/tools/dart_scheduler.m b/matlab/algorithms/DART/tools/dart_scheduler.m
new file mode 100644
index 0000000..007fcaa
--- /dev/null
+++ b/matlab/algorithms/DART/tools/dart_scheduler.m
@@ -0,0 +1,53 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
+function dart_scheduler(D_tmpl, iterations, settings)
+
+
+for base_index = 1:numel(settings.base_loop)
+
+ % create new DART object
+ D = DART(settings.base_loop{base_index});
+
+ dart_scheduler1D(D, iterations, settings.parameter1);
+
+ % copy from templates
+ D.tomography = D_tmpl.tomography;
+ D.smoothing = D_tmpl.smoothing;
+ D.segmentation = D_tmpl.segmentation;
+ D.masking = D_tmpl.masking;
+ D.statistics = D_tmpl.statistics;
+ D.output = D_tmpl.output;
+
+ % set output options
+ D.output = OutputScheduler();
+ D.output.directory = output_folder{base_index};
+
+ % run DART
+ D = D.initialize();
+ D = D.iterate(iterations);
+
+end
+
+end
+
+function dart_scheduler1d(D, iterations, parameter_loop1)
+
+end
+
+
+
+
+
+
+
+
+
+
diff --git a/matlab/algorithms/DART/tools/rNMPOptimFunc.m b/matlab/algorithms/DART/tools/rNMPOptimFunc.m
new file mode 100644
index 0000000..7649d7f
--- /dev/null
+++ b/matlab/algorithms/DART/tools/rNMPOptimFunc.m
@@ -0,0 +1,35 @@
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2014, iMinds-Vision Lab, University of Antwerp
+% 2014, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+
+classdef rNMPOptimFunc < handle
+
+ %----------------------------------------------------------------------
+ properties (SetAccess=private, GetAccess=public)
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ function rnmp = calculate(~, D, ~)
+ if isfield(D.stats,'rnmp');
+ rnmp = D.stats.rnmp;
+ else
+ rnmp = compute_rnmp(D.base.phantom, D.S);
+ end
+ end
+
+ end
+ %----------------------------------------------------------------------
+
+
+end
+
+