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