summaryrefslogtreecommitdiffstats
path: root/matlab/algorithms/DART
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <WillemJan.Palenstijn@uantwerpen.be>2013-07-01 22:34:11 +0000
committerwpalenst <WillemJan.Palenstijn@uantwerpen.be>2013-07-01 22:34:11 +0000
commitb2fc6c70434674d74551c3a6c01ffb3233499312 (patch)
treeb17f080ebc504ab85ebb7c3d89f917fd87ce9e00 /matlab/algorithms/DART
downloadastra-b2fc6c70434674d74551c3a6c01ffb3233499312.tar.gz
astra-b2fc6c70434674d74551c3a6c01ffb3233499312.tar.bz2
astra-b2fc6c70434674d74551c3a6c01ffb3233499312.tar.xz
astra-b2fc6c70434674d74551c3a6c01ffb3233499312.zip
Update version to 1.3
Diffstat (limited to 'matlab/algorithms/DART')
-rw-r--r--matlab/algorithms/DART/DARTalgorithm.m229
-rw-r--r--matlab/algorithms/DART/IterativeTomography.m455
-rw-r--r--matlab/algorithms/DART/IterativeTomography3D.m433
-rw-r--r--matlab/algorithms/DART/Kernels.m68
-rw-r--r--matlab/algorithms/DART/MaskingDefault.m212
-rw-r--r--matlab/algorithms/DART/MaskingGPU.m94
-rw-r--r--matlab/algorithms/DART/OutputDefault.m173
-rw-r--r--matlab/algorithms/DART/SegmentationDefault.m55
-rw-r--r--matlab/algorithms/DART/SmoothingDefault.m179
-rw-r--r--matlab/algorithms/DART/SmoothingGPU.m119
-rw-r--r--matlab/algorithms/DART/StatisticsDefault.m72
-rw-r--r--matlab/algorithms/DART/TomographyDefault.m73
-rw-r--r--matlab/algorithms/DART/TomographyDefault3D.m73
-rw-r--r--matlab/algorithms/DART/examples/cylinders.pngbin0 -> 3934 bytes
-rw-r--r--matlab/algorithms/DART/examples/example1.m79
-rw-r--r--matlab/algorithms/DART/examples/example2.m80
-rw-r--r--matlab/algorithms/DART/examples/example3.m79
-rw-r--r--matlab/algorithms/DART/examples/phantom3d.matbin0 -> 242654 bytes
18 files changed, 2473 insertions, 0 deletions
diff --git a/matlab/algorithms/DART/DARTalgorithm.m b/matlab/algorithms/DART/DARTalgorithm.m
new file mode 100644
index 0000000..b7e63b5
--- /dev/null
+++ b/matlab/algorithms/DART/DARTalgorithm.m
@@ -0,0 +1,229 @@
+% 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 DARTalgorithm < matlab.mixin.Copyable
+
+ % Algorithm class for Discrete Algebraic Reconstruction Technique (DART).
+
+ % todo: reset()
+ % todo: fixed random seed
+ % todo: initialize from settings (?)
+
+ %----------------------------------------------------------------------
+ properties (GetAccess=public, SetAccess=public)
+
+ tomography = TomographyDefault(); % POLICY: Tomography object.
+ segmentation = SegmentationDefault(); % POLICY: Segmentation object.
+ smoothing = SmoothingDefault(); % POLICY: Smoothing object.
+ masking = MaskingDefault(); % POLICY: Masking object.
+ output = OutputDefault(); % POLICY: Output object.
+ statistics = StatisticsDefault(); % POLICY: Statistics object.
+
+ base = struct(); % DATA(set): base structure, should contain: 'sinogram', 'proj_geom', 'phantom' (optional).
+
+ memory = 'yes'; % SETTING: reduce memory usage? (disables some features)
+
+ testdata = struct();
+
+ end
+ %----------------------------------------------------------------------
+ properties (GetAccess=public, SetAccess=private)
+ V0 = []; % DATA(get): Initial reconstruction.
+ V = []; % DATA(get): Reconstruction.
+ S = []; % DATA(get): Segmentation.
+ R = []; % DATA(get): Residual projection data.
+ Mask = []; % DATA(get): Reconstruction Mask.
+ stats = struct(); % Structure containing various statistics.
+ iterationcount = 0; % Number of performed iterations.
+ start_tic = 0;
+ initialized = 0; % Is initialized?
+ end
+ %----------------------------------------------------------------------
+ properties (Access=private)
+ adaptparam_name = {};
+ adaptparam_values = {};
+ adaptparam_iters = {};
+ end
+
+ %----------------------------------------------------------------------
+ methods
+
+ %------------------------------------------------------------------
+ function this = DARTalgorithm(base)
+ % 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);
+ else
+ this.base = base;
+ end
+ end
+
+
+ %------------------------------------------------------------------
+ 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)
+ if isa(this.(props{i}), 'handle')
+ D.(props{i}) = copy(this.(props{i}));
+ end
+ end
+
+ end
+
+ %------------------------------------------------------------------
+ function this = initialize(this)
+ % Initializes this object.
+ % >> D.initialize();
+
+ % Initialize tomography part
+ this.tomography.initialize(this);
+
+ % 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.output.post_initial_iteration(this);
+ end
+ this.V = this.V0;
+ if strcmp(this.memory,'yes')
+ this.base.V0 = [];
+ this.V0 = [];
+ end
+ this.initialized = 1;
+ end
+
+ %------------------------------------------------------------------
+ % iterate
+ function this = iterate(this, iters)
+ % Perform several iterations of the DART algorithm.
+ % >> D.iterate(iterations);
+
+ this.start_tic = tic;
+
+ for iteration = 1:iters
+
+ this.iterationcount = this.iterationcount + 1;
+
+ %----------------------------------------------------------
+ % Initial 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
+
+ %----------------------------------------------------------
+ % 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 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
+ this.V = this.smoothing.apply(this, this.V);
+ %this.testdata.V3{iteration} = this.V;
+
+ %----------------------------------------------------------
+ % Calculate Statistics
+ this.stats = this.statistics.apply(this);
+
+ %----------------------------------------------------------
+ % Output
+ this.output.post_iteration(this);
+
+ end % end iteration loop
+
+ %test = this.testdata;
+ %save('testdata.mat','test');
+
+ end
+
+ %------------------------------------------------------------------
+ % get data
+ function data = getdata(this, string)
+ if numel(this.(string)) == 1
+ data = astra_mex_data2d('get',this.(string));
+ else
+ data = this.(string);
+ end
+ end
+
+ %------------------------------------------------------------------
+ % add adaptive parameter
+ function this = adaptiveparameter(this, name, values, iterations)
+ this.adaptparam_name{end+1} = name;
+ this.adaptparam_values{end+1} = values;
+ this.adaptparam_iters{end+1} = iterations;
+ end
+
+ %------------------------------------------------------------------
+ function settings = getsettings(this)
+ % Returns a structure containing all settings of this object.
+ % >> settings = tomography.getsettings();
+
+ settings.tomography = this.tomography.getsettings();
+ settings.smoothing = this.smoothing.getsettings();
+ settings.masking = this.masking.getsettings();
+ settings.segmentation = this.segmentation.getsettings();
+ end
+ %------------------------------------------------------------------
+
+ end % methods
+
+end % class
+
+
diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
new file mode 100644
index 0000000..b30640e
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography.m
@@ -0,0 +1,455 @@
+% 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
new file mode 100644
index 0000000..3f89457
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography3D.m
@@ -0,0 +1,433 @@
+% 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
new file mode 100644
index 0000000..b5e3134
--- /dev/null
+++ b/matlab/algorithms/DART/Kernels.m
@@ -0,0 +1,68 @@
+% 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 Kernels
+ %KERNELS Summary of this class goes here
+ % Detailed explanation goes here
+
+ properties
+
+ end
+
+ methods(Static)
+
+ function K = BinaryPixelKernel(radius, conn)
+
+ if nargin < 2
+ conn = 8;
+ end
+
+ % 2D, 4conn
+ if conn == 4
+ K = [0 1 0; 1 1 1; 0 1 0];
+ for i = 2:radius
+ K = conv2(K,K);
+ end
+ K = double(K >= 1);
+
+ % 2D, 8conn
+ elseif conn == 8
+ K = ones(2*radius+1, 2*radius+1);
+
+ % 3D, 6conn
+ elseif conn == 6
+ K = zeros(3,3,3);
+ K(:,:,1) = [0 0 0; 0 1 0; 0 0 0];
+ K(:,:,2) = [0 1 0; 1 1 1; 0 1 0];
+ K(:,:,3) = [0 0 0; 0 1 0; 0 0 0];
+ for i = 2:radius
+ K = convn(K,K);
+ end
+ K = double(K >= 1);
+
+ % 2D, 27conn
+ elseif conn == 26
+ K = ones(2*radius+1, 2*radius+1, 2*radius+1);
+
+ else
+ disp('Invalid conn')
+ end
+ end
+
+
+
+
+
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/MaskingDefault.m b/matlab/algorithms/DART/MaskingDefault.m
new file mode 100644
index 0000000..6bd25a5
--- /dev/null
+++ b/matlab/algorithms/DART/MaskingDefault.m
@@ -0,0 +1,212 @@
+% 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 MaskingDefault < matlab.mixin.Copyable
+
+ % Default policy class for masking for DART.
+
+ %----------------------------------------------------------------------
+ properties (Access=public)
+
+ radius = 1; % SETTING: Radius of masking kernel.
+ conn = 8; % SETTING: Connectivity window. For 2D: 4 or 8. For 3D: 6 or 26.
+ edge_threshold = 1; % SETTING: Number of pixels in the window that should be different.
+ random = 0.1; % SETTING: Percentage of random points. Between 0 and 1.
+ gpu = 'yes'; % SETTING: Use gpu? {'yes', 'no'}
+ gpu_core = 0; % SETTING: Which gpu core to use, only when gpu='yes'.
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function settings = getsettings(this)
+ % Returns a structure containing all settings of this object.
+ % >> settings = DART.masking.getsettings();
+ settings.radius = this.radius;
+ settings.conn = this.conn;
+ settings.edge_threshold = this.edge_threshold;
+ settings.random = this.random;
+ end
+
+ %------------------------------------------------------------------
+ function Mask = apply(this, ~, S)
+ % Applies masking.
+ % >> Mask = DART.segmentation.apply(DART, S);
+
+ % 2D, one slice
+ if size(S,3) == 1
+ if strcmp(this.gpu,'yes')
+ Mask = this.apply_2D_gpu(S);
+ else
+ Mask = this.apply_2D(S);
+ end
+
+ % 3D, slice by slice
+ elseif this.conn == 4 || this.conn == 8
+ Mask = zeros(size(S));
+ for slice = 1:size(S,3)
+ if strcmp(this.gpu,'yes')
+ Mask(:,:,slice) = this.apply_2D_gpu(S(:,:,slice));
+ else
+ Mask(:,:,slice) = this.apply_2D(S(:,:,slice));
+ end
+ end
+
+ % 3D, full
+ else
+ if strcmp(this.gpu,'yes')
+ Mask = this.apply_3D_gpu(S);
+ else
+ Mask = this.apply_3D(S);
+ end
+ end
+
+ end
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=protected)
+
+ %------------------------------------------------------------------
+ function Mask = apply_2D_gpu(this, S)
+
+ vol_geom = astra_create_vol_geom(size(S));
+ data_id = astra_mex_data2d('create', '-vol', vol_geom, S);
+ mask_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+ cfg = astra_struct('DARTMASK_CUDA');
+ cfg.SegmentationDataId = data_id;
+ cfg.MaskDataId = mask_id;
+ cfg.option.GPUindex = this.gpu_core;
+ cfg.option.Connectivity = this.conn;
+ cfg.option.Threshold = this.edge_threshold;
+ cfg.option.Radius = this.radius;
+
+ alg_id = astra_mex_algorithm('create',cfg);
+ astra_mex_algorithm('iterate',alg_id,1);
+ Mask = astra_mex_data2d('get', mask_id);
+
+ astra_mex_algorithm('delete', alg_id);
+ astra_mex_data2d('delete', data_id, mask_id);
+
+ % random
+ RandomField = double(rand(size(S)) < this.random);
+
+ % combine
+ Mask = or(Mask, RandomField);
+
+ end
+
+ %------------------------------------------------------------------
+ function Mask = apply_2D(this, S)
+
+ r = this.radius;
+ w = 2 * r + 1;
+
+ kernel = Kernels.BinaryPixelKernel(r, this.conn);
+
+ % edges
+ Xlarge = zeros(size(S,1)+w-1, size(S,2)+w-1);
+ Xlarge(1+r:end-r, 1+r:end-r) = S;
+
+ Edges = zeros(size(S));
+ for s = -r:r
+ for t = -r:r
+ if kernel(s+r+1, t+r+1) == 0
+ continue
+ end
+ Temp = abs(Xlarge(1+r:end-r, 1+r:end-r) - Xlarge(1+r+s:end-r+s, 1+r+t:end-r+t));
+ Edges(Temp > eps) = Edges(Temp > eps) + 1;
+ end
+ end
+
+ Edges = Edges > this.edge_threshold;
+
+ % random
+ RandomField = double(rand(size(S)) < this.random);
+
+ % combine
+ Mask = or(Edges, RandomField);
+
+ end
+
+ %------------------------------------------------------------------
+ function Mask = apply_3D(this, S)
+
+ r = this.radius;
+ w = 2 * r + 1;
+
+ kernel = Kernels.BinaryPixelKernel(r, this.conn);
+
+ % edges
+ Xlarge = zeros(size(S,1)+w-1, size(S,2)+w-1, size(S,3)+w-1);
+ Xlarge(1+r:end-r, 1+r:end-r, 1+r:end-r) = S;
+
+ Edges = zeros(size(S));
+ for s = -r:r
+ for t = -r:r
+ for u = -r:r
+ if kernel(s+r+1, t+r+1, u+r+1) == 0
+ continue
+ end
+ Temp = abs(Xlarge(1+r:end-r, 1+r:end-r, 1+r:end-r) - Xlarge(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u));
+ Edges(Temp > eps) = 1;
+ end
+ end
+ end
+
+ clear Xlarge;
+
+ % random
+ RandomField = double(rand(size(S)) < this.random);
+
+ % combine
+ Mask = or(Edges, RandomField);
+
+ end
+
+ %------------------------------------------------------------------
+ function Mask = apply_3D_gpu(this, S)
+
+ vol_geom = astra_create_vol_geom(size(S));
+ data_id = astra_mex_data3d('create', '-vol', vol_geom, S);
+
+ cfg = astra_struct('DARTMASK3D_CUDA');
+ cfg.SegmentationDataId = data_id;
+ cfg.MaskDataId = data_id;
+ cfg.option.GPUindex = this.gpu_core;
+ cfg.option.Connectivity = this.conn;
+ cfg.option.Threshold = this.edge_threshold;
+ cfg.option.Radius = this.radius;
+
+ alg_id = astra_mex_algorithm('create',cfg);
+ astra_mex_algorithm('iterate',alg_id,1);
+ Mask = astra_mex_data3d('get', data_id);
+
+ astra_mex_algorithm('delete', alg_id);
+ astra_mex_data3d('delete', data_id);
+
+ % random
+ RandomField = double(rand(size(S)) < this.random);
+
+ % combine
+ Mask = or(Mask, RandomField);
+
+ end
+
+ %------------------------------------------------------------------
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/MaskingGPU.m b/matlab/algorithms/DART/MaskingGPU.m
new file mode 100644
index 0000000..c4ef2b7
--- /dev/null
+++ b/matlab/algorithms/DART/MaskingGPU.m
@@ -0,0 +1,94 @@
+% 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 MaskingGPU < matlab.mixin.Copyable
+
+ % Policy class for masking for DART with GPU accelerated code (deprecated).
+
+ %----------------------------------------------------------------------
+ properties (Access=public)
+
+ radius = 1; % SETTING: Radius of masking kernel.
+ conn = 8; % SETTING: Connectivity window. For 2D: 4 or 8. For 3D: 6 or 26.
+ edge_threshold = 1; % SETTING: Number of pixels in the window that should be different.
+ gpu_core = 0; % SETTING:
+ random = 0.1; % SETTING: Percentage of random points. Between 0 and 1.
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function settings = getsettings(this)
+ % Returns a structure containing all settings of this object.
+ % >> settings = DART.masking.getsettings();
+ settings.radius = this.radius;
+ settings.conn = this.conn;
+ settings.edge_threshold = this.edge_threshold;
+ settings.random = this.random;
+ end
+
+ %------------------------------------------------------------------
+ function Mask = apply(this, ~, V_in)
+ % Applies masking.
+ % >> Mask = DART.segmentation.apply(DART, V_in);
+
+ % 2D, one slice
+ if size(V_in,3) == 1
+ Mask = this.apply_2D(V_in);
+
+ % 3D, slice by slice
+ elseif this.conn == 4 || this.conn == 8
+ Mask = zeros(size(V_in));
+ for slice = 1:size(V_in,3)
+ Mask(:,:,slice) = this.apply_2D(V_in(:,:,slice));
+ end
+
+ % 3D, full
+ else
+ error('Full 3D masking on GPU not implemented.')
+ end
+
+ end
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=protected)
+
+ %------------------------------------------------------------------
+ function Mask = apply_2D(this, S)
+
+ vol_geom = astra_create_vol_geom(size(S));
+ data_id = astra_mex_data2d('create', '-vol', vol_geom, S);
+ mask_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+ cfg = astra_struct('DARTMASK_CUDA');
+ cfg.SegmentationDataId = data_id;
+ cfg.MaskDataId = mask_id;
+ cfg.option.GPUindex = this.gpu_core;
+ %cfg.option.Connectivity = this.conn;
+
+ alg_id = astra_mex_algorithm('create',cfg);
+ astra_mex_algorithm('iterate',alg_id,1);
+ Mask = astra_mex_data2d('get', mask_id);
+
+ astra_mex_algorithm('delete', alg_id);
+ astra_mex_data2d('delete', data_id, mask_id);
+
+ end
+ end
+
+
+
+end
+
diff --git a/matlab/algorithms/DART/OutputDefault.m b/matlab/algorithms/DART/OutputDefault.m
new file mode 100644
index 0000000..a34a430
--- /dev/null
+++ b/matlab/algorithms/DART/OutputDefault.m
@@ -0,0 +1,173 @@
+% 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 OutputDefault < matlab.mixin.Copyable
+
+ % Default policy class for output for DART.
+
+ properties (Access=public)
+
+ directory = ''; % SETTING: Directory to save output.
+ pre = ''; % SETTING: Prefix of output.
+
+ save_images = 'no'; % SETTING: Save the images. 'no', 'yes' (='S') OR {'S', 'I', 'Mask', 'P'}
+ save_results = 'no'; % SETTING: Save the results. 'yes', 'no' OR {'base', 'stats', 'settings', 'S', 'V', 'V0'}
+ save_object = 'no'; % SETTING: Save the DART object. {'no','yes'}
+
+ save_interval = 1; % SETTING: # DART iteration between saves.
+ save_images_interval = []; % SETTING: Overwrite interval for save_images.
+ save_results_interval = []; % SETTING: Overwrite interval for save_results.
+ save_object_interval = []; % SETTING: Overwrite interval for save_object.
+
+ slices = 1; % SETTING: In case of 3D, which slices to save?
+
+ verbose = 'no'; % SETTING: Verbose? {'no','yes'}
+
+ end
+
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function pre_initial_iteration(this, ~)
+ if strcmp(this.verbose,'yes')
+ tic
+ fprintf(1, 'initial iteration...');
+ end
+ end
+
+ %------------------------------------------------------------------
+ function post_initial_iteration(this, ~)
+ if strcmp(this.verbose,'yes')
+ t = toc;
+ fprintf(1, 'done in %f s.\n', t);
+ end
+ end
+
+ %------------------------------------------------------------------
+ function pre_iteration(this, DART)
+ if strcmp(this.verbose,'yes')
+ tic;
+ fprintf(1, '%s dart iteration %d...', this.pre, DART.iterationcount);
+ end
+ end
+
+ %------------------------------------------------------------------
+ function post_iteration(this, DART)
+
+ % print output
+ if strcmp(this.verbose,'yes')
+ t = toc;
+ s = DART.statistics.tostring(DART.stats);
+ fprintf(1, 'done in %0.2fs %s.\n', t, s);
+ end
+
+ % save DART object
+ if do_object(this, DART)
+ save(sprintf('%s%sobject_%i.mat', this.directory, this.pre, DART.iterationcount), '-v7.3', 'DART');
+ end
+
+ % save .mat
+ if do_results(this, DART)
+ base = DART.base;
+ stats = DART.stats;
+ S = DART.S;
+ V = DART.V;
+ V0 = DART.V0;
+ settings = DART.getsettings();
+ if ~iscell(this.save_results)
+ save(sprintf('%s%sresults_%i.mat', this.directory, this.pre, DART.iterationcount), '-v7.3', 'base', 'stats', 'S', 'V', 'V0', 'settings');
+ else
+ string = [];
+ for i = 1:numel(this.save_results)
+ string = [string this.save_results{i} '|'];
+ end
+ save(sprintf('%s%sresults_%i.mat', this.directory, this.pre, DART.iterationcount), '-v7.3', '-regexp', string(1:end-1));
+ end
+ end
+
+ % save images
+ if do_images(this, DART)
+
+ if ~iscell(this.save_images) && strcmp(this.save_images, 'yes')
+ output_image(this, DART, 'S')
+ elseif iscell(this.save_images)
+ for i = 1:numel(this.save_images)
+ output_image(this, DART, this.save_images{i});
+ end
+ end
+
+ end
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=private)
+
+ function output_image(this, DART, data)
+ % 2D
+ if numel(size(DART.S)) == 2
+ eval(['imwritesc(DART.' data ', sprintf(''%s%s' data '_%i.png'', this.directory, this.pre, DART.iterationcount))']);
+ % 3D
+ elseif numel(size(DART.S)) == 3
+ for slice = this.slices
+ eval(['imwritesc(DART.' data '(:,:,slice), sprintf(''%s%s' data '_%i_slice%i.png'', this.directory, this.pre, DART.iterationcount, slice))']);
+ end
+ end
+ end
+
+ %------------------------------------------------------------------
+ function out = do_object(this, DART)
+ if strcmp(this.save_object,'no')
+ out = 0;
+ return
+ end
+ if numel(this.save_object_interval) == 0 && mod(DART.iterationcount, this.save_interval) == 0
+ out = 1;
+ elseif mod(DART.iterationcount, this.save_object_interval) == 0
+ out = 1;
+ else
+ out = 0;
+ end
+ end
+ %------------------------------------------------------------------
+ function out = do_results(this, DART)
+ if strcmp(this.save_results,'no')
+ out = 0;
+ return
+ end
+ if numel(this.save_results_interval) == 0 && mod(DART.iterationcount, this.save_interval) == 0
+ out = 1;
+ elseif mod(DART.iterationcount, this.save_results_interval) == 0
+ out = 1;
+ else
+ out = 0;
+ end
+ end
+
+ %------------------------------------------------------------------
+ function out = do_images(this, DART)
+ if numel(this.save_images_interval) == 0 && mod(DART.iterationcount, this.save_interval) == 0
+ out = 1;
+ elseif mod(DART.iterationcount, this.save_images_interval) == 0
+ out = 1;
+ else
+ out = 0;
+ end
+ end
+ %------------------------------------------------------------------
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/SegmentationDefault.m b/matlab/algorithms/DART/SegmentationDefault.m
new file mode 100644
index 0000000..c1d7d99
--- /dev/null
+++ b/matlab/algorithms/DART/SegmentationDefault.m
@@ -0,0 +1,55 @@
+% 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 SegmentationDefault < matlab.mixin.Copyable
+
+ % Default policy class for segmentation for DART.
+
+ %----------------------------------------------------------------------
+ properties (Access=public)
+ rho = []; % SETTING: Grey levels.
+ tau = []; % SETTING: Threshold values.
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function settings = getsettings(this)
+ % Returns a structure containing all settings of this object.
+ % >> settings = DART.segmentation.getsettings();
+ settings.rho = this.rho;
+ settings.tau = this.tau;
+ end
+
+ %------------------------------------------------------------------
+ function this = estimate_grey_levels(this, ~, ~)
+ % Estimates grey levels
+ % >> DART.segmentation.estimate_grey_levels();
+ end
+
+ %------------------------------------------------------------------
+ function V_out = apply(this, ~, V_in)
+ % Applies segmentation.
+ % >> V_out = DART.segmentation.apply(DART, V_in);
+
+ V_out = ones(size(V_in)) * this.rho(1);
+ for n = 2:length(this.rho)
+ V_out(this.tau(n-1) < V_in) = this.rho(n);
+ end
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/SmoothingDefault.m b/matlab/algorithms/DART/SmoothingDefault.m
new file mode 100644
index 0000000..58a8baa
--- /dev/null
+++ b/matlab/algorithms/DART/SmoothingDefault.m
@@ -0,0 +1,179 @@
+% 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 SmoothingDefault < matlab.mixin.Copyable
+
+ % Default policy class for smoothing for DART.
+
+ %----------------------------------------------------------------------
+ properties (Access=public)
+ radius = 1; % SETTING: Radius of smoothing kernel.
+ b = 0.1; % SETTING: Intensity of smoothing. Between 0 and 1.
+ full3d = 'yes'; % SETTING: smooth in 3D? {'yes','no'}
+ gpu = 'yes'; % SETTING: Use gpu? {'yes', 'no'}
+ gpu_core = 0; % SETTING: Which gpu core to use, only when gpu='yes'.
+ end
+
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function settings = getsettings(this)
+ % Returns a structure containing all settings of this object.
+ % >> settings = DART.smoothing.getsettings();
+ settings.radius = this.radius;
+ settings.b = this.b;
+ settings.full3d = this.full3d;
+ end
+
+ %------------------------------------------------------------------
+ function V_out = apply(this, ~, V_in)
+ % Applies smoothing.
+ % >> V_out = DART.smoothing.apply(DART, V_in);
+
+ % 2D, one slice
+ if size(V_in,3) == 1
+ if strcmp(this.gpu,'yes')
+ V_out = this.apply_2D_gpu(V_in);
+ else
+ V_out = this.apply_2D(V_in);
+ end
+
+ % 3D, slice by slice
+ elseif ~strcmp(this.full3d,'yes')
+ V_out = zeros(size(V_in));
+ for slice = 1:size(V_in,3)
+ if strcmp(this.gpu,'yes')
+ V_out(:,:,slice) = this.apply_2D_gpu(V_in(:,:,slice));
+ else
+ V_out(:,:,slice) = this.apply_2D(V_in(:,:,slice));
+ end
+ end
+
+ % 3D, full
+ else
+ if strcmp(this.gpu,'yes')
+ V_out = this.apply_3D_gpu(V_in);
+ else
+ V_out = this.apply_3D(V_in);
+ end
+ end
+
+ end
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=protected)
+
+ %------------------------------------------------------------------
+ function V_out = apply_2D(this, V_in)
+
+ r = this.radius;
+ w = 2 * r + 1;
+
+ % Set Kernel
+ K = ones(w) * this.b / (w.^2-1); % edges
+ K(r+1,r+1) = 1 - this.b; % center
+
+ % output window
+ V_out = zeros(size(V_in,1) + w-1, size(V_in,2) + w - 1);
+
+ % blur convolution
+ for s = -r:r
+ for t = -r:r
+ V_out(1+r+s:end-r+s, 1+r+t:end-r+t) = V_out(1+r+s:end-r+s, 1+r+t:end-r+t) + K(r+1+s, r+1+t) * V_in;
+ end
+ end
+
+ % shrink output window
+ V_out = V_out(1+r:end-r, 1+r:end-r);
+
+ end
+
+ %------------------------------------------------------------------
+ function V_out = apply_2D_gpu(this, V_in)
+
+ vol_geom = astra_create_vol_geom(size(V_in));
+ in_id = astra_mex_data2d('create', '-vol', vol_geom, V_in);
+ out_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+ cfg = astra_struct('DARTSMOOTHING_CUDA');
+ cfg.InDataId = in_id;
+ cfg.OutDataId = out_id;
+ cfg.option.Intensity = this.b;
+ cfg.option.Radius = this.radius;
+ cfg.option.GPUindex = this.gpu_core;
+
+ alg_id = astra_mex_algorithm('create',cfg);
+ astra_mex_algorithm('iterate',alg_id,1);
+ V_out = astra_mex_data2d('get', out_id);
+
+ astra_mex_algorithm('delete', alg_id);
+ astra_mex_data2d('delete', in_id, out_id);
+
+ end
+
+ %------------------------------------------------------------------
+ function I_out = apply_3D(this, I_in)
+
+ r = this.radius;
+ w = 2 * r + 1;
+
+ % Set Kernel
+ K = ones(w,w,w) * this.b / (w.^3-1); % edges
+ K(r+1,r+1,r+1) = 1 - this.b; % center
+
+ % output window
+ I_out = zeros(size(I_in,1)+w-1, size(I_in,2)+w-1, size(I_in,3)+w-1);
+
+ % blur convolution
+ for s = -r:r
+ for t = -r:r
+ for u = -r:r
+ I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) = I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) + K(r+1+s, r+1+t, r+1+u) * I_in;
+ end
+ end
+ end
+
+ % shrink output window
+ I_out = I_out(1+r:end-r, 1+r:end-r, 1+r:end-r);
+
+ end
+
+ %------------------------------------------------------------------
+ function V_out = apply_3D_gpu(this, V_in)
+
+ vol_geom = astra_create_vol_geom(size(V_in));
+ data_id = astra_mex_data3d('create', '-vol', vol_geom, V_in);
+
+ cfg = astra_struct('DARTSMOOTHING3D_CUDA');
+ cfg.InDataId = data_id;
+ cfg.OutDataId = data_id;
+ cfg.option.Intensity = this.b;
+ cfg.option.Radius = this.radius;
+ cfg.option.GPUindex = this.gpu_core;
+
+ alg_id = astra_mex_algorithm('create',cfg);
+ astra_mex_algorithm('iterate', alg_id, 1);
+ V_out = astra_mex_data3d('get', data_id);
+
+ astra_mex_algorithm('delete', alg_id);
+ astra_mex_data3d('delete', data_id);
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/SmoothingGPU.m b/matlab/algorithms/DART/SmoothingGPU.m
new file mode 100644
index 0000000..857da37
--- /dev/null
+++ b/matlab/algorithms/DART/SmoothingGPU.m
@@ -0,0 +1,119 @@
+% 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 SmoothingGPU < matlab.mixin.Copyable
+
+ % Default policy class for smoothing for DART.
+
+ %----------------------------------------------------------------------
+ properties (Access=public)
+ radius = 1; % SETTING: Radius of smoothing kernel.
+ b = 0.1; % SETTING: Intensity of smoothing. Between 0 and 1.
+ full3d = 'yes'; % SETTING: smooth in 3D? {'yes','no'}
+ gpu_core = 0; % SETTING:
+ end
+
+
+ %----------------------------------------------------------------------
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function settings = getsettings(this)
+ % Returns a structure containing all settings of this object.
+ % >> settings = DART.smoothing.getsettings();
+ settings.radius = this.radius;
+ settings.b = this.b;
+ settings.full3d = this.full3d;
+ end
+
+ %------------------------------------------------------------------
+ function V_out = apply(this, ~, V_in)
+ % Applies smoothing.
+ % >> V_out = DART.smoothing.apply(DART, V_in);
+
+ % 2D, one slice
+ if size(V_in,3) == 1
+ V_out = this.apply_2D(V_in);
+
+ % 3D, slice by slice
+ elseif ~strcmp(this.full3d,'yes')
+ V_out = zeros(size(V_in));
+ for slice = 1:size(V_in,3)
+ V_out(:,:,slice) = this.apply_2D(V_in(:,:,slice));
+ end
+
+ % 3D, full
+ else
+ V_out = this.apply_3D(V_in);
+ end
+
+ end
+
+ end
+
+ %----------------------------------------------------------------------
+ methods (Access=protected)
+
+ %------------------------------------------------------------------
+ function V_out = apply_2D(this, V_in)
+
+ vol_geom = astra_create_vol_geom(size(V_in));
+ in_id = astra_mex_data2d('create', '-vol', vol_geom, V_in);
+ out_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+ cfg = astra_struct('DARTSMOOTHING_CUDA');
+ cfg.InDataId = in_id;
+ cfg.OutDataId = out_id;
+ cfg.Intensity = this.b;
+ cfg.option.GPUindex = this.gpu_core;
+
+ alg_id = astra_mex_algorithm('create',cfg);
+ astra_mex_algorithm('iterate',alg_id,1);
+ V_out = astra_mex_data2d('get', out_id);
+
+ astra_mex_algorithm('delete', alg_id);
+ astra_mex_data2d('delete', in_id, out_id);
+
+
+ end
+
+ %------------------------------------------------------------------
+ function I_out = apply_3D(this, I_in)
+
+ r = this.radius;
+ w = 2 * r + 1;
+
+ % Set Kernel
+ K = ones(w,w,w) * this.b / (w.^3-1); % edges
+ K(r+1,r+1,r+1) = 1 - this.b; % center
+
+ % output window
+ I_out = zeros(size(I_in,1)+w-1, size(I_in,2)+w-1, size(I_in,3)+w-1);
+
+ % blur convolution
+ for s = -r:r
+ for t = -r:r
+ for u = -r:r
+ I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) = I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) + K(r+1+s, r+1+t, r+1+u) * I_in;
+ end
+ end
+ end
+
+ % shrink output window
+ I_out = I_out(1+r:end-r, 1+r:end-r, 1+r:end-r);
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/StatisticsDefault.m b/matlab/algorithms/DART/StatisticsDefault.m
new file mode 100644
index 0000000..7822c5f
--- /dev/null
+++ b/matlab/algorithms/DART/StatisticsDefault.m
@@ -0,0 +1,72 @@
+% 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 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'}
+ end
+
+
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function stats = apply(this, DART)
+ % Applies statistics.
+ % >> stats = DART.statistics.apply(DART);
+
+ stats = DART.stats;
+
+ % timing
+ if strcmp(this.timing, 'yes')
+ stats.timing(DART.iterationcount) = toc(DART.start_tic);
+ end
+
+ % pixel error
+ if strcmp(this.pixel_error, 'yes') && isfield(DART.base,'phantom')
+ [stats.rnmp, stats.nmp] = compute_rnmp(DART.base.phantom, DART.S);
+ stats.rnmp_hist(DART.iterationcount) = stats.rnmp;
+ stats.nmp_hist(DART.iterationcount) = stats.nmp;
+ end
+
+ % projection difference
+ if strcmp(this.proj_diff, 'yes')
+ new_sino = DART.tomography.createForwardProjection(DART, DART.S);
+ stats.proj_diff = sum((new_sino(:) - DART.base.sinogram(:)) .^2 ) ./ (sum(DART.base.sinogram(:)) );
+ stats.proj_diff_hist(DART.iterationcount) = stats.proj_diff;
+ end
+
+ end
+
+ %------------------------------------------------------------------
+ function s = tostring(~, stats)
+ % To string.
+ % >> stats = DART.statistics.apply(stats);
+
+ s = '';
+ if isfield(stats, 'nmp')
+ s = sprintf('%s [%d]', s, stats.nmp);
+ end
+ if isfield(stats, 'proj_diff')
+ s = sprintf('%s {%0.2d}', s, stats.proj_diff);
+ end
+
+ end
+ %------------------------------------------------------------------
+
+ end
+
+end
+
diff --git a/matlab/algorithms/DART/TomographyDefault.m b/matlab/algorithms/DART/TomographyDefault.m
new file mode 100644
index 0000000..4db3905
--- /dev/null
+++ b/matlab/algorithms/DART/TomographyDefault.m
@@ -0,0 +1,73 @@
+% 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
new file mode 100644
index 0000000..2be1b17
--- /dev/null
+++ b/matlab/algorithms/DART/TomographyDefault3D.m
@@ -0,0 +1,73 @@
+% 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/cylinders.png b/matlab/algorithms/DART/examples/cylinders.png
new file mode 100644
index 0000000..8d1c0b2
--- /dev/null
+++ b/matlab/algorithms/DART/examples/cylinders.png
Binary files differ
diff --git a/matlab/algorithms/DART/examples/example1.m b/matlab/algorithms/DART/examples/example1.m
new file mode 100644
index 0000000..daa3ce8
--- /dev/null
+++ b/matlab/algorithms/DART/examples/example1.m
@@ -0,0 +1,79 @@
+clear all;
+
+addpath('..');
+
+%
+% Example 1: parallel beam, three slices.
+%
+
+% Configuration
+proj_count = 20;
+slice_count = 3;
+dart_iterations = 20;
+filename = 'cylinders.png';
+outdir = './';
+prefix = 'example1';
+rho = [0, 1];
+tau = 0.5;
+gpu_core = 0;
+
+% Load phantom.
+I = double(imread(filename)) / 255;
+
+% 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);
+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);
+
+%
+% DART
+%
+
+base.sinogram = sinogram;
+base.proj_geom = proj_geom;
+
+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.segmentation.rho = rho;
+D.segmentation.tau = tau;
+
+D.smoothing.b = 0.1;
+D.smoothing.full3d = 'yes';
+D.smoothing.gpu_core = gpu_core;
+
+D.masking.random = 0.1;
+D.masking.conn = 6;
+D.masking.gpu_core = gpu_core;
+
+D.output.directory = outdir;
+D.output.pre = [prefix '_'];
+D.output.save_images = 'no';
+D.output.save_results = {'stats', 'settings', 'S', 'V'};
+D.output.save_interval = dart_iterations;
+D.output.verbose = 'yes';
+
+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']);
diff --git a/matlab/algorithms/DART/examples/example2.m b/matlab/algorithms/DART/examples/example2.m
new file mode 100644
index 0000000..8ee8cba
--- /dev/null
+++ b/matlab/algorithms/DART/examples/example2.m
@@ -0,0 +1,80 @@
+clear all;
+
+addpath('..');
+
+%
+% Example 2: cone beam, full cube.
+%
+
+% Configuration
+det_count = 128;
+proj_count = 45;
+slice_count = det_count;
+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.
+[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(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.segmentation.rho = rho;
+D.segmentation.tau = tau;
+
+D.smoothing.b = 0.1;
+D.smoothing.full3d = 'yes';
+D.smoothing.gpu_core = gpu_core;
+
+D.masking.random = 0.1;
+D.masking.conn = 6;
+D.masking.gpu_core = gpu_core;
+
+D.output.directory = outdir;
+D.output.pre = [prefix '_'];
+D.output.save_images = 'no';
+D.output.save_results = {'stats', 'settings', 'S', 'V'};
+D.output.save_interval = dart_iterations;
+D.output.verbose = 'yes';
+
+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']);
diff --git a/matlab/algorithms/DART/examples/example3.m b/matlab/algorithms/DART/examples/example3.m
new file mode 100644
index 0000000..f6e360e
--- /dev/null
+++ b/matlab/algorithms/DART/examples/example3.m
@@ -0,0 +1,79 @@
+clear all;
+
+addpath('..');
+
+%
+% Example 3: parallel beam, 2D
+%
+
+% Configuration
+proj_count = 30;
+dart_iterations = 20;
+filename = 'cylinders.png';
+outdir = './';
+prefix = 'example3';
+rho = [0, 1];
+tau = 0.5;
+gpu_core = 0;
+
+% Load phantom.
+I = double(imread(filename)) / 255;
+
+% 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);
+
+% Create sinogram.
+[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;
+
+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.segmentation.rho = rho;
+D.segmentation.tau = tau;
+
+D.smoothing.b = 0.1;
+D.smoothing.full3d = 'yes';
+D.smoothing.gpu_core = gpu_core;
+
+D.masking.random = 0.1;
+D.masking.conn = 6;
+D.masking.gpu_core = gpu_core;
+
+D.output.directory = outdir;
+D.output.pre = [prefix '_'];
+D.output.save_images = 'no';
+D.output.save_results = {'stats', 'settings', 'S', 'V'};
+D.output.save_interval = dart_iterations;
+D.output.verbose = 'yes';
+
+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']);
diff --git a/matlab/algorithms/DART/examples/phantom3d.mat b/matlab/algorithms/DART/examples/phantom3d.mat
new file mode 100644
index 0000000..6d70c16
--- /dev/null
+++ b/matlab/algorithms/DART/examples/phantom3d.mat
Binary files differ