diff --git a/MATLAB/Algorithms/ASD_POCS.m b/MATLAB/Algorithms/ASD_POCS.m index 8c314f99..51b56a6c 100644 --- a/MATLAB/Algorithms/ASD_POCS.m +++ b/MATLAB/Algorithms/ASD_POCS.m @@ -1,425 +1,75 @@ -function [ f,qualMeasOut ] = ASD_POCS(proj,geo,angles,maxiter,varargin) -%ASD_POCS Solves the ASD_POCS total variation constrained image in 3D -% tomography. -% -% ASD_POCS(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem -% using the projection data PROJ taken over ALPHA angles, corresponding -% to the geometry described in GEO, using NITER iterations. -% -% ASD_POCS(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The -% possible options in OPT are: -% -% -% 'lambda': Sets the value of the hyperparameter for the SART iterations. -% Default is 1 -% -% 'lambdared': Reduction of lambda. Every iteration -% lambda = lambdared * lambda. Default is 0.99 -% -% 'init': Describes different initialization techniques. -% • 'none' : Initializes the image to zeros (default) -% -% • 'FDK' : Initializes image to FDK reconstruction -% 'TViter': Defines the amount of TV iterations performed per SART -% iteration. Default is 20 -% -% 'alpha': Defines the TV hyperparameter. Default is 0.002 -% -% 'alpha_red': Defines the reduction rate of the TV hyperparameter -% -% 'Ratio': The maximum allowed image/TV update ratio. If the TV -% update changes the image more than this, the parameter -% will be reduced. Default is 0.95 -% 'maxL2err' Maximum L2 error to accept an image as valid. This -% parameter is crucial for the algorithm, determines at -% what point an image should not be updated further. -% Default is 20% of the FDK L2 norm. -% 'Verbose' 1 or 0. Default is 1. Gives information about the -% progress of the algorithm. -% -% 'OrderStrategy' Chooses the subset ordering strategy. Options are -% 'ordered' : uses them in the input order, but divided -% 'random' : orders them randomly -% 'angularDistance': chooses the next subset with the -% biggest angular distance with the ones used. -% 'redundancy_weighting': true or false. Default is true. Applies data -% redundancy weighting to projections in the update step -% (relevant for offset detector geometry) -% 'groundTruth' an image as ground truth, to be used if quality measures -% are requested, to plot their change w.r.t. this known -% data. -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and -% CERN-European Organization for Nuclear Research -% All rights reserved. -% -% License: Open Source under BSD. -% See the full license at -% https://github.com/CERN/TIGRE/blob/master/LICENSE -% -% Contact: tigre.toolbox@gmail.com -% Codes: https://github.com/CERN/TIGRE/ -% Coded by: Ander Biguri and Manasavee Lohvithee -%-------------------------------------------------------------------------- - - - -%% parse inputs -blocksize=1; -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); - -measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); -if ~any(isnan(gt(:))) - QualMeasOpts{end+1}='error_norm'; - res_prev=gt; - clear gt -end -if nargout<2 && measurequality - warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them") - measurequality=false; -end -qualMeasOut=zeros(length(QualMeasOpts),maxiter); - - -[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); -angles_reorder=cell2mat(alphablocks); -index_angles=cell2mat(orig_index); - -% does detector rotation exists? -if ~isfield(geo,'rotDetector') - geo.rotDetector=[0;0;0]; -end -%% Create weighting matrices for the SART step -% the reason we do this, instead of calling the SART function is not to -% recompute the weights every ASD-POCS iteration, thus effectively doubling -% the computational time -% Projection weight, W -W=computeW(geo,angles,gpuids); - - -% Back-Projection weight, V -V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids); - -if redundancy_weights - % Data redundancy weighting, W_r implemented using Wang weighting - % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc - - num_frames = size(proj,3); - W_r = redundancy_weighting(geo); - W_r = repmat(W_r,[1,1,num_frames]); - % disp('Size of redundancy weighting matrix'); - % disp(size(W_r)); - W = W.*W_r; % include redundancy weighting in W -end - - -%% -stop_criteria=0; -iter=0; -offOrigin=geo.offOrigin; -offDetector=geo.offDetector; -rotDetector=geo.rotDetector; -DSD=geo.DSD; -DSO=geo.DSO; -while ~stop_criteria %POCS - % If quality is going to be measured, then we need to save previous image - if measurequality && ~strcmp(QualMeasOpts,'error_norm') - res_prev = f; % only store if necessary - end - - f0=f; - if (iter==0 && verbose==1);tic;end - iter=iter+1; - - for jj=1:size(angles,2) - if size(offOrigin,2)==size(angles,2) - geo.offOrigin=offOrigin(:,index_angles(:,jj)); - end - if size(offDetector,2)==size(angles,2) - geo.offDetector=offDetector(:,index_angles(:,jj)); - end - if size(rotDetector,2)==size(angles,2) - geo.rotDetector=rotDetector(:,index_angles(:,jj)); - end - if size(DSD,2)==size(angles,2) - geo.DSD=DSD(jj); - end - if size(DSO,2)==size(angles,2) - geo.DSO=DSO(jj); - end - % proj_err=proj(:,:,jj)-Ax(f,geo,angles(:,jj)); % (b-Ax) - % weighted_err=W(:,:,jj).*proj_err; % W^-1 * (b-Ax) - % backprj=Atb(weighted_err,geo,angles(:,jj)); % At * W^-1 * (b-Ax) - % weight_backprj=bsxfun(@times,1./V(:,:,jj),backprj); % V * At * W^-1 * (b-Ax) - % f=f+beta*weight_backprj; % x= x + lambda * V * At * W^-1 * (b-Ax) - % Enforce positivity - f=f+beta* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(f,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids)); - % non-negativity constrain - if nonneg - f=max(f,0); - end - end - - geo.offDetector=offDetector; - geo.offOrigin=offOrigin; - geo.DSD=DSD; - geo.DSO=DSO; - geo.rotDetector=rotDetector; - % Save copy of image. - if measurequality - qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); - end - % compute L2 error of actual image. Ax-b - dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2'); - % compute change in the image after last SART iteration - dp_vec=(f-f0); - dp=im3Dnorm(dp_vec,'L2'); - - if iter==1 - dtvg=alpha*dp; - %Convert the steepest-descent step-size from a fraction of a - %step-size to an absolute image distance on the first iteration. - end - f0=f; - - % TV MINIMIZATION - % ========================================================================= - % Call GPU to minimize TV - f=minimizeTV(f0,dtvg,ng,'gpuids',gpuids); % This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays - % for ii=1:ng - % % Steepest descend of TV norm - % tv(ng*(iter-1)+ii)=im3Dnorm(f,'TV','forward'); - % df=weighted_gradientTVnorm2(f,0.002); - % df=df./im3Dnorm(df,'L2'); - % f=f-dtvg.*df; - % end - - % update parameters - % ========================================================================== - - % compute change by TV min - dg_vec=(f-f0); - dg=im3Dnorm(dg_vec,'L2'); - % if change in TV is bigger than the change in SART AND image error is still bigger than acceptable - if dg>rmax*dp && dd>epsilon - dtvg=dtvg*alpha_red; - end - % reduce SART step - beta=beta*beta_red; - % Check convergence criteria - % ========================================================================== - - %Define c_alpha as in equation 21 in the journal - c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2)); - %This c is examined to see if it is close to -1.0 - - if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>=maxiter - if verbose - disp(['Stopping criteria met']); - disp([' c = ' num2str(c), '(Desired: c<-0.99)']); - disp([' beta = ' num2str(beta), '(Desired: beta<0.005)']); - disp([' iter = ' num2str(iter), ]); - end - stop_criteria=true; - end - if (iter==1 && verbose==1) - expected_time=toc*maxiter; - disp('ADS_POCS'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end - -end - - - -end - -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin) - -opts= {'lambda','lambda_red','init','tviter','verbose','alpha','alpha_red','ratio','maxl2err','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'}; -defaults=ones(length(opts),1); -% Check inputs -nVarargs = length(argin); -if mod(nVarargs,2) - error('TIGRE:ASD_POCS:InvalidInput','Invalid number of inputs') -end - -% check if option has been passed as input -for ii=1:2:nVarargs - ind=find(ismember(opts,lower(argin{ii}))); - if ~isempty(ind) - defaults(ind)=0; - else - error('TIGRE:ASD_POCS:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); - end -end - -for ii=1:length(opts) - opt=opts{ii}; - default=defaults(ii); - % if one option isnot default, then extract value from input - if default==0 - ind=double.empty(0,1);jj=1; - while isempty(ind) - ind=find(isequal(opt,lower(argin{jj}))); - jj=jj+1; - end - if isempty(ind) - error('TIGRE:ASD_POCS:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); - end - val=argin{jj}; - end - % parse inputs - switch opt - % Verbose - % ========================================================================= - case 'verbose' - if default - verbose=1; - else - verbose=val; - end - if ~is2014bOrNewer - warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); - verbose=false; - end - % Lambda - % ========================================================================= - % Its called beta in ASD-POCS - case 'lambda' - if default - beta=1; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:ASD_POCS:InvalidInput','Invalid lambda') - end - beta=val; - end - % Lambda reduction - % ========================================================================= - case 'lambda_red' - if default - beta_red=0.99; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:ASD_POCS:InvalidInput','Invalid lambda') - end - beta_red=val; - end - % Initial image - % ========================================================================= - case 'init' - if default || strcmp(val,'none') - f0=zeros(geo.nVoxel','single'); - - else - if strcmp(val,'FDK') - f0=FDK(proj, geo, angles); - else - error('TIGRE:ASD_POCS:InvalidInput','Invalid init') - - end - end - % Number of iterations of TV - % ========================================================================= - case 'tviter' - if default - ng=20; - else - ng=val; - end - % TV hyperparameter - % ========================================================================= - case 'alpha' - if default - alpha=0.002; % 0.2 - else - alpha=val; - end - % TV hyperparameter redution - % ========================================================================= - case 'alpha_red' - if default - alpha_red=0.95; - else - alpha_red=val; - end - % Maximum update ratio - % ========================================================================= - case 'ratio' - if default - rmax=0.95; - else - rmax=val; - end - % Maximum L2 error to have a "good image" - % ========================================================================= - case 'maxl2err' - if default - epsilon=im3Dnorm(Ax(FDK(proj,geo,angles),geo,angles)-proj,'L2')*0.2; %heuristic - else - epsilon=val; - end - % Order strategy - % ========================================================================= - case 'orderstrategy' - if default - OrderStrategy='random'; - else - OrderStrategy=val; - end - % Image Quality Measure - % ========================================================================= - case 'qualmeas' - if default - QualMeasOpts={}; - else - if iscellstr(val) - QualMeasOpts=val; - else - error('TIGRE:ASD_POCS:InvalidInput','Invalid quality measurement parameters'); - end - end - % Non negative - % ========================================================================= - case 'nonneg' - if default - nonneg=true; - else - nonneg=val; - end - % GPU Ids - % ========================================================================= - case 'gpuids' - if default - gpuids = GpuIds(); - else - gpuids = val; - end - % Data Redundancy weighting - case 'redundancy_weighting' - if default - redundancy_weights = true; - else - redundancy_weights = val; - end - case 'groundtruth' - if default - gt=nan; - else - gt=val; - end - otherwise - error('TIGRE:ASD_POCS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in ASD_POCS()']); - - end -end - -end - - +function [ f,qualMeasOut ] = ASD_POCS(proj,geo,angles,maxiter,varargin) +%ASD_POCS Solves the ASD_POCS total variation constrained image in 3D +% tomography. +% +% ASD_POCS(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ALPHA angles, corresponding +% to the geometry described in GEO, using NITER iterations. +% +% ASD_POCS(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'lambda': Sets the value of the hyperparameter for the SART iterations. +% Default is 1 +% +% 'lambdared': Reduction of lambda. Every iteration +% lambda = lambdared * lambda. Default is 0.99 +% +% 'init': Describes different initialization techniques. +% • 'none' : Initializes the image to zeros (default) +% +% • 'FDK' : Initializes image to FDK reconstruction +% 'TViter': Defines the amount of TV iterations performed per SART +% iteration. Default is 20 +% +% 'alpha': Defines the TV hyperparameter. Default is 0.002 +% +% 'alpha_red': Defines the reduction rate of the TV hyperparameter +% +% 'Ratio': The maximum allowed image/TV update ratio. If the TV +% update changes the image more than this, the parameter +% will be reduced. Default is 0.95 +% 'maxL2err' Maximum L2 error to accept an image as valid. This +% parameter is crucial for the algorithm, determines at +% what point an image should not be updated further. +% Default is 20% of the FDK L2 norm. +% 'Verbose' 1 or 0. Default is 1. Gives information about the +% progress of the algorithm. +% +% 'OrderStrategy' Chooses the subset ordering strategy. Options are +% 'ordered' : uses them in the input order, but divided +% 'random' : orders them randomly +% 'angularDistance': chooses the next subset with the +% biggest angular distance with the ones used. +% 'redundancy_weighting': true or false. Default is true. Applies data +% redundancy weighting to projections in the update step +% (relevant for offset detector geometry) +% 'groundTruth' an image as ground truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the TIGRE Toolbox +% +% Copyright (c) 2015, University of Bath and +% CERN-European Organization for Nuclear Research +% All rights reserved. +% +% License: Open Source under BSD. +% See the full license at +% https://github.com/CERN/TIGRE/blob/master/LICENSE +% +% Contact: tigre.toolbox@gmail.com +% Codes: https://github.com/CERN/TIGRE/ +% Coded by: Ander Biguri and Manasavee Lohvithee +%-------------------------------------------------------------------------- + + + +%% parse inputs + if nargin < 4 + error('TIGRE:ASD_POCS:InvalidInput','Not enough inputs') + end + [ f,qualMeasOut ] = OS_ASD_POCS(proj,geo,angles,maxiter, varargin{:}, 'BlockSize', 1); +end diff --git a/MATLAB/Algorithms/AwASD_POCS.m b/MATLAB/Algorithms/AwASD_POCS.m index 6cac5e24..b8094a9a 100644 --- a/MATLAB/Algorithms/AwASD_POCS.m +++ b/MATLAB/Algorithms/AwASD_POCS.m @@ -1,423 +1,80 @@ -function [ f,qualMeasOut]= AwASD_POCS(proj,geo,angles,maxiter,varargin) -%AwASD_POCS Solves the 3D tomography problem using the adaptive-weighted -%ASD_POCS algorithm which extends from the method ASD_POCS available in the -%TIGRE toolbox by adding weight equation to better preserve the edge of the -%reconstructed image -% -% AwASD_POCS(PROJ,GEO,ALPHA,NITER) added adaptive-weighted solves the reconstruction problem -% using the projection data PROJ taken over ALPHA angles, corresponding -% to the geometry described in GEO, using NITER iterations. -% -% AwASD_POCS(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The -% possible options in OPT are: -% -% -% 'lambda': Sets the value of the hyperparameter for the SART iterations. -% Default is 1 -% -% 'lambdared': Reduction of lambda.Every iteration -% lambda=lambdared*lambda. Default is 0.99 -% -% 'init': Describes different initialization techniques. -% • 'none' : Initializes the image to zeros (default) - -% • 'FDK' : initializes image to FDK reconstruction -% -% 'TViter': Defines the amount of TV iterations performed per SART -% iteration. Default is 20 -% -% 'alpha': Defines the TV hyperparameter. Default is 0.002 -% -% 'alpha_red': Defines the reduction rate of the TV hyperparameter -% -% 'Ratio': The maximum allowed image/TV update ratio. If the TV -% update changes the image more than this, the parameter -% will be reduced. Default is 0.95 -% 'maxL2err' Maximum L2 error to accept an image as valid. This -% parameter is crucial for the algorithm, determines at -% what point an image should not be updated further. -% Default is 20% of the FDK L2 norm. -% 'Verbose' 1 or 0. Default is 1. Gives information about the -% progress of the algorithm. -% 'delta' Defines Parameter to control the amount of smoothing -% for pixels at the edges. A large 'delta' is not able to -% differentiate image gradients at different pixels. A -% small 'delta' give low weights to almost every pixels, -% making the algorithm inefficient in removing noise or -% straking artifacts. Default is -0.00055 -% 'OrderStrategy' Chooses the subset ordering strategy. Options are -% 'ordered' :uses them in the input order, but divided -% 'random' : orders them randomly -% 'angularDistance': chooses the next subset with the -% biggest angular distance with the ones used. -% 'redundancy_weighting': true or false. Default is true. Applies data -% redundancy weighting to projections in the update step -% (relevant for offset detector geometry) -% 'groundTruth' an image as ground truth, to be used if quality measures -% are requested, to plot their change w.r.t. this known -% data. -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and -% CERN-European Organization for Nuclear Research -% All rights reserved. -% -% License: Open Source under BSD. -% See the full license at -% https://github.com/CERN/TIGRE/blob/master/LICENSE -% -% Contact: tigre.toolbox@gmail.com -% Codes: https://github.com/CERN/TIGRE/ -% Coded by: Ander Biguri and Manasavee Lohvithee - -%% parse inputs -blocksize=1; -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); - -measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); -if ~any(isnan(gt(:))) - QualMeasOpts{end+1}='error_norm'; - res_prev=gt; - clear gt -end -if nargout<2 && measurequality - warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them") - measurequality=false; -end -qualMeasOut=zeros(length(QualMeasOpts),niter); - - -[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); - -angles=cell2mat(alphablocks); -index_angles=cell2mat(orig_index); - -% does detector rotation exists? -if ~isfield(geo,'rotDetector') - geo.rotDetector=[0;0;0]; -end - -%% Create weighting matrices for the SART step -% the reason we do this, instead of calling the SART function is not to -% recompute the weights every AwASD-POCS iteration, thus effectively doubling -% the computational time -% Projection weight, W -W=computeW(geo,angles,gpuids); - - -% Back-Projection weight, V -V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids); - -if redundancy_weights - % Data redundancy weighting, W_r implemented using Wang weighting - % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc - - num_frames = size(proj,3); - W_r = redundancy_weighting(geo); - W_r = repmat(W_r,[1,1,num_frames]); - % disp('Size of redundancy weighting matrix'); - % disp(size(W_r)); - W = W.*W_r; % include redundancy weighting in W -end - -stop_criteria=0; -iter=0; -offOrigin=geo.offOrigin; -offDetector=geo.offDetector; -rotDetector=geo.rotDetector; -DSD=geo.DSD; -DSO=geo.DSO; - -while ~stop_criteria %POCS - % If quality is going to be measured, then we need to save previous image - if measurequality && ~strcmp(QualMeasOpts,'error_norm') - res_prev = f; % only store if necessary - end - f0=f; - if (iter==0 && verbose==1);tic;end - iter=iter+1; - - for jj=1:size(angles,2) - if size(offOrigin,2)==size(angles,2) - geo.offOrigin=offOrigin(:,index_angles(:,jj)); - end - if size(offDetector,2)==size(angles,2) - geo.offDetector=offDetector(:,index_angles(:,jj)); - end - if size(rotDetector,2)==size(angles,2) - geo.rotDetector=rotDetector(:,index_angles(:,jj)); - end - if size(DSD,2)==size(angles,2) - geo.DSD=DSD(jj); - end - if size(DSO,2)==size(angles,2) - geo.DSO=DSO(jj); - end - f=f+beta* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,jj).*(proj(:,:,index_angles(:,jj))-Ax(f,geo,angles(:,jj),'gpuids',gpuids)),geo,angles(:,jj),'gpuids',gpuids)); - if nonneg - f(f<0)=0; - end - end - - geo.offDetector=offDetector; - geo.offOrigin=offOrigin; - geo.DSD=DSD; - geo.DSO=DSO; - geo.rotDetector=rotDetector; - if measurequality - qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); - end - - % compute L2 error of actual image. Ax-b - dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2'); - % compute change in the image after last SART iteration - dp_vec=(f-f0); - dp=im3Dnorm(dp_vec,'L2'); - - if iter==1 - dtvg=alpha*dp; - %Convert the steepest-descent step-size from a fraction of a - %step-size to an absolute image distance on the first iteration. - end - f0=f; - - % TV MINIMIZATION - % ========================================================================= - % Call GPU to minimize TV - f=minimizeAwTV(f0,dtvg,ng,delta,'gpuids',gpuids); % This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays - % for ii=1:ng - % % Steepest descend of TV norm - % tv(ng*(iter-1)+ii)=im3Dnorm(f,'TV','forward'); - % df=weighted_gradientTVnorm(f,delta); - % df=df./im3Dnorm(df,'L2'); - % f=f-dtvg.*df; - % end - - % update parameters - % ========================================================================== - - % compute change by TV min - dg_vec=(f-f0); - dg=im3Dnorm(dg_vec,'L2'); - % if change in TV is bigger than the change in SART AND image error is still bigger than acceptable - if dg>rmax*dp && dd>epsilon - dtvg=dtvg*alpha_red; - end - % reduce SART step - beta=beta*beta_red; - % Check convergence criteria - % ========================================================================== - - %Define c_alpha as in equation 21 in the journal - c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2)); - %This c is examined to see if it is close to -1.0 - - if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>maxiter - if verbose - disp(['Stopping criteria met']); - disp([' c = ' num2str(c)]); - disp([' beta = ' num2str(beta)]); - disp([' iter = ' num2str(iter)]); - end - stop_criteria=true; - end - - if (iter==1 && verbose==1) - expected_time=toc*maxiter; - disp('AwASD_POCS'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end - -end - -end - -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','alpha','alpha_red','ratio','maxl2err','delta','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'}; -defaults=ones(length(opts),1); -% Check inputs -nVarargs = length(argin); -if mod(nVarargs,2) - error('TIGRE:AwASD_POCS:InvalidInput','Invalid number of inputs') -end - -% check if option has been passed as input -for ii=1:2:nVarargs - ind=find(ismember(opts,lower(argin{ii}))); - if ~isempty(ind) - defaults(ind)=0; - else - error('TIGRE:AwASD_POCS:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); - end -end - -for ii=1:length(opts) - opt=opts{ii}; - default=defaults(ii); - % if one option isnot default, then extract value from input - if default==0 - ind=double.empty(0,1);jj=1; - while isempty(ind) - ind=find(isequal(opt,lower(argin{jj}))); - jj=jj+1; - end - if isempty(ind) - error('TIGRE:AwASD_POCS:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); - end - val=argin{jj}; - end - % parse inputs - switch opt - % Verbose - % ========================================================================= - case 'verbose' - if default - verbose=1; - else - verbose=val; - end - if ~is2014bOrNewer - warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); - verbose=false; - end - % Lambda - % ========================================================================= - case 'lambda' - if default - beta=1; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:AwASD_POCS:InvalidInput','Invalid lambda') - end - beta=val; - end - % Lambda reduction - % ========================================================================= - case 'lambda_red' - if default - beta_red=0.99; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:AwASD_POCS:InvalidInput','Invalid lambda') - end - beta_red=val; - end - % Initial image - % ========================================================================= - case 'init' - if default || strcmp(val,'none') - f0=zeros(geo.nVoxel','single'); - - else - if strcmp(val,'FDK') - f0=FDK(proj, geo, angles); - else - error('TIGRE:AwASD_POCS:InvalidInput','Invalid init') - end - end - % Number of iterations of TV - % ========================================================================= - case 'tviter' - if default - ng=20; - else - ng=val; - end - % TV hyperparameter - % ========================================================================= - case 'alpha' - if default - alpha=0.002; % 0.2 - else - alpha=val; - end - % TV hyperparameter redution - % ========================================================================= - case 'alpha_red' - if default - alpha_red=0.95; - else - alpha_red=val; - end - % Maximum update ratio - % ========================================================================= - case 'ratio' - if default - rmax=0.95; - else - rmax=val; - end - % Maximum L2 error to have a "good image" - % ========================================================================= - case 'maxl2err' - if default - epsilon=im3Dnorm(FDK(proj,geo,angles),'L2')*0.2; %heuristic - else - epsilon=val; - end - %Parameter to control the amount of smoothing for pixels at the - %edges - % ========================================================================= - case 'delta' - if default - delta=-0.005; - else - delta=val; - end - case 'orderstrategy' - if default - OrderStrategy='random'; - else - OrderStrategy=val; - end - % Image Quality Measure - % ========================================================================= - case 'qualmeas' - if default - QualMeasOpts={}; - else - if iscellstr(val) - QualMeasOpts=val; - else - error('TIGRE:AwASD_POCS:InvalidInput','Invalid quality measurement parameters'); - end - end - % Non negative - % ========================================================================= - case 'nonneg' - if default - nonneg=true; - else - nonneg=val; - end - % GPU Ids - % ========================================================================= - case 'gpuids' - if default - gpuids = GpuIds(); - else - gpuids = val; - end - case 'redundancy_weighting' - if default - redundancy_weights = true; - else - redundancy_weights = val; - end - case 'groundtruth' - if default - gt=nan; - else - gt=val; - end - otherwise - error('TIGRE:AwASD_POCS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in AwASD_POCS()']); - - end -end - -end +function [ f,qualMeasOut]= AwASD_POCS(proj,geo,angles,maxiter,varargin) +%AwASD_POCS Solves the 3D tomography problem using the adaptive-weighted +%ASD_POCS algorithm which extends from the method ASD_POCS available in the +%TIGRE toolbox by adding weight equation to better preserve the edge of the +%reconstructed image +% +% AwASD_POCS(PROJ,GEO,ALPHA,NITER) added adaptive-weighted solves the reconstruction problem +% using the projection data PROJ taken over ALPHA angles, corresponding +% to the geometry described in GEO, using NITER iterations. +% +% AwASD_POCS(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'lambda': Sets the value of the hyperparameter for the SART iterations. +% Default is 1 +% +% 'lambdared': Reduction of lambda.Every iteration +% lambda=lambdared*lambda. Default is 0.99 +% +% 'init': Describes different initialization techniques. +% • 'none' : Initializes the image to zeros (default) + +% • 'FDK' : initializes image to FDK reconstruction +% +% 'TViter': Defines the amount of TV iterations performed per SART +% iteration. Default is 20 +% +% 'alpha': Defines the TV hyperparameter. Default is 0.002 +% +% 'alpha_red': Defines the reduction rate of the TV hyperparameter +% +% 'Ratio': The maximum allowed image/TV update ratio. If the TV +% update changes the image more than this, the parameter +% will be reduced. Default is 0.95 +% 'maxL2err' Maximum L2 error to accept an image as valid. This +% parameter is crucial for the algorithm, determines at +% what point an image should not be updated further. +% Default is 20% of the FDK L2 norm. +% 'Verbose' 1 or 0. Default is 1. Gives information about the +% progress of the algorithm. +% 'delta' Defines Parameter to control the amount of smoothing +% for pixels at the edges. A large 'delta' is not able to +% differentiate image gradients at different pixels. A +% small 'delta' give low weights to almost every pixels, +% making the algorithm inefficient in removing noise or +% straking artifacts. Default is -0.00055 +% 'OrderStrategy' Chooses the subset ordering strategy. Options are +% 'ordered' :uses them in the input order, but divided +% 'random' : orders them randomly +% 'angularDistance': chooses the next subset with the +% biggest angular distance with the ones used. +% 'redundancy_weighting': true or false. Default is true. Applies data +% redundancy weighting to projections in the update step +% (relevant for offset detector geometry) +% 'groundTruth' an image as ground truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the TIGRE Toolbox +% +% Copyright (c) 2015, University of Bath and +% CERN-European Organization for Nuclear Research +% All rights reserved. +% +% License: Open Source under BSD. +% See the full license at +% https://github.com/CERN/TIGRE/blob/master/LICENSE +% +% Contact: tigre.toolbox@gmail.com +% Codes: https://github.com/CERN/TIGRE/ +% Coded by: Ander Biguri and Manasavee Lohvithee + +%% parse inputs + if nargin < 4 + error('TIGRE:AwASD_POCS:InvalidInput','Not enough inputs') + end + [ f,qualMeasOut] = OS_AwASD_POCS(proj,geo,angles,maxiter, varargin{:}, 'BlockSize', 1); +end diff --git a/MATLAB/Algorithms/AwPCSD.m b/MATLAB/Algorithms/AwPCSD.m index cd5d01a1..3f3054ea 100644 --- a/MATLAB/Algorithms/AwPCSD.m +++ b/MATLAB/Algorithms/AwPCSD.m @@ -1,407 +1,66 @@ -function [ f,resL2,qualMeasOut] = AwPCSD(proj,geo,angles,maxiter,varargin) -%AwPCSD solves the reconstruction problem using adaptive-weighted -%projection-controlled steepest descent method -% -% AwPCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using -% the projection data PROJ taken over ALPHA angles, corresponding to the -% geometry described in GEO, using NITER iterations. -% -% AwPCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The -% possible options in OPT are: -% -% -% 'lambda': Sets the value of the hyperparameter for the SART iterations. -% Default is 1 -% -% 'lambdared': Reduction of lambda.Every iteration -% lambda=lambdared*lambda. Default is 0.99 -% -% 'init': Describes different initialization techniques. -% • 'none' : Initializes the image to zeros(default) - -% • 'FDK' : initializes image to FDK reconstruction -% -% 'TViter': Defines the amount of TV iterations performed per SART -% iteration. Default is 20 -% -% 'maxL2err' Maximum L2 error to accept an image as valid. This -% parameter is crucial for the algorithm, determines at -% what point an image should not be updated further. -% Default is 20% of the FDK L2 norm. -% 'Verbose' 1 or 0. Default is 1. Gives information about the -% progress of the algorithm. -% 'delta' Defines Parameter to control the amount of smoothing -% for pixels at the edges. A large 'delta' is not able to -% differentiate image gradients at different pixels. A -% small 'delta' give low weights to almost every pixels, -% making the algorithm inefficient in removing noise or -% straking artifacts. Default is -0.00055 -% 'redundancy_weighting': true or false. Default is true. Applies data -% redundancy weighting to projections in the update step -% (relevant for offset detector geometry) -% 'groundTruth' an image as ground truth, to be used if quality measures -% are requested, to plot their change w.r.t. this known -% data. -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and -% CERN-European Organization for Nuclear Research -% All rights reserved. -% -% License: Open Source under BSD. -% See the full license at -% https://github.com/CERN/TIGRE/blob/master/LICENSE -% -% Contact: tigre.toolbox@gmail.com -% Codes: https://github.com/CERN/TIGRE/ -% Coded by: Ander Biguri and Manasavee Lohvithee -%-------------------------------------------------------------------------- -%% parse inputs -[beta,beta_red,f,ng,verbose,epsilon,delta,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); -if ~any(isnan(gt(:))) - QualMeasOpts{end+1}='error_norm'; - res_prev=gt; - clear gt -end -if nargout<3 && measurequality - warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them") - measurequality=false; -end -qualMeasOut=zeros(length(QualMeasOpts),niter); - -if nargout>1 - computeL2=true; -else - computeL2=false; -end -resL2=zeros(1,niter); - - - - -% [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); -% -% angles=cell2mat(alphablocks); -% index_angles=cell2mat(orig_index); - - -% does detector rotation exists? -if ~isfield(geo,'rotDetector') - geo.rotDetector=[0;0;0]; -end - -%% Create weighting matrices for the SART step -% the reason we do this, instead of calling the SART function is not to -% recompute the weights every AwASD-POCS iteration, thus effectively doubling -% the computational time -% Projection weight, W -W=computeW(geo,angles,gpuids); - -% Back-Projection weight, V -V=computeV(geo,angles,num2cell(angles),num2cell(1:length(angles)),'gpuids',gpuids); - -if redundancy_weights - % Data redundancy weighting, W_r implemented using Wang weighting - % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc - - num_frames = size(proj,3); - W_r = redundancy_weighting(geo); - W_r = repmat(W_r,[1,1,num_frames]); - % disp('Size of redundancy weighting matrix'); - % disp(size(W_r)); - W = W.*W_r; % include redundancy weighting in W -end -%Initialize image. -%f=zeros(geo.nVoxel','single'); - -iter=0; -offOrigin=geo.offOrigin; -offDetector=geo.offDetector; -rotDetector=geo.rotDetector; -stop_criteria=0; -DSD=geo.DSD; -DSO=geo.DSO; -while ~stop_criteria %POCS - % If quality is going to be measured, then we need to save previous image - if measurequality && ~strcmp(QualMeasOpts,'error_norm') - res_prev = f; % only store if necessary - end - f0=f; - if (iter==0 && verbose==1);tic;end - iter=iter+1; - - %Estimation error in the projection domain - est_proj=Ax(f,geo,angles,'interpolated'); - delta_p=im3Dnorm(est_proj-proj,'L2'); - - %Enforcing ART along all projections if squared delta_p > epsilon - if (delta_p^2)>epsilon - for jj=1:size(angles,2) - if size(offOrigin,2)==size(angles,2) - geo.offOrigin=offOrigin(:,jj); - end - if size(offDetector,2)==size(angles,2) - geo.offDetector=offDetector(:,jj); - end - if size(rotDetector,2)==size(angles,2) - geo.rotDetector=rotDetector(:,jj); - end - if size(DSD,2)==size(angles,2) - geo.DSD=DSD(jj); - end - if size(DSO,2)==size(angles,2) - geo.DSO=DSO(jj); - end - - f=f+beta* bsxfun(@times,1./sum(V(:,:,jj),3),Atb(W(:,:,jj).*(proj(:,:,jj)-Ax(f,geo,angles(:,jj),'gpuids',gpuids)),geo,angles(:,jj),'gpuids',gpuids)); - - end - end - - %Non-negativity projection on all pixels - if nonneg - f=max(f,0); - end - - geo.offDetector=offDetector; - geo.offOrigin=offOrigin; - - if measurequality - qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); - end - - % Compute L2 error of actual image. Ax-b - dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2'); - - % Compute change in the image after last SART iteration - dp_vec=(f-f0); - - if iter==1 - step=1; - else - step=delta_p/delta_p_first; - end - f0=f; - % TV MINIMIZATION - % ========================================================================= - % Call GPU to minimize AwTV - f=minimizeAwTV(f0,step,ng,delta,'gpuids',gpuids); % This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays - % for ii=1:ng - % %delta=-0.00038 for thorax phantom - % df=weighted_gradientTVnorm(f,delta); - % df=df./im3Dnorm(df,'L2'); - % f=f-(step.*df); - % end - % Compute change by TV min - dg_vec=(f-f0); - - if iter==1 - delta_p_first=im3Dnorm((Ax(f0,geo,angles,'interpolated','gpuids',gpuids))-proj,'L2'); - end - - % Reduce SART step - beta=beta*beta_red; - - - % Check convergence criteria - % ========================================================================== - - %Define c_alpha as in equation 21 in the journal - c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2)); - %This c is examined to see if it is close to -1.0 - - if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>maxiter - if verbose - disp(['Stopping criteria met']); - disp([' c = ' num2str(c)]); - disp([' beta = ' num2str(beta)]); - disp([' iter = ' num2str(iter)]); - end - stop_criteria=true; - end - - if computeL2 - geo.offOrigin=offOrigin; - geo.offDetector=offDetector; - errornow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax - % If the error is not minimized. - if iter~=1 && errornow>errorL2(end) - if verbose - disp(['Convergence criteria met, exiting on iteration number:', num2str(iter)]); - end - return; - end - errorL2=[errorL2 errornow]; - end - - - if (iter==1 && verbose==1) - expected_time=toc*maxiter; - disp('AwPCSD'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end - -end -end - - -function [beta,beta_red,f0,ng,verbose,epsilon,delta,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','maxl2err','delta','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'}; -defaults=ones(length(opts),1); -% Check inputs -nVarargs = length(argin); -if mod(nVarargs,2) - error('TIGRE:AwPCSD:InvalidInput','Invalid number of inputs') -end - -% check if option has been passed as input -for ii=1:2:nVarargs - ind=find(ismember(opts,lower(argin{ii}))); - if ~isempty(ind) - defaults(ind)=0; - else - error('TIGRE:AwPCSD:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); - end -end - -for ii=1:length(opts) - opt=opts{ii}; - default=defaults(ii); - % if one option isnot default, then extract value from input - if default==0 - ind=double.empty(0,1);jj=1; - while isempty(ind) - ind=find(isequal(opt,lower(argin{jj}))); - jj=jj+1; - end - if isempty(ind) - error('TIGRE:AwPCSD:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); - end - val=argin{jj}; - end - % parse inputs - switch opt - % Verbose - % ========================================================================= - case 'verbose' - if default - verbose=1; - else - verbose=val; - end - if ~is2014bOrNewer - warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); - verbose=false; - end - % Lambda - % ========================================================================= - case 'lambda' - if default - beta=1; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:AwPCSD:InvalidInput','Invalid lambda') - end - beta=val; - end - % Lambda reduction - % ========================================================================= - case 'lambda_red' - if default - beta_red=0.99; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:AwPCSD:InvalidInput','Invalid lambda') - end - beta_red=val; - end - % Initial image - % ========================================================================= - case 'init' - if default || strcmp(val,'none') - f0=zeros(geo.nVoxel','single'); - - else - if strcmp(val,'FDK') - f0=FDK(proj, geo, angles); - else - error('TIGRE:AwPCSD:InvalidInput','Invalid init') - end - end - % Number of iterations of TV - % ========================================================================= - case 'tviter' - if default - ng=20; - else - ng=val; - end - % Maximum L2 error to have a "good image" - % ========================================================================= - case 'maxl2err' - if default - epsilon=im3Dnorm(FDK(proj,geo,angles))*0.2; %heuristic - else - epsilon=val; - end - %Parameter to control the amount of smoothing for pixels at the - %edges - % ========================================================================= - case 'delta' - if default - delta=-0.005; - else - delta=val; - end - % Image Quality Measure - % ========================================================================= - case 'qualmeas' - if default - QualMeasOpts={}; - else - if iscellstr(val) - QualMeasOpts=val; - else - error('TIGRE:AwPCSD:InvalidInput','Invalid quality measurement parameters'); - end - end - % Non negative - % ========================================================================= - case 'nonneg' - if default - nonneg=true; - else - nonneg=val; - end - % GPU Ids - % ========================================================================= - case 'gpuids' - if default - gpuids = GpuIds(); - else - gpuids = val; - end - case 'redundancy_weighting' - if default - redundancy_weights = true; - else - redundancy_weights = val; - end - case 'groundtruth' - if default - gt=nan; - else - gt=val; - end - otherwise - error('TIGRE:AwPCSD:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in AwPCSD()']); - - end -end - -end +function [ f,resL2,qualMeasOut] = AwPCSD(proj,geo,angles,maxiter,varargin) +%AwPCSD solves the reconstruction problem using adaptive-weighted +%projection-controlled steepest descent method +% +% AwPCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using +% the projection data PROJ taken over ALPHA angles, corresponding to the +% geometry described in GEO, using NITER iterations. +% +% AwPCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'lambda': Sets the value of the hyperparameter for the SART iterations. +% Default is 1 +% +% 'lambdared': Reduction of lambda.Every iteration +% lambda=lambdared*lambda. Default is 0.99 +% +% 'init': Describes different initialization techniques. +% • 'none' : Initializes the image to zeros(default) + +% • 'FDK' : initializes image to FDK reconstruction +% +% 'TViter': Defines the amount of TV iterations performed per SART +% iteration. Default is 20 +% +% 'maxL2err' Maximum L2 error to accept an image as valid. This +% parameter is crucial for the algorithm, determines at +% what point an image should not be updated further. +% Default is 20% of the FDK L2 norm. +% 'Verbose' 1 or 0. Default is 1. Gives information about the +% progress of the algorithm. +% 'delta' Defines Parameter to control the amount of smoothing +% for pixels at the edges. A large 'delta' is not able to +% differentiate image gradients at different pixels. A +% small 'delta' give low weights to almost every pixels, +% making the algorithm inefficient in removing noise or +% straking artifacts. Default is -0.00055 +% 'redundancy_weighting': true or false. Default is true. Applies data +% redundancy weighting to projections in the update step +% (relevant for offset detector geometry) +% 'groundTruth' an image as ground truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the TIGRE Toolbox +% +% Copyright (c) 2015, University of Bath and +% CERN-European Organization for Nuclear Research +% All rights reserved. +% +% License: Open Source under BSD. +% See the full license at +% https://github.com/CERN/TIGRE/blob/master/LICENSE +% +% Contact: tigre.toolbox@gmail.com +% Codes: https://github.com/CERN/TIGRE/ +% Coded by: Ander Biguri and Manasavee Lohvithee +%-------------------------------------------------------------------------- +%% parse inputs + if nargin < 4 + error('TIGRE:AwPCSD:InvalidInput','Not enough inputs') + end + [ f,resL2,qualMeasOut] = OS_AwPCSD(proj,geo,angles,maxiter, varargin{:}, 'BlockSize', 1); +end diff --git a/MATLAB/Algorithms/PCSD.m b/MATLAB/Algorithms/PCSD.m index c39dcc5b..373a4267 100644 --- a/MATLAB/Algorithms/PCSD.m +++ b/MATLAB/Algorithms/PCSD.m @@ -1,364 +1,59 @@ -function [ f,qualMeasOut] = PCSD(proj,geo,angles,maxiter,varargin) -%PCSD solves the reconstruction problem using projection-controlled steepest descent method -% -% PCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using -% the projection data PROJ taken over ALPHA angles, corresponding to the -% geometry described in GEO, using NITER iterations. -% -% PCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The -% possible options in OPT are: -% -% -% 'lambda': Sets the value of the hyperparameter for the SART iterations. -% Default is 1 -% -% 'lambdared': Reduction of lambda.Every iteration -% lambda=lambdared*lambda. Default is 0.99 -% -% 'init': Describes different initialization techniques. -% • 'none' : Initializes the image to zeros (default) -% • 'FDK' : initializes image to FDK reconstruction -% -% 'TViter': Defines the amount of TV iterations performed per SART -% iteration. Default is 20 -% -% 'maxL2err' Maximum L2 error to accept an image as valid. This -% parameter is crucial for the algorithm, determines at -% what point an image should not be updated further. -% Default is 20% of the FDK L2 norm. -% 'Verbose' 1 or 0. Default is 1. Gives information about the -% progress of the algorithm. -% 'redundancy_weighting': true or false. Default is true. Applies data -% redundancy weighting to projections in the update step -% (relevant for offset detector geometry) -% 'groundTruth' an image as ground truth, to be used if quality measures -% are requested, to plot their change w.r.t. this known -% data. -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and -% CERN-European Organization for Nuclear Research -% All rights reserved. -% -% License: Open Source under BSD. -% See the full license at -% https://github.com/CERN/TIGRE/blob/master/LICENSE -% -% Contact: tigre.toolbox@gmail.com -% Codes: https://github.com/CERN/TIGRE/ -% Coded by: Ander Biguri and Manasavee Lohvithee -%-------------------------------------------------------------------------- - -%% parse inputs -[beta,beta_red,f,ng,verbose,epsilon,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); - -measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); -if ~any(isnan(gt(:))) - QualMeasOpts{end+1}='error_norm'; - res_prev=gt; - clear gt -end -if nargout<2 && measurequality - warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them") - measurequality=false; -end -qualMeasOut=zeros(length(QualMeasOpts),niter); - -% does detector rotation exists? -if ~isfield(geo,'rotDetector') - geo.rotDetector=[0;0;0]; -end - -%% Create weighting matrices for the SART step -% the reason we do this, instead of calling the SART function is not to -% recompute the weights every AwASD-POCS iteration, thus effectively doubling -% the computational time -% Projection weight, W - -geoaux=geo; -geoaux.sVoxel([1 2])=geo.sVoxel([1 2])*1.1; % a Bit bigger, to avoid numerical division by zero (small number) -geoaux.sVoxel(3)=max(geo.sDetector(2),geo.sVoxel(3)); % make sure lines are not cropped. One is for when image is bigger than detector and viceversa -geoaux.nVoxel=[2,2,2]'; % accurate enough? -geoaux.dVoxel=geoaux.sVoxel./geoaux.nVoxel; -W=Ax(ones(geoaux.nVoxel','single'),geoaux,angles,'Siddon','gpuids',gpuids); -W(W epsilon - if (delta_p^2)>epsilon - for jj=1:size(angles,2) - if size(offOrigin,2)==size(angles,2) - geo.offOrigin=offOrigin(:,jj); - end - if size(offDetector,2)==size(angles,2) - geo.offDetector=offDetector(:,jj); - end - if size(rotDetector,2)==size(angles,2) - geo.rotDetector=rotDetector(:,jj); - end - if size(DSD,2)==size(angles,2) - geo.DSD=DSD(jj); - end - if size(DSO,2)==size(angles,2) - geo.DSO=DSO(jj); - end - f=f+beta* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,jj).*(proj(:,:,jj)-Ax(f,geo,angles(:,jj),'gpuids',gpuids)),geo,angles(:,jj),'gpuids',gpuids)); - - end - end - - %Non-negativity projection on all pixels - if nonneg - f=max(f,0); - end - - geo.offDetector=offDetector; - geo.offOrigin=offOrigin; - geo.DSD=DSD; - geo.DSO=DSO; - geo.rotDetector=rotDetector; - if measurequality - qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); - end - - % Compute L2 error of actual image. Ax-b - dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2'); - % Compute change in the image after last SART iteration - dp_vec=(f-f0); - - if iter==1 - step=1; - else - step=delta_p/delta_p_first; - end - f0=f; - % TV MINIMIZATION - % ========================================================================= - % Call GPU to minimize TV - f=minimizeTV(f0,step,ng,'gpuids',gpuids); % This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays - % for ii=1:ng - % %delta=-0.00038 for thorax phantom - % df=weighted_gradientTVnorm(f,delta); - % df=df./im3Dnorm(df,'L2'); - % f=f-(step.*df); - % end - - % Compute change by TV min - dg_vec=(f-f0); - - if iter==1 - delta_p_first=im3Dnorm((Ax(f0,geo,angles,'interpolated','gpuids',gpuids))-proj,'L2'); - end - - % Reduce SART step - beta=beta*beta_red; - - % Check convergence criteria - % ========================================================================== - - %Define c_alpha as in equation 21 in the journal - c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2)); - %This c is examined to see if it is close to -1.0 - - if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>maxiter - if verbose - disp(['Stopping criteria met']); - disp([' c = ' num2str(c)]); - disp([' beta = ' num2str(beta)]); - disp([' iter = ' num2str(iter)]); - end - stop_criteria=true; - end - - if (iter==1 && verbose==1) - expected_time=toc*maxiter; - disp('PCSD'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end - -end -end - - -function [beta,beta_red,f0,ng,verbose,epsilon,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','maxl2err','qualmeas','nonneg','gpuids','redundancy_weighting'}; -defaults=ones(length(opts),1); -% Check inputs -nVarargs = length(argin); -if mod(nVarargs,2) - error('CBCT:PCSD:InvalidInput','Invalid number of inputs') -end - -% check if option has been passed as input -for ii=1:2:nVarargs - ind=find(ismember(opts,lower(argin{ii}))); - if ~isempty(ind) - defaults(ind)=0; - else - error('CBCT:PCSD:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); - end -end - -for ii=1:length(opts) - opt=opts{ii}; - default=defaults(ii); - % if one option isnot default, then extract value from input - if default==0 - ind=double.empty(0,1);jj=1; - while isempty(ind) - ind=find(isequal(opt,lower(argin{jj}))); - jj=jj+1; - end - if isempty(ind) - error('CBCT:PCSD:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); - end - val=argin{jj}; - end - % parse inputs - switch opt - % Verbose - % ========================================================================= - case 'verbose' - if default - verbose=1; - else - verbose=val; - end - if ~is2014bOrNewer - warning('Verbose mode not available for older versions than MATLAB R2014b'); - verbose=false; - end - % Lambda - % ========================================================================= - case 'lambda' - if default - beta=1; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:PCSD:InvalidInput','Invalid lambda') - end - beta=val; - end - % Lambda reduction - % ========================================================================= - case 'lambda_red' - if default - beta_red=0.99; - else - if length(val)>1 || ~isnumeric( val) - error('TIGRE:PCSD:InvalidInput','Invalid lambda') - end - beta_red=val; - end - % Initial image - % ========================================================================= - case 'init' - if default || strcmp(val,'none') - f0=zeros(geo.nVoxel','single'); - - else - if strcmp(val,'FDK') - f0=FDK(proj, geo, angles); - else - error('TIGRE:PCSD:InvalidInput','Invalid init') - end - end - % Number of iterations of TV - % ========================================================================= - case 'tviter' - if default - ng=20; - else - ng=val; - end - % Maximum L2 error to have a "good image" - % ========================================================================= - case 'maxl2err' - if default - epsilon=im3Dnorm(FDK(proj,geo,angles))*0.2; %heuristic - else - epsilon=val; - end - %Image Quality Measure - % ========================================================================= - case 'qualmeas' - if default - QualMeasOpts={}; - else - if iscellstr(val) - QualMeasOpts=val; - else - error('TIGRE:PCSD:InvalidInput','Invalid quality measurement parameters'); - end - end - % Non negative - % ========================================================================= - case 'nonneg' - if default - nonneg=true; - else - nonneg=val; - end - % GPU Ids - % ========================================================================= - case 'gpuids' - if default - gpuids = GpuIds(); - else - gpuids = val; - end - case 'redundancy_weighting' - if default - redundancy_weights = true; - else - redundancy_weights = val; - end - otherwise - error('TIGRE:PCSD:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in PCSD()']); - - end -end - -end +function [ f,qualMeasOut] = PCSD(proj,geo,angles,maxiter,varargin) +%PCSD solves the reconstruction problem using projection-controlled steepest descent method +% +% PCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using +% the projection data PROJ taken over ALPHA angles, corresponding to the +% geometry described in GEO, using NITER iterations. +% +% PCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'lambda': Sets the value of the hyperparameter for the SART iterations. +% Default is 1 +% +% 'lambdared': Reduction of lambda.Every iteration +% lambda=lambdared*lambda. Default is 0.99 +% +% 'init': Describes different initialization techniques. +% • 'none' : Initializes the image to zeros (default) +% • 'FDK' : initializes image to FDK reconstruction +% +% 'TViter': Defines the amount of TV iterations performed per SART +% iteration. Default is 20 +% +% 'maxL2err' Maximum L2 error to accept an image as valid. This +% parameter is crucial for the algorithm, determines at +% what point an image should not be updated further. +% Default is 20% of the FDK L2 norm. +% 'Verbose' 1 or 0. Default is 1. Gives information about the +% progress of the algorithm. +% 'redundancy_weighting': true or false. Default is true. Applies data +% redundancy weighting to projections in the update step +% (relevant for offset detector geometry) +% 'groundTruth' an image as ground truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the TIGRE Toolbox +% +% Copyright (c) 2015, University of Bath and +% CERN-European Organization for Nuclear Research +% All rights reserved. +% +% License: Open Source under BSD. +% See the full license at +% https://github.com/CERN/TIGRE/blob/master/LICENSE +% +% Contact: tigre.toolbox@gmail.com +% Codes: https://github.com/CERN/TIGRE/ +% Coded by: Ander Biguri and Manasavee Lohvithee +%-------------------------------------------------------------------------- + +%% parse inputs + if nargin < 4 + error('TIGRE:PCSD:InvalidInput','Not enough inputs') + end + [ f,qualMeasOut] = OS_PCSD(proj,geo,angles,maxiter, varargin{:}, 'BlockSize', 1); +end diff --git a/MATLAB/Algorithms/SART.m b/MATLAB/Algorithms/SART.m index 9916fd20..3f60baf1 100644 --- a/MATLAB/Algorithms/SART.m +++ b/MATLAB/Algorithms/SART.m @@ -1,402 +1,68 @@ -function [res,resL2,qualMeasOut]=SART(proj,geo,angles,niter,varargin) -% SART solves Cone Beam CT image reconstruction using Oriented Subsets -% Simultaneous Algebraic Reconstruction Technique algorithm -% -% SART(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem -% using the projection data PROJ taken over ALPHA angles, corresponding -% to the geometry described in GEO, using NITER iterations. -% -% SART(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The -% possible options in OPT are: -% -% -% 'lambda': Sets the value of the hyperparameter. Default is 1 -% -% 'lambda_red': Reduction of lambda. Every iteration -% lambda=lambdared*lambda. Default is 0.99 -% -% 'Init': Describes different initialization techniques. -% 'none' : Initializes the image to zeros (default) -% 'FDK' : Initializes image to FDK reconstruction -% 'multigrid': Initializes image by solving the problem in -% small scale and increasing it when relative -% convergence is reached. -% 'image' : Initialization using a user specified -% image. Not recommended unless you really -% know what you are doing. -% 'InitImg' an image for the 'image' initialization. Avoid. -% -% 'Verbose' 1 or 0. Default is 1. Gives information about the -% progress of the algorithm. -% 'QualMeas' Asks the algorithm for a set of quality measurement -% parameters. Input should contain a cell array of desired -% quality measurement names. Example: {'CC','RMSE','MSSIM'} -% These will be computed in each iteration. -% 'OrderStrategy' Chooses the subset ordering strategy. Options are -% 'ordered' : uses them in the input order, but divided -% 'random' : orders them randomly -% 'angularDistance': chooses the next subset with the -% biggest angular distance with the ones used. -% 'redundancy_weighting': true or false. Default is true. Applies data -% redundancy weighting to projections in the update step -% (relevant for offset detector geometry) -% 'groundTruth' an image as ground truth, to be used if quality measures -% are requested, to plot their change w.r.t. this known -% data. -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and -% CERN-European Organization for Nuclear Research -% All rights reserved. -% -% License: Open Source under BSD. -% See the full license at -% https://github.com/CERN/TIGRE/blob/master/LICENSE -% -% Contact: tigre.toolbox@gmail.com -% Codes: https://github.com/CERN/TIGRE/ -% Coded by: Ander Biguri -%-------------------------------------------------------------------------- - -%% Deal with input parameters -blocksize=1; -[lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); - -measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); -if ~any(isnan(gt(:))) - QualMeasOpts{end+1}='error_norm'; - res_prev=gt; - clear gt -end -if nargout<3 && measurequality - warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them") - measurequality=false; -end -qualMeasOut=zeros(length(QualMeasOpts),niter); - -resL2=zeros(1,niter); -if nargout>1 - computeL2=true; -else - computeL2=false; -end - -[alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); -index_angles=cell2mat(orig_index); -angles_reorder=cell2mat(alphablocks); - -% does detector rotation exists? -if ~isfield(geo,'rotDetector') - geo.rotDetector=[0;0;0]; -end -%% Create weighting matrices - -% Projection weight, W - -% Projection weight, W -W=computeW(geo,angles,gpuids); - -% Back-Projection weight, V -V=computeV(geo,angles,alphablocks,orig_index,'gpuids',gpuids); - -if redundancy_weights - % Data redundancy weighting, W_r implemented using Wang weighting - % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc - - num_frames = size(proj,3); - W_r = redundancy_weighting(geo); - W_r = repmat(W_r,[1,1,num_frames]); - % disp('Size of redundancy weighting matrix'); - % disp(size(W_r)); - W = W.*W_r; % include redundancy weighting in W -end - -clear A x y dx dz; -%% hyperparameter stuff -nesterov=false; -if ischar(lambda)&&strcmp(lambda,'nesterov') - nesterov=true; - lambda=(1+sqrt(1+4))/2; - gamma=0; - ynesterov=zeros(size(res),'single'); - ynesterov_prev=ynesterov; -end -%% Iterate -offOrigin=geo.offOrigin; -offDetector=geo.offDetector; -rotDetector=geo.rotDetector; -DSD=geo.DSD; -DSO=geo.DSO; -% TODO : Add options for Stopping criteria -for ii=1:niter - if (ii==1 && verbose==1);tic;end - % If quality is going to be measured, then we need to save previous image - % THIS TAKES MEMORY! - if measurequality && ~strcmp(QualMeasOpts,'error_norm') - res_prev = res; % only store if necessary - end - - % reorder angles - - for jj=index_angles - if size(offOrigin,2)==size(angles,2) - geo.offOrigin=offOrigin(:,jj); - end - if size(offDetector,2)==size(angles,2) - geo.offDetector=offDetector(:,jj); - end - if size(rotDetector,2)==size(angles,2) - geo.rotDetector=rotDetector(:,jj); - end - if size(DSD,2)==size(angles,2) - geo.DSD=DSD(jj); - end - if size(DSO,2)==size(angles,2) - geo.DSO=DSO(jj); - end - % --------- Memory expensive----------- - - % proj_err=proj(:,:,jj)-Ax(res,geo,angles(jj)); % (b-Ax) - % weighted_err=W(:,:,jj).*proj_err; % W^-1 * (b-Ax) - % backprj=Atb(weighted_err,geo,angles(jj)); % At * W^-1 * (b-Ax) - % weight_backprj=bsxfun(@times,1./V(:,:,jj),backprj); % V * At * W^-1 * (b-Ax) - % res=res+lambda*weight_backprj; % x= x + lambda * V * At * W^-1 * (b-Ax) - %------------------------------------ - %--------- Memory cheap(er)----------- - if nesterov - % The nesterov update is quite similar to the normal update, it - % just uses this update, plus part of the last one. - ynesterov_prev = ynesterov; - ynesterov=res+ bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids)); - res=(1-gamma)*ynesterov+gamma*ynesterov_prev; - else - res=res+lambda* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids)); - end - if nonneg - res=max(res,0); - end - end - - % If quality is being measured - if measurequality - qualMeasOut(:,ii)=Measure_Quality(res,res_prev,QualMeasOpts); - end - - if nesterov - gamma=(1-lambda); - lambda=(1+sqrt(1+4*lambda^2))/2; - gamma=gamma/lambda; - else - lambda=lambda*lambdared; - end - if computeL2 || nesterov - geo.offOrigin=offOrigin; - geo.offDetector=offDetector; - geo.DSD=DSD; - geo.rotDetector=rotDetector; - resL2(ii)=im3Dnorm(proj-Ax(res,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax - % If the error is not minimized. - if ii~=1 && resL2(ii)>resL2(ii-1) - if verbose - disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]); - end - return - end - end - - if (ii==1 && verbose==1) - expected_time=toc*niter; - disp('SART'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end -end - - - - - -end - -function initres=init_multigrid(proj,geo,alpha) - -finalsize=geo.nVoxel; -% start with 64 -geo.nVoxel=[64;64;64]; -geo.dVoxel=geo.sVoxel./geo.nVoxel; -if any(finalsizefinalsize)=finalsize(geo.nVoxel>finalsize); - geo.dVoxel=geo.sVoxel./geo.nVoxel; - % Upsample! - % (hopefully computer has enough memory............) - [y, x, z]=ndgrid(linspace(1,size(initres,1),geo.nVoxel(1)),... - linspace(1,size(initres,2),geo.nVoxel(2)),... - linspace(1,size(initres,3),geo.nVoxel(3))); - initres=interp3(initres,x,y,z); - clear x y z -end -end - - -function [lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin) -opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','orderstrategy','nonneg','gpuids','redundancy_weighting','groundtruth'}; -defaults=ones(length(opts),1); -% Check inputs -nVarargs = length(argin); -if mod(nVarargs,2) - error('TIGRE:SART:InvalidInput','Invalid number of inputs') -end - -% check if option has been passed as input -for ii=1:2:nVarargs - ind=find(ismember(opts,lower(argin{ii}))); - if ~isempty(ind) - defaults(ind)=0; - else - error('TIGRE:SART:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); - end -end - -for ii=1:length(opts) - opt=opts{ii}; - default=defaults(ii); - % if one option is not default, then extract value from input - if default==0 - ind=double.empty(0,1);jj=1; - while isempty(ind) - ind=find(isequal(opt,lower(argin{jj}))); - jj=jj+1; - end - if isempty(ind) - error('TIGRE:SART:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); - end - val=argin{jj}; - end - - switch opt - % % % % % % % Verbose - case 'verbose' - if default - verbose=1; - else - verbose=val; - end - if ~is2014bOrNewer - warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b'); - verbose=false; - end - % % % % % % % hyperparameter, LAMBDA - case 'lambda' - if default - lambda=1; - elseif ischar(val)&&strcmpi(val,'nesterov') - lambda='nesterov'; % just for lowercase/uppercase - elseif length(val)>1 || ~isnumeric(val) - error('TIGRE:SART:InvalidInput','Invalid lambda') - else - lambda=val; - end - case 'lambda_red' - if default - lambdared=1; - else - if length(val)>1 || ~isnumeric(val) - error('TIGRE:SART:InvalidInput','Invalid lambda') - end - lambdared=val; - end - case 'init' - res=[]; - if default || strcmp(val,'none') - res=zeros(geo.nVoxel','single'); - continue - end - if strcmp(val,'FDK') - res=FDK(proj,geo,alpha); - continue - end - if strcmp(val,'multigrid') - res=init_multigrid(proj,geo,alpha); - continue - end - if strcmp(val,'image') - initwithimage=1; % it is used (10 lines below) - continue - end - if isempty(res) - error('TIGRE:SART:InvalidInput','Invalid Init option') - end - % % % % % % % ERROR - case 'initimg' - if default - continue - end - if exist('initwithimage','var') - if isequal(size(val),geo.nVoxel') - res=single(val); - else - error('TIGRE:SART:InvalidInput','Invalid image for initialization'); - end - end - case 'qualmeas' - if default - QualMeasOpts={}; - else - if iscellstr(val) - QualMeasOpts=val; - else - error('TIGRE:SART:InvalidInput','Invalid quality measurement parameters'); - end - end - case 'orderstrategy' - if default - OrderStrategy='random'; - else - OrderStrategy=val; - end - case 'nonneg' - if default - nonneg=true; - else - nonneg=val; - end - case 'gpuids' - if default - gpuids = GpuIds(); - else - gpuids = val; - end - case 'redundancy_weighting' - if default - redundancy_weights = true; - else - redundancy_weights = val; - end - case 'groundtruth' - if default - gt=nan; - else - gt=val; - end - otherwise - error('TIGRE:SART:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option']); - end -end - -end \ No newline at end of file +function [res,resL2,qualMeasOut]=SART(proj,geo,angles,niter,varargin) +% SART solves Cone Beam CT image reconstruction using Oriented Subsets +% Simultaneous Algebraic Reconstruction Technique algorithm +% +% SART(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ALPHA angles, corresponding +% to the geometry described in GEO, using NITER iterations. +% +% SART(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'lambda': Sets the value of the hyperparameter. Default is 1 +% +% 'lambda_red': Reduction of lambda. Every iteration +% lambda=lambdared*lambda. Default is 0.99 +% +% 'Init': Describes different initialization techniques. +% 'none' : Initializes the image to zeros (default) +% 'FDK' : Initializes image to FDK reconstruction +% 'multigrid': Initializes image by solving the problem in +% small scale and increasing it when relative +% convergence is reached. +% 'image' : Initialization using a user specified +% image. Not recommended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% +% 'Verbose' 1 or 0. Default is 1. Gives information about the +% progress of the algorithm. +% 'QualMeas' Asks the algorithm for a set of quality measurement +% parameters. Input should contain a cell array of desired +% quality measurement names. Example: {'CC','RMSE','MSSIM'} +% These will be computed in each iteration. +% 'OrderStrategy' Chooses the subset ordering strategy. Options are +% 'ordered' : uses them in the input order, but divided +% 'random' : orders them randomly +% 'angularDistance': chooses the next subset with the +% biggest angular distance with the ones used. +% 'redundancy_weighting': true or false. Default is true. Applies data +% redundancy weighting to projections in the update step +% (relevant for offset detector geometry) +% 'groundTruth' an image as ground truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the TIGRE Toolbox +% +% Copyright (c) 2015, University of Bath and +% CERN-European Organization for Nuclear Research +% All rights reserved. +% +% License: Open Source under BSD. +% See the full license at +% https://github.com/CERN/TIGRE/blob/master/LICENSE +% +% Contact: tigre.toolbox@gmail.com +% Codes: https://github.com/CERN/TIGRE/ +% Coded by: Ander Biguri +%-------------------------------------------------------------------------- + +%% Deal with input parameters + if nargin < 4 + error('TIGRE:SART:InvalidInput','Not enough inputs') + end + [res,resL2,qualMeasOut] = OS_SART(proj,geo,angles,niter, varargin{:}, 'BlockSize', 1); +end diff --git a/MATLAB/Algorithms/SIRT.m b/MATLAB/Algorithms/SIRT.m index d79fc913..145ecdb3 100644 --- a/MATLAB/Algorithms/SIRT.m +++ b/MATLAB/Algorithms/SIRT.m @@ -1,364 +1,64 @@ -function [res,resL2,qualMeasOut]=SIRT(proj,geo,angles,niter,varargin) -% SIRT solves Cone Beam CT image reconstruction using Oriented Subsets -% Simultaneous Algebraic Reconstruction Technique algorithm -% -% SIRT(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem -% using the projection data PROJ taken over ALPHA angles, corresponding -% to the geometry described in GEO, using NITER iterations. -% -% SIRT(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The -% possible options in OPT are: -% -% -% 'lambda': Sets the value of the hyperparameter. Default is 1 -% -% 'lambda_red': Reduction of lambda. Every iteration -% lambda=lambdared*lambda. Default is 0.95 -% -% 'Init': Describes different initialization techniques. -% 'none' : Initializes the image to zeros (default) -% 'FDK' : Initializes image to FDK reconstruction -% 'multigrid': Initializes image by solving the problem in -% small scale and increasing it when relative -% convergence is reached. -% 'image' : Initialization using a user specified -% image. Not recommended unless you really -% know what you are doing. -% 'InitImg' an image for the 'image' initialization. Avoid. -% -% 'Verbose' 1 or 0. Default is 1. Gives information about the -% progress of the algorithm. -% 'QualMeas' Asks the algorithm for a set of quality measurement -% parameters. Input should contain a cell array of desired -% quality measurement names. Example: {'CC','RMSE','MSSIM'} -% These will be computed in each iteration. -% 'redundancy_weighting': true or false. Default is true. Applies data -% redundancy weighting to projections in the update step -% (relevant for offset detector geometry) -% 'groundTruth' an image as ground truth, to be used if quality measures -% are requested, to plot their change w.r.t. this known -% data. -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and -% CERN-European Organization for Nuclear Research -% All rights reserved. -% -% License: Open Source under BSD. -% See the full license at -% https://github.com/CERN/TIGRE/blob/master/LICENSE -% -% Contact: tigre.toolbox@gmail.com -% Codes: https://github.com/CERN/TIGRE/ -% Coded by: Ander Biguri -%-------------------------------------------------------------------------- - -%% Deal with input parameters - -[lambda,res,lambdared,verbose,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); - -measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); -if ~any(isnan(gt(:))) - QualMeasOpts{end+1}='error_norm'; - res_prev=gt; - clear gt -end -if nargout<3 && measurequality - warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them") - measurequality=false; -end -qualMeasOut=zeros(length(QualMeasOpts),niter); - -resL2=zeros(1,niter); -if nargout>1 - computeL2=true; -else - computeL2=false; -end - -%% Create weighting matrices - -% Projection weight, W -% Projection weight, W -W=computeW(geo,angles,gpuids); - -% disp('Size of W matrix'); -% disp(size(W)); -clear geoaux - -% Back-Projection weight, V -V=computeV(geo,angles,{angles},{1:length(angles)},'gpuids',gpuids); - -if redundancy_weights - % Data redundancy weighting, W_r implemented using Wang weighting - % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc - - num_frames = size(proj,3); - W_r = redundancy_weighting(geo); - W_r = repmat(W_r,[1,1,num_frames]); - % disp('Size of redundancy weighting matrix'); - % disp(size(W_r)); - W = W.*W_r; % include redundancy weighting in W -end - -%% hyperparameter stuff -nesterov=false; -if ischar(lambda)&&strcmp(lambda,'nesterov') - nesterov=true; - lambda=(1+sqrt(1+4))/2; - gamma=0; - ynesterov=zeros(size(res),'single'); - ynesterov_prev=ynesterov; -end -%% Iterate - -% TODO : Add options for Stopping criteria -for ii=1:niter - if (ii==1 && verbose==1);tic;end - % If quality is going to be measured, then we need to save previous image - % THIS TAKES MEMORY! - if measurequality && ~strcmp(QualMeasOpts,'error_norm') - res_prev = res; % only store if necessary - end - % --------- Memory expensive----------- - % proj_err=proj-Ax(res,geo,angles); % (b-Ax) - % weighted_err=W.*proj_err; % W^-1 * (b-Ax) - % backprj=Atb(weighted_err,geo,angles); % At * W^-1 * (b-Ax) - % weight_backprj=bsxfun(@times,1./V,backprj); % V * At * W^-1 * (b-Ax) - % res=res+lambda*weight_backprj; % x= x + lambda * V * At * W^-1 * (b-Ax) - % ------------------------------------ - % --------- Memory cheap(er)----------- - - if nesterov - % The nesterov update is quite similar to the normal update, it - % just uses this update, plus part of the last one. - ynesterov_prev = ynesterov; - ynesterov=res + bsxfun(@times,1./V,Atb(W.*(proj-Ax(res,geo,angles,'gpuids',gpuids)),geo,angles,'gpuids',gpuids)); - res=(1-gamma)*ynesterov+gamma*ynesterov_prev; - else - res=res+lambda*bsxfun(@times,1./V,Atb(W.*(proj-Ax(res,geo,angles,'gpuids',gpuids)),geo,angles,'gpuids',gpuids)); % x= x + lambda * V * At * W^-1 * (b-Ax) - end - % ------------------------------------ - - if nonneg - res=max(res,0); - end - if nesterov - gamma=(1-lambda); - lambda=(1+sqrt(1+4*lambda^2))/2; - gamma=gamma/lambda; - else - lambda=lambda*lambdared; - end - % If quality is being measured - if measurequality - % HERE GOES - qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts); - end - - if computeL2 || nesterov - resL2(ii)=im3Dnorm(proj-Ax(res,geo,angles,'gpuids',gpuids),'L2','gpuids',gpuids); % Compute error norm2 of b-Ax - % If the error is not minimized. - if ii~=1 && resL2(ii)>resL2(ii-1) - if verbose - disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]); - end - return - end - - end - - - if (ii==1 && verbose==1) - expected_time=toc*niter; - disp('SIRT'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end -end - - - - - -end - - -function initres=init_multigrid(proj,geo,alpha) - -finalsize=geo.nVoxel; -% start with 64 -geo.nVoxel=[64;64;64]; -geo.dVoxel=geo.sVoxel./geo.nVoxel; -if any(finalsizefinalsize)=finalsize(geo.nVoxel>finalsize); - geo.dVoxel=geo.sVoxel./geo.nVoxel; - % Upsample! - % (hopefully computer has enough memory............) - [y, x, z]=ndgrid(linspace(1,size(initres,1),geo.nVoxel(1)),... - linspace(1,size(initres,2),geo.nVoxel(2)),... - linspace(1,size(initres,3),geo.nVoxel(3))); - initres=interp3(initres,x,y,z); - clear x y z -end - - -end - - -function [lambda,res,lambdared,verbose,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin) -opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'}; -defaults=ones(length(opts),1); -% Check inputs -nVarargs = length(argin); -if mod(nVarargs,2) - error('TIGRE:SIRT:InvalidInput','Invalid number of inputs') -end - -% check if option has been passed as input -for ii=1:2:nVarargs - ind=find(ismember(opts,lower(argin{ii}))); - if ~isempty(ind) - defaults(ind)=0; - else - error('TIGRE:SIRT:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); - end -end - -for ii=1:length(opts) - opt=opts{ii}; - default=defaults(ii); - % if one option is not default, then extract value from input - if default==0 - ind=double.empty(0,1);jj=1; - while isempty(ind) - ind=find(isequal(opt,lower(argin{jj}))); - jj=jj+1; - end - if isempty(ind) - error('TIGRE:SIRT:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); - end - val=argin{jj}; - end - - switch opt - % % % % % % % Verbose - case 'verbose' - if default - verbose=1; - else - verbose=val; - end - if ~is2014bOrNewer - warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b'); - verbose=false; - end - % % % % % % % hyperparameter, LAMBDA - case 'lambda' - if default - lambda=1; - elseif ischar(val)&&strcmpi(val,'nesterov') - lambda='nesterov'; % just for lowercase/uppercase - elseif length(val)>1 || ~isnumeric(val) - error('TIGRE:SIRT:InvalidInput','Invalid lambda') - else - lambda=val; - end - case 'lambda_red' - if default - lambdared=1; - else - if length(val)>1 || ~isnumeric(val) - error('TIGRE:SIRT:InvalidInput','Invalid lambda') - end - lambdared=val; - end - case 'init' - res=[]; - if default || strcmp(val,'none') - res=zeros(geo.nVoxel','single'); - continue - end - if strcmp(val,'FDK') - res=FDK(proj,geo,alpha); - continue - end - if strcmp(val,'multigrid') - res=init_multigrid(proj,geo,alpha); - continue - end - if strcmp(val,'image') - initwithimage=1; - continue - end - if isempty(res) - error('TIGRE:SIRT:InvalidInput','Invalid Init option') - end - % % % % % % % ERROR - case 'initimg' - if default - continue - end - if exist('initwithimage','var') - if isequal(size(val),geo.nVoxel') - res=single(val); - else - error('TIGRE:SIRT:InvalidInput','Invalid image for initialization'); - end - end - case 'qualmeas' - if default - QualMeasOpts={}; - else - if iscellstr(val) - QualMeasOpts=val; - else - error('TIGRE:SIRT:InvalidInput','Invalid quality measurement parameters'); - end - end - case 'nonneg' - if default - nonneg=true; - else - nonneg=val; - end - case 'gpuids' - if default - gpuids = GpuIds(); - else - gpuids = val; - end - case 'redundancy_weighting' - if default - redundancy_weights = true; - else - redundancy_weights = val; - end - case 'groundtruth' - if default - gt=nan; - else - gt=val; - end - otherwise - error('TIGRE:SIRT:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in SIRT()']); - end -end - -end \ No newline at end of file +function [res,resL2,qualMeasOut]=SIRT(proj,geo,angles,niter,varargin) +% SIRT solves Cone Beam CT image reconstruction using Oriented Subsets +% Simultaneous Algebraic Reconstruction Technique algorithm +% +% SIRT(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ALPHA angles, corresponding +% to the geometry described in GEO, using NITER iterations. +% +% SIRT(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'lambda': Sets the value of the hyperparameter. Default is 1 +% +% 'lambda_red': Reduction of lambda. Every iteration +% lambda=lambdared*lambda. Default is 0.95 +% +% 'Init': Describes different initialization techniques. +% 'none' : Initializes the image to zeros (default) +% 'FDK' : Initializes image to FDK reconstruction +% 'multigrid': Initializes image by solving the problem in +% small scale and increasing it when relative +% convergence is reached. +% 'image' : Initialization using a user specified +% image. Not recommended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% +% 'Verbose' 1 or 0. Default is 1. Gives information about the +% progress of the algorithm. +% 'QualMeas' Asks the algorithm for a set of quality measurement +% parameters. Input should contain a cell array of desired +% quality measurement names. Example: {'CC','RMSE','MSSIM'} +% These will be computed in each iteration. +% 'redundancy_weighting': true or false. Default is true. Applies data +% redundancy weighting to projections in the update step +% (relevant for offset detector geometry) +% 'groundTruth' an image as ground truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the TIGRE Toolbox +% +% Copyright (c) 2015, University of Bath and +% CERN-European Organization for Nuclear Research +% All rights reserved. +% +% License: Open Source under BSD. +% See the full license at +% https://github.com/CERN/TIGRE/blob/master/LICENSE +% +% Contact: tigre.toolbox@gmail.com +% Codes: https://github.com/CERN/TIGRE/ +% Coded by: Ander Biguri +%-------------------------------------------------------------------------- + +%% Deal with input parameters + + if nargin < 4 + error('TIGRE:SIRT:InvalidInput','Not enough inputs') + end + [res,resL2,qualMeasOut] = OS_SART(proj,geo,angles,niter, varargin{:}, 'BlockSize', length(angles)); +end