summaryrefslogtreecommitdiffstats
path: root/matlab
diff options
context:
space:
mode:
authorWim van Aarle <wimvanaarle@gmail.com>2015-02-23 16:27:33 +0100
committerWim van Aarle <wimvanaarle@gmail.com>2015-02-23 16:27:33 +0100
commitb16c26d1fd7a50fda648dd63076e011d122543b8 (patch)
treebb5580d86fddfc2c62919b27cbb668ed007115b2 /matlab
parent1727dee5b0a1c5a932e42af0770a0702c1b1b740 (diff)
downloadastra-b16c26d1fd7a50fda648dd63076e011d122543b8.tar.gz
astra-b16c26d1fd7a50fda648dd63076e011d122543b8.tar.bz2
astra-b16c26d1fd7a50fda648dd63076e011d122543b8.tar.xz
astra-b16c26d1fd7a50fda648dd63076e011d122543b8.zip
fixed matlab DART implementation
some files were missing, cleaned up the rest
Diffstat (limited to 'matlab')
-rw-r--r--matlab/algorithms/DART/DARTalgorithm.m11
-rw-r--r--matlab/algorithms/DART/IterativeTomography.m435
-rw-r--r--matlab/algorithms/DART/IterativeTomography3D.m405
-rw-r--r--matlab/algorithms/DART/examples/example1.m40
-rw-r--r--matlab/algorithms/DART/examples/example2.m55
-rw-r--r--matlab/algorithms/DART/examples/example3.m21
6 files changed, 885 insertions, 82 deletions
diff --git a/matlab/algorithms/DART/DARTalgorithm.m b/matlab/algorithms/DART/DARTalgorithm.m
index 85a3ca0..fc707dd 100644
--- a/matlab/algorithms/DART/DARTalgorithm.m
+++ b/matlab/algorithms/DART/DARTalgorithm.m
@@ -9,13 +9,8 @@
%--------------------------------------------------------------------------
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)
@@ -78,7 +73,6 @@ classdef DARTalgorithm < matlab.mixin.Copyable
error('invalid arguments')
end
end
-
%------------------------------------------------------------------
function D = deepcopy(this)
@@ -100,7 +94,6 @@ classdef DARTalgorithm < matlab.mixin.Copyable
% Initialize tomography part
if ~this.tomography.initialized
- this.tomography.sinogram = this.base.sinogram;
this.tomography.proj_geom = this.base.proj_geom;
this.tomography.initialize();
end
@@ -110,7 +103,7 @@ classdef DARTalgorithm < matlab.mixin.Copyable
this.V0 = this.base.V0;
else
this.output.pre_initial_iteration(this);
- this.V0 = this.tomography.reconstruct2(this.base.sinogram, [], this.t0);
+ this.V0 = this.tomography.reconstruct(this.base.sinogram, this.t0);
this.output.post_initial_iteration(this);
end
this.V = this.V0;
@@ -163,7 +156,7 @@ classdef DARTalgorithm < matlab.mixin.Copyable
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);
+ this.V = this.tomography.reconstruct_mask(this.R, this.V, this.Mask, this.t);
% blur
this.V = this.smoothing.apply(this, this.V);
diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
new file mode 100644
index 0000000..3875e6b
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography.m
@@ -0,0 +1,435 @@
+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)
+ proj_geom = []; % DATA: Projection geometry.
+ vol_geom = []; % DATA: Volume geometry.
+ end
+ %----------------------------------------------------------------------
+ properties (SetAccess=private, GetAccess=public)
+ initialized = 0; % Is this object initialized?
+ end
+ %----------------------------------------------------------------------
+ properties (SetAccess=protected, GetAccess=protected)
+ cfg_base = struct(); % PROTECTED: base configuration structure for the reconstruction algorithm.
+ 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')
+ end
+ %----------------------------------------------------------------------
+
+ methods (Access=public)
+
+ %------------------------------------------------------------------
+ function this = IterativeTomography(varargin)
+ % Constructor
+ % >> tomography = IterativeTomography(proj_geom);
+ % >> tomography = IterativeTomography(proj_geom, vol_geom);
+
+ % Input: IterativeTomography(proj_geom)
+ if nargin == 1
+ this.proj_geom = varargin{1};
+
+ % Input: IterativeTomography(proj_geom, vol_geom)
+ elseif nargin == 2
+ this.proj_geom = varargin{1};
+ this.vol_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
+ if this.superresolution > 1
+ this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);
+ else
+ this.proj_geom_sr = this.proj_geom;
+ end
+
+ % 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;
+ end
+ this.cfg_base.option.DetectorSuperSampling = this.superresolution;
+ if strcmp(this.gpu,'yes')
+ this.cfg_base.option.GPUindex = this.gpu_core;
+ end
+ if this.use_minc
+ this.cfg_base.option.MinConstraint = 0;
+ end
+
+ this.initialized = 1;
+ ok = this.initialized;
+ end
+
+ %------------------------------------------------------------------
+ function projections = project(this, volume)
+ % Compute forward projection.
+ % >> projections = tomography.project(volume);
+
+ if ~this.initialized
+ this.initialize();
+ end
+
+ % project
+ projections = this.project_c(volume);
+ end
+
+ %------------------------------------------------------------------
+ function reconstruction = reconstruct(this, varargin)
+ % Compute reconstruction.
+ % Uses tomography.sinogram
+ % Initial solution (if available) should be stored in tomography.V
+ % >> reconstruction = tomography.reconstruct(projections, iterations);
+ % >> reconstruction = tomography.reconstruct(projections, volume0, iterations);
+
+ if ~this.initialized
+ this.initialize();
+ end
+
+ if numel(varargin) == 2
+ reconstruction = this.reconstruct_c(varargin{1}, [], [], varargin{2});
+ elseif numel(varargin) == 3
+ reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, [], varargin{3});
+ else
+ error('invalid parameter list')
+ end
+
+ if strcmp(this.inner_circle,'yes')
+ reconstruction = this.selectROI(reconstruction);
+ end
+ end
+
+ %------------------------------------------------------------------
+ function reconstruction = reconstruct_mask(this, varargin)
+ % Compute reconstruction with mask.
+ % Uses tomography.sinogram
+ % Initial solution (if available) should be stored in tomography.V
+ % >> reconstruction = tomography.reconstructMask(projections, mask, iterations);
+ % >> reconstruction = tomography.reconstructMask(projections, volume0, mask, iterations);
+
+ if ~this.initialized
+ this.initialize();
+ end
+
+ if numel(varargin) == 3
+ reconstruction = this.reconstruct_c(varargin{1}, [], varargin{2}, varargin{3});
+ elseif numel(varargin) == 4
+ reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, varargin{3}, varargin{4});
+ else
+ error('invalid parameter list')
+ end
+
+ if strcmp(this.inner_circle,'yes')
+ reconstruction = this.selectROI(reconstruction);
+ end
+
+ 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);
+ 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)
+ this.cfg_base.option.ProjectionOrder = 'random';
+ 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..29b963f
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography3D.m
@@ -0,0 +1,405 @@
+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 = IterativeTomography3D(varargin)
+ % Constructor
+ % >> tomography = IterativeTomography3D(proj_geom);
+ % >> tomography = IterativeTomography3D(proj_geom, vol_geom);
+
+ % Input: IterativeTomography(proj_geom)
+ if nargin == 1
+ this.proj_geom = varargin{1};
+
+ % Input: IterativeTomography(proj_geom, vol_geom)
+ elseif nargin == 2
+ this.proj_geom = varargin{1};
+ this.vol_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;
+ 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 projections = project(this, volume)
+ % Compute forward projection.
+ % >> projections = tomography.project(volume);
+
+ if ~this.initialized
+ this.initialize();
+ end
+
+ % project
+ projections = this.project_c(volume);
+ end
+
+ %------------------------------------------------------------------
+ function reconstruction = reconstruct(this, varargin)
+ % Compute reconstruction.
+ % Uses tomography.sinogram
+ % Initial solution (if available) should be stored in tomography.V
+ % >> reconstruction = tomography.reconstruct(projections, iterations);
+ % >> reconstruction = tomography.reconstruct(projections, volume0, iterations);
+
+ if ~this.initialized
+ this.initialize();
+ end
+
+ if numel(varargin) == 2
+ reconstruction = this.reconstruct_c(varargin{1}, [], [], varargin{2});
+ elseif numel(varargin) == 3
+ reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, [], varargin{3});
+ else
+ error('invalid parameter list')
+ end
+
+ if strcmp(this.inner_circle,'yes')
+ reconstruction = this.selectROI(reconstruction);
+ end
+ end
+
+ %------------------------------------------------------------------
+ function reconstruction = reconstruct_mask(this, varargin)
+ % Compute reconstruction with mask.
+ % Uses tomography.sinogram
+ % Initial solution (if available) should be stored in tomography.V
+ % >> reconstruction = tomography.reconstructMask(projections, mask, iterations);
+ % >> reconstruction = tomography.reconstructMask(projections, volume0, mask, iterations);
+
+ if ~this.initialized
+ this.initialize();
+ end
+
+ if numel(varargin) == 3
+ reconstruction = this.reconstruct_c(varargin{1}, [], varargin{2}, varargin{3});
+ elseif numel(varargin) == 4
+ reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, varargin{3}, varargin{4});
+ else
+ error('invalid parameter list')
+ end
+
+ if strcmp(this.inner_circle,'yes')
+ reconstruction = this.selectROI(reconstruction);
+ end
+
+ 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);
+ 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/examples/example1.m b/matlab/algorithms/DART/examples/example1.m
index 9a836f8..cb02e0f 100644
--- a/matlab/algorithms/DART/examples/example1.m
+++ b/matlab/algorithms/DART/examples/example1.m
@@ -8,15 +8,11 @@
% Website: http://sf.net/projects/astra-toolbox
%--------------------------------------------------------------------------
-clear all;
-
addpath('../');
-%%
-% Example 1: 2D parallel beam, cuda
-%
+%% Example 1: 2D parallel beam, cuda
-% Configuration
+% configuration
proj_count = 20;
dart_iterations = 20;
filename = 'cylinders.png';
@@ -26,26 +22,20 @@ rho = [0, 255];
tau = 128;
gpu_core = 0;
-% Load phantom
+% load phantom
I = imreadgs(filename);
-% Create projection and volume geometries
+% create projection and volume geometries
det_count = size(I, 1);
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);
+vol_geom = astra_create_vol_geom(det_count, det_count);
-% Create sinogram.
+% 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(sinogram, proj_geom);
D.t0 = 100;
D.t = 10;
@@ -63,21 +53,19 @@ D.smoothing.gpu_core = gpu_core;
D.masking.random = 0.1;
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.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.statistics.proj_diff = 'no';
D.initialize();
D.iterate(dart_iterations);
-%%
-% Convert middle slice of final iteration to png.
-%
+% save the reconstruction and the segmentation to file
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 7fe6988..89660a5 100644
--- a/matlab/algorithms/DART/examples/example2.m
+++ b/matlab/algorithms/DART/examples/example2.m
@@ -8,15 +8,11 @@
% Website: http://sf.net/projects/astra-toolbox
%--------------------------------------------------------------------------
-clear all;
-
addpath('../');
-%%
-% Example Z: 3D parallel beam, cuda
-%
+%% Example 2: 3D parallel beam, cuda
-% Configuration
+% configuration
proj_count = 20;
dart_iterations = 20;
outdir = './';
@@ -25,50 +21,47 @@ rho = [0, 0.5, 1];
tau = [0.25, 0.75];
gpu_core = 0;
-% Load phantom
+% load phantom
load('phantom3d');
-% Create projection and volume geometries
+% 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
+% create sinogram
[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
astra_mex_data3d('delete', sinogram_id);
-%%
% DART
-%
-
D = DARTalgorithm(sinogram, proj_geom);
D.t0 = 100;
D.t = 10;
D.tomography = IterativeTomography3D();
-D.tomography.method = 'SIRT3D_CUDA';
-D.tomography.gpu_core = gpu_core;
-D.tomography.use_minc = 'yes';
+D.tomography.method = 'SIRT3D_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.full3d = 'yes';
+D.smoothing.gpu_core = gpu_core;
-D.masking.random = 0.1;
-D.masking.conn = 4;
-D.masking.gpu_core = gpu_core;
+D.masking.random = 0.1;
+D.masking.conn = 4;
+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.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';
@@ -76,8 +69,6 @@ D.initialize();
D.iterate(dart_iterations);
-%%
-% Convert output of final iteration to png.
-%
+% save the central slice of the reconstruction and the segmentation to file
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 895630b..cc80b0f 100644
--- a/matlab/algorithms/DART/examples/example3.m
+++ b/matlab/algorithms/DART/examples/example3.m
@@ -8,15 +8,11 @@
% Website: http://sf.net/projects/astra-toolbox
%--------------------------------------------------------------------------
-clear all;
-
addpath('../');
-%%
-% Example 3: 3D cone beam, cuda
-%
+%% Example 3: 3D cone beam, cuda
-% Configuration
+% configuration
proj_count = 20;
dart_iterations = 20;
outdir = './';
@@ -25,24 +21,21 @@ rho = [0, 0.5, 1];
tau = [0.25, 0.75];
gpu_core = 0;
-% Load phantom
+% load phantom
load('phantom3d');
-% Create projection and volume geometries
+% 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('cone', 1, 1, slice_count, det_count, angles, 500, 0);
vol_geom = astra_create_vol_geom(size(I));
-% Create sinogram
+% create sinogram
[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
astra_mex_data3d('delete', sinogram_id);
-%%
% DART
-%
-
D = DARTalgorithm(sinogram, proj_geom);
D.t0 = 100;
D.t = 10;
@@ -76,8 +69,6 @@ D.initialize();
D.iterate(dart_iterations);
-%%
-% Convert output of final iteration to png.
-%
+% save the central slice of the reconstruction and the segmentation to file
imwritesc(D.S(:,:,64), [outdir '/' prefix '_S_slice_64.png']);
imwritesc(D.V(:,:,64), [outdir '/' prefix '_V_slice_64.png']);