diff options
Diffstat (limited to 'matlab/algorithms/DART')
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 Binary files differnew file mode 100644 index 0000000..8d1c0b2 --- /dev/null +++ b/matlab/algorithms/DART/examples/cylinders.png 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 Binary files differnew file mode 100644 index 0000000..6d70c16 --- /dev/null +++ b/matlab/algorithms/DART/examples/phantom3d.mat |
