diff options
| -rw-r--r-- | matlab/algorithms/DART/DARTalgorithm.m | 11 | ||||
| -rw-r--r-- | matlab/algorithms/DART/IterativeTomography.m | 435 | ||||
| -rw-r--r-- | matlab/algorithms/DART/IterativeTomography3D.m | 405 | ||||
| -rw-r--r-- | matlab/algorithms/DART/examples/example1.m | 40 | ||||
| -rw-r--r-- | matlab/algorithms/DART/examples/example2.m | 55 | ||||
| -rw-r--r-- | matlab/algorithms/DART/examples/example3.m | 21 |
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']); |