diff --git a/.gitignore b/.gitignore index 390d9a52..c0f350d9 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,8 @@ so.* *.jpg *.bmp *.tiff +*.png +*.mat ## Matlab ignore *.asv diff --git a/MATLAB/Algorithms/AB_GMRES.m b/MATLAB/Algorithms/AB_GMRES.m new file mode 100644 index 00000000..ee583d73 --- /dev/null +++ b/MATLAB/Algorithms/AB_GMRES.m @@ -0,0 +1,257 @@ +function [x,resL2,qualMeasOut]= AB_GMRES(proj,geo,angles,niter,varargin) + +% AB_GMRES solves the CBCT problem using AB_GMRES. +% This is a stable Krylov method for when the backprojector is not adjoint +% +% AB_GMRES(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% AB_GMRES(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% 'groundTruth' an image as grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +% 'backprojector' Descrives which backprojector to use +% +% *'FDK': This will use the FDK algorithm as a backprojector +% NOTE: not the backprojector TIGRE calls "FDK", but +% the actual algorithhm +% *'matched': (DEFAULT) uses the pseudo-matched backprojectors +% +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- + +%% + +[verbose,x,QualMeasOpts,gpuids,gt,bp]=parse_inputs(proj,geo,angles,varargin); + +if strcmpi(bp,'FDK') + backproject=@(proj,geo,angles,gpuids)FDK(proj,geo,angles,'gpuids',gpuids); +else + backproject=@(proj,geo,angles,gpuids)Atb(proj,geo,angles,'matched','gpuids',gpuids); +end + +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x0=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 +if nargout==1 + warning("Divergence caused by numerical issues not checked, due to its high computational cost for AB-GMRES") +else + warning("AB-GMRES takes very long to compute residual norms and quality measures and error norms"); +end + +qualMeasOut=zeros(length(QualMeasOpts),niter); +resL2=zeros(1,niter); + + +% Per Cristian Hansen: +% GMRES methods for tomographic reconstruction with an unmatched back projector + +w=zeros(numel(proj),niter+1,'single'); +r=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); +w(:,1) = r(:)/norm(r(:),2); +h=zeros(niter+1,niter); +for k=1:niter + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + x0 = x; % only store if necesary + end + if (k==1 && verbose);tic;end + + qk=Ax(backproject(reshape(w(:,k),geo.nDetector(1),geo.nDetector(2),length(angles)),geo,angles,gpuids),geo,angles,'Siddon','gpuids',gpuids); + e1=zeros(k+1,1); + e1(1)=1; + for ii=1:k + h(ii,k)=sum(qk(:).*w(:,ii)); + qk(:)=qk(:)-h(ii,k)*w(:,ii); + end + h(k+1,k)=norm(qk(:),2); + w(:,k+1)=qk(:)/h(k+1,k); + y=h(1:k+1,1:k)\(e1*norm(r(:),2)); + if measurequality + qualMeasOut(:,k)=Measure_Quality(x0,compute_res(x,w(:,1:k),y,geo,angles,backproject,gpuids) ,QualMeasOpts); + end + + if nargout>1 + resL2(k)=im3Dnorm(proj-Ax(compute_res(x,w(:,1:k),y,geo,angles,backproject,gpuids),geo,angles,'Siddon','gpuids',gpuids),'L2'); + if k>1 && resL2(k)>resL2(k-1) + x=compute_res(x,w(:,1:end-1),y(1),geo,angles,backproject,gpuids); + disp(['Algorithm stoped in iteration ', num2str(k),' due to loss of ortogonality.']) + return; + end + + end + + if (k==1 && verbose) + expected_time=toc*niter; + disp('AB-GMRES'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end +end +x=compute_res(x,w(:,1:end-1),y,geo,angles,backproject,gpuids); + + + +end + +% x is not explicit in this algorith, and its quite expensive to compute +function x=compute_res(x,w,y,geo,angles,backproject,gpuIds) + +for ii=1:size(w,2) + x=x+backproject(reshape(w(:,ii),geo.nDetector(1),geo.nDetector(2),length(angles)),geo,angles,gpuIds)*y(ii); +end + +end + +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids,gt,bp]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','groundtruth','backprojector'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:AB-GMRES: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:AB-GMRES: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 extranc 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:AB-GMRES:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:AB-GMRES:InvalidInput','Invalid Init option') + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:AB-GMRES:InvalidInput','Invalid image for initialization'); + end + end + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:AB-GMRES:InvalidInput','Invalid quality measurement parameters'); + end + end + case 'verbose' + if default + verbose=1; + else + verbose=val; + end + if ~is2014bOrNewer + warning('TIGRE:AB-GMRES:Verbose mode not available for older versions than MATLAB R2014b'); + verbose=false; + end + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + case 'backprojector' + if default + bp='matched'; + else + bp=val; + end + otherwise + error('TIGRE:AB-GMRES:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in CGLS()']); + end +end + + +end \ No newline at end of file diff --git a/MATLAB/Algorithms/ASD_POCS.m b/MATLAB/Algorithms/ASD_POCS.m index 66e2a6b1..0afc1141 100644 --- a/MATLAB/Algorithms/ASD_POCS.m +++ b/MATLAB/Algorithms/ASD_POCS.m @@ -18,7 +18,7 @@ % % 'init': Describes diferent initialization techniques. % • 'none' : Initializes the image to zeros (default) - +% % • 'FDK' : intializes image to FDK reconstrucition % 'TViter': Defines the amount of TV iterations performed per SART % iteration. Default is 20 @@ -45,6 +45,9 @@ % '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 grounf 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 @@ -66,17 +69,25 @@ %% parse inputs blocksize=1; -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[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); -if measurequality - qualMeasOut=zeros(length(QualMeasOpts),maxiter); -end - % does detector rotation exists? if ~isfield(geo,'rotDetector') geo.rotDetector=[0;0;0]; @@ -114,6 +125,11 @@ 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 necesary + end + f0=f; if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -154,7 +170,7 @@ geo.rotDetector=rotDetector; % Save copy of image. if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + 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'); @@ -223,9 +239,9 @@ end -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) +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'}; +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); @@ -272,9 +288,9 @@ warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= - % Its called beta in ASD-POCS + % Lambda + % ========================================================================= + % Its called beta in ASD-POCS case 'lambda' if default beta=1; @@ -284,8 +300,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -295,70 +311,70 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % 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 - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % TV hyperparameter - % ========================================================================= + % TV hyperparameter + % ========================================================================= case 'alpha' if default alpha=0.002; % 0.2 else alpha=val; end - % TV hyperparameter redution - % ========================================================================= + % TV hyperparameter redution + % ========================================================================= case 'alpha_red' if default alpha_red=0.95; else alpha_red=val; end - % Maximum update ratio - % ========================================================================= + % Maximum update ratio + % ========================================================================= case 'ratio' if default rmax=0.95; else rmax=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % ========================================================================= + % Order strategy + % ========================================================================= case 'orderstrategy' if default OrderStrategy='random'; else OrderStrategy=val; end - % Image Quality Measure - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -369,29 +385,35 @@ error('TIGRE:ASD_POCS:InvalidInput','Invalid quality measurement parameters'); end end - % Non negative - % ========================================================================= + % Non negative + % ========================================================================= case 'nonneg' if default nonneg=true; else nonneg=val; end - % GPU Ids - % ========================================================================= + % GPU Ids + % ========================================================================= case 'gpuids' if default gpuids = GpuIds(); else gpuids = val; end - % Data Redundancy weighting + % 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()']); diff --git a/MATLAB/Algorithms/AwASD_POCS.m b/MATLAB/Algorithms/AwASD_POCS.m index 88fc2d5b..923023d1 100644 --- a/MATLAB/Algorithms/AwASD_POCS.m +++ b/MATLAB/Algorithms/AwASD_POCS.m @@ -53,6 +53,9 @@ % '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 grounf 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 @@ -71,12 +74,20 @@ %% parse inputs blocksize=1; -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[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); -if measurequality - qualMeasOut=zeros(length(QualMeasOpts),maxiter); +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); @@ -120,6 +131,10 @@ 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 necesary + end f0=f; if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -141,7 +156,9 @@ 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)); - f(f<0)=0; + if nonneg + f(f<0)=0; + end end geo.offDetector=offDetector; @@ -150,7 +167,7 @@ geo.DSO=DSO; geo.rotDetector=rotDetector; if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); end % compute L2 error of actual image. Ax-b @@ -219,8 +236,8 @@ end -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=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'}; +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); @@ -267,8 +284,8 @@ warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= + % Lambda + % ========================================================================= case 'lambda' if default beta=1; @@ -278,8 +295,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -289,12 +306,12 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % Initial image + % ========================================================================= case 'init' if default || strcmp(val,'none') f0=zeros(geo.nVoxel','single'); - + else if strcmp(val,'FDK') f0=FDK(proj, geo, angles); @@ -302,49 +319,49 @@ error('TIGRE:AwASD_POCS:InvalidInput','Invalid init') end end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % TV hyperparameter - % ========================================================================= + % TV hyperparameter + % ========================================================================= case 'alpha' if default alpha=0.002; % 0.2 else alpha=val; end - % TV hyperparameter redution - % ========================================================================= + % TV hyperparameter redution + % ========================================================================= case 'alpha_red' if default alpha_red=0.95; else alpha_red=val; end - % Maximum update ratio - % ========================================================================= + % Maximum update ratio + % ========================================================================= case 'ratio' if default rmax=0.95; else rmax=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % ========================================================================= + %Parameter to control the amount of smoothing for pixels at the + %edges + % ========================================================================= case 'delta' if default delta=-0.005; @@ -357,8 +374,8 @@ else OrderStrategy=val; end - % Image Quality Measure - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -369,16 +386,16 @@ error('TIGRE:AwASD_POCS:InvalidInput','Invalid quality measurement parameters'); end end - % Non negative - % ========================================================================= + % Non negative + % ========================================================================= case 'nonneg' if default nonneg=true; else nonneg=val; end - % GPU Ids - % ========================================================================= + % GPU Ids + % ========================================================================= case 'gpuids' if default gpuids = GpuIds(); @@ -391,6 +408,12 @@ 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()']); diff --git a/MATLAB/Algorithms/AwPCSD.m b/MATLAB/Algorithms/AwPCSD.m index ede7717b..a16769e5 100644 --- a/MATLAB/Algorithms/AwPCSD.m +++ b/MATLAB/Algorithms/AwPCSD.m @@ -1,4 +1,4 @@ -function [ f,errorSART,errorTV,errorL2,qualMeasOut] = AwPCSD(proj,geo,angles,maxiter,varargin) +function [ f,resL2,qualMeasOut] = AwPCSD(proj,geo,angles,maxiter,varargin) %AwPCSD solves the reconstruction problem using adaptive-weighted %projection-controlled steepest descent method % @@ -39,6 +39,9 @@ % '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 grounf 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 @@ -56,19 +59,28 @@ % Coded by: Ander Biguri and Manasavee Lohvithee %-------------------------------------------------------------------------- %% parse inputs -[beta,beta_red,f,ng,verbose,epsilon,delta,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[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 -errorL2=[]; +resL2=zeros(1,niter); + -errorSART=[]; -errorTV=[]; % [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); % @@ -113,6 +125,10 @@ 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 necesary + end f0=f; if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -146,24 +162,20 @@ end %Non-negativity projection on all pixels - f=max(f,0); + if nonneg + f=max(f,0); + end geo.offDetector=offDetector; geo.offOrigin=offOrigin; if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + 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 errorSART - errorSARTnow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2'); - errorSART=[errorSART errorSARTnow]; - - % Compute change in the image after last SART iteration dp_vec=(f-f0); @@ -190,11 +202,6 @@ delta_p_first=im3Dnorm((Ax(f0,geo,angles,'interpolated','gpuids',gpuids))-proj,'L2'); end - %Compute errorTV - errorTVnow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2'); - errorTV=[errorTV errorTVnow]; - - % Reduce SART step beta=beta*beta_red; @@ -243,8 +250,8 @@ end -function [beta,beta_red,f0,ng,verbose,epsilon,delta,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','maxl2err','delta','qualmeas','nonneg','gpuids','redundancy_weighting'}; +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); @@ -291,8 +298,8 @@ warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= + % Lambda + % ========================================================================= case 'lambda' if default beta=1; @@ -302,8 +309,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -313,12 +320,12 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % Initial image + % ========================================================================= case 'init' if default || strcmp(val,'none') f0=zeros(geo.nVoxel','single'); - + else if strcmp(val,'FDK') f0=FDK(proj, geo, angles); @@ -326,33 +333,33 @@ error('TIGRE:AwPCSD:InvalidInput','Invalid init') end end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % ========================================================================= + %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 - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -363,16 +370,16 @@ error('TIGRE:AwPCSD:InvalidInput','Invalid quality measurement parameters'); end end - % Non negative - % ========================================================================= + % Non negative + % ========================================================================= case 'nonneg' if default nonneg=true; else nonneg=val; end - % GPU Ids - % ========================================================================= + % GPU Ids + % ========================================================================= case 'gpuids' if default gpuids = GpuIds(); @@ -385,6 +392,12 @@ 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()']); diff --git a/MATLAB/Algorithms/BA_GMRES.m b/MATLAB/Algorithms/BA_GMRES.m new file mode 100644 index 00000000..e12db5ab --- /dev/null +++ b/MATLAB/Algorithms/BA_GMRES.m @@ -0,0 +1,253 @@ +function [x,resL2,qualMeasOut]= BA_GMRES(proj,geo,angles,niter,varargin) + +% BA_GMRES solves the CBCT problem using BA_GMRES. +% This is a stable Krylov method for when the backprojector is not adjoint +% +% BA_GMRES(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% BA_GMRES(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% 'groundTruth' an image as grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +% 'backprojector' Descrives which backprojector to use +% +% *'FDK': This will use the FDK algorithm as a backprojector +% NOTE: not the backprojector TIGRE calls "FDK", but +% the actual algorithhm +% *'matched': (DEFAULT) uses the pseudo-matched backprojectors +% +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- + +%% + +[verbose,x,QualMeasOpts,gpuids,gt,bp]=parse_inputs(proj,geo,angles,varargin); + +if strcmpi(bp,'FDK') + backproject=@(proj,geo,angles,gpuids)FDK(proj,geo,angles,'gpuids',gpuids); +else + backproject=@(proj,geo,angles,gpuids)Atb(proj,geo,angles,'matched','gpuids',gpuids); +end + +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x0=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); + + +% Per Cristian Hansen: +% GMRES methods for tomographic reconstruction with an unmatched back projector + +w=zeros(prod(geo.nVoxel),niter+1,'single'); +r=backproject(proj,geo,angles,gpuids)-backproject(Ax(x,geo,angles,'Siddon','gpuids',gpuids),geo,angles,gpuids); +w(:,1) = r(:)/norm(r(:),2); + +h=zeros(niter+1,niter); +for k=1:niter + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + x0 = x; % only store if necesary + end + if (k==1 && verbose);tic;end + + qk=backproject(Ax(reshape(w(:,k),geo.nVoxel.'),geo,angles,'Siddon','gpuids',gpuids),geo,angles,gpuids); + e1=zeros(k+1,1); + e1(1)=1; + for ii=1:k + h(ii,k)=sum(qk(:).*w(:,ii)); + qk(:)=qk(:)-h(ii,k)*w(:,ii); + end + h(k+1,k)=norm(qk(:),2); + w(:,k+1)=qk(:)/h(k+1,k); + y=h(1:k+1,1:k)\(e1*norm(r(:),2)); + if measurequality + qualMeasOut(:,k)=Measure_Quality(x0,compute_res(x,w(:,1:k),y,geo) ,QualMeasOpts); + end + + if nargout>1 + aux=proj-Ax(compute_res(x,w(:,1:k),y,geo),geo,angles,'Siddon','gpuids',gpuids); + resL2(k)=im3Dnorm(aux,'L2'); + if k>1 && resL2(k)>resL2(k-1) + x=compute_res(x,w(:,1:k),y,geo); + disp(['Algorithm stoped in iteration ', num2str(k),' due to loss of ortogonality.']) + return + end + + end + + if (k==1 && verbose) + expected_time=toc*niter; + disp('BA-GMRES'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end +end +x=compute_res(x,w(:,end-1),y,geo); + + + +end + +function x=compute_res(x,w,y,geo) + +for ii=1:size(w,2) + x=x+reshape(w(:,ii),geo.nVoxel.')*y(ii); +end + +end + +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids,gt,bp]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','groundtruth','backprojector'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:BA-GMRES: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:BA-GMRES: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 extranc 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:BA-GMRES:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:BA-GMRES:InvalidInput','Invalid Init option') + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:BA-GMRES:InvalidInput','Invalid image for initialization'); + end + end + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:BA-GMRES:InvalidInput','Invalid quality measurement parameters'); + end + end + case 'verbose' + if default + verbose=1; + else + verbose=val; + end + if ~is2014bOrNewer + warning('TIGRE:BA-GMRES:Verbose mode not available for older versions than MATLAB R2014b'); + verbose=false; + end + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + case 'backprojector' + if default + bp='matched'; + else + bp=val; + end + otherwise + error('TIGRE:BA-GMRES:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in CGLS()']); + end +end + + +end \ No newline at end of file diff --git a/MATLAB/Algorithms/B_ASD_POCS_beta.m b/MATLAB/Algorithms/B_ASD_POCS_beta.m index 897533b2..5e8794a9 100644 --- a/MATLAB/Algorithms/B_ASD_POCS_beta.m +++ b/MATLAB/Algorithms/B_ASD_POCS_beta.m @@ -61,6 +61,9 @@ % '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 grounf 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 @@ -83,17 +86,24 @@ %% parse inputs blocksize=1; -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,bregman,bregman_red,bregman_iter,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,bregman,bregman_red,bregman_iter,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); angles_reorder=cell2mat(alphablocks); index_angles=cell2mat(orig_index); -if measurequality - qualMeasOut=zeros(length(QualMeasOpts),maxiter); +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); % does detector rotation exists? @@ -136,6 +146,10 @@ 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 necesary + end f0=f; if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -167,11 +181,11 @@ f=max(f,0); end end - + if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); end - + geo.offDetector=offDetector; geo.offOrigin=offOrigin; geo.DSD=DSD; @@ -249,9 +263,9 @@ end -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,bregman,bregman_red,bregman_iter,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) +function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,bregman,bregman_red,bregman_iter,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','alpha','alpha_red','ratio','maxl2err','beta','beta_red','bregman_iter','orderstrategy','nonneg','qualmeas','gpuids','redundancy_weighting'}; +opts= {'lambda','lambda_red','init','tviter','verbose','alpha','alpha_red','ratio','maxl2err','beta','beta_red','bregman_iter','orderstrategy','nonneg','qualmeas','gpuids','redundancy_weighting','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -298,9 +312,9 @@ warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= - % Its called beta in ASD-POCS + % Lambda + % ========================================================================= + % Its called beta in ASD-POCS case 'lambda' if default beta=1; @@ -310,8 +324,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -321,70 +335,70 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % 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:B_ASD_POCS_beta:InvalidInput','Invalid init') - + end end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % TV hyperparameter - % ========================================================================= + % TV hyperparameter + % ========================================================================= case 'alpha' if default alpha=0.002; % 0.2 else alpha=val; end - % TV hyperparameter redution - % ========================================================================= + % TV hyperparameter redution + % ========================================================================= case 'alpha_red' if default alpha_red=0.95; else alpha_red=val; end - % Maximum update ratio - % ========================================================================= + % Maximum update ratio + % ========================================================================= case 'ratio' if default rmax=0.95; else rmax=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % TV bregman hyperparameter - % ========================================================================= + % TV bregman hyperparameter + % ========================================================================= case 'beta' if default bregman=1; else bregman=val; end - % TV bregman hyperparameter redution - % ========================================================================= + % TV bregman hyperparameter redution + % ========================================================================= case 'beta_red' if default bregman_red=0.75; @@ -409,8 +423,8 @@ else nonneg=val; end - % Image Quality Measure - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -419,7 +433,7 @@ QualMeasOpts=val; else error('TIGRE:B_ASD_POCS_beta:InvalidInput','Invalid quality measurement parameters'); - + end end case 'gpuids' @@ -434,6 +448,12 @@ else redundancy_weights = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:B_ASD_POCS_beta:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option']); diff --git a/MATLAB/Algorithms/CGLS.m b/MATLAB/Algorithms/CGLS.m index 30162b99..db45931d 100644 --- a/MATLAB/Algorithms/CGLS.m +++ b/MATLAB/Algorithms/CGLS.m @@ -1,15 +1,15 @@ -function [x,errorL2,qualMeasOut]= CGLS(proj,geo,angles,niter,varargin) -% CGLS_CBCT solves the CBCT problem using the conjugate gradient least +function [x,resL2,qualMeasOut]= CGLS(proj,geo,angles,niter,varargin) +% CGLS solves the CBCT problem using the conjugate gradient least % squares -% -% CGLS_CBCT(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem -% using the projection data PROJ taken over ALPHA angles, corresponding +% +% CGLS(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding % to the geometry descrived in GEO, using NITER iterations. -% -% CGLS_CBCT(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% +% CGLS(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The % possible options in OPT are: -% -% +% +% % 'Init' Describes diferent initialization techniques. % * 'none' : Initializes the image to zeros (default) % * 'FDK' : intializes image to FDK reconstrucition @@ -20,107 +20,123 @@ % image. Not recomended unless you really % know what you are doing. % 'InitImg' an image for the 'image' initialization. Avoid. -% '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 grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +% 'restart' true or false. By default the algorithm will restart when +% loss of ortogonality is found. %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- % This file is part of the TIGRE Toolbox -% -% Copyright (c) 2015, University of Bath and +% +% Copyright (c) 2015, University of Bath and % CERN-European Organization for Nuclear Research % All rights reserved. % -% License: Open Source under BSD. +% 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 +% Coded by: Ander Biguri %-------------------------------------------------------------------------- %% -[verbose,x,QualMeasOpts,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[verbose,x,QualMeasOpts,gpuids,gt,restart]=parse_inputs(proj,geo,angles,varargin); +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x0=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 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)); -end +resL2=zeros(1,niter); % //doi: 10.1088/0031-9155/56/13/004 - -r=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); -if redundancy_weights - p=Atb(W_r .* r,geo,angles,'matched','gpuids',gpuids); -else +iter=0; +remember=0; +while iter1 && errorL2(ii)>errorL2(ii-1) - % OUT! - x=x-alpha*p; - if verbose - disp(['CGLS stoped in iteration N', num2str(ii),' due to divergence.']) - end - return; - end - % If step is adecuatem, then continue withg CGLS - r=r-alpha*q; - - if redundancy_weights - s=Atb(W_r .* r,geo,angles,'matched','gpuids',gpuids); - else + gamma=norm(p(:),2)^2; + for ii=iter:niter + iter=iter+1; + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + x0 = x; % only store if necesary + end + if (iter==1 && verbose);tic;end + + q=Ax(p,geo,angles,'Siddon','gpuids',gpuids); + alpha=gamma/norm(q(:),2)^2; + x=x+alpha*p; + + + + if measurequality + qualMeasOut(:,iter)=Measure_Quality(x0,x,QualMeasOpts); + end + + % The following should never happen, but the reallity is that if we use + % the residual from the algorithm, it starts diverging from this explicit residual value. + % This is an interesting fact that I believe may be caused either by + % the mismatch of the backprojection w.r.t the real adjoint, or + % numerical issues related to doing several order of magnitude + % difference operations on single precission numbers. + aux=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); + resL2(iter)=im3Dnorm(aux,'L2'); + if iter>1 && resL2(iter)>resL2(iter-1) + % we lost orthogonality, lets restart the algorithm unless the + % user asked us not to. + + % undo bad step. + x=x-alpha*p; + % if the restart didn't work. + if remember==iter || ~restart + disp(['Algorithm stoped in iteration ', num2str(iter),' due to loss of ortogonality.']) + return; + end + remember=iter; + iter=iter-1; + if verbose + disp(['Orthogonality lost, restarting at iteration ', num2str(iter) ]) + end + break + + end + % If step is adecuate, then continue withg CGLS + r=r-alpha*q; s=Atb(r,geo,angles,'matched','gpuids',gpuids); + + gamma1=norm(s(:),2)^2; + beta=gamma1/gamma; + gamma=gamma1; + p=s+beta*p; + + + if (iter==1 && verbose) + expected_time=toc*niter; + disp('CGLS'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end end - gamma1=norm(s(:),2)^2; - beta=gamma1/gamma; - gamma=gamma1; - p=s+beta*p; - - - if (ii==1 && verbose) - expected_time=toc*niter; - disp('CGLS'); - disp(['Expected duration : ',secs2hms(expected_time)]); - disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); - disp(''); - end end - end %% parse inputs' -function [verbose,x,QualMeasOpts,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) -opts= {'init','initimg','verbose','qualmeas','gpuids','redundancy_weighting'}; +function [verbose,x,QualMeasOpts,gpuids,gt,restart]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','groundtruth','restart'}; defaults=ones(length(opts),1); % Check inputs @@ -135,7 +151,7 @@ if ~isempty(ind) defaults(ind)=0; else - error('TIGRE:CGLS:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); + error('TIGRE:CGLS:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]); end end @@ -143,14 +159,14 @@ opt=opts{ii}; default=defaults(ii); % if one option isnot default, then extranc value from input - if default==0 + 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:CGLS:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + error('TIGRE:CGLS:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); end val=argin{jj}; end @@ -175,7 +191,7 @@ continue; end if isempty(x) - error('TIGRE:CGLS:InvalidInput','Invalid Init option') + error('TIGRE:CGLS:InvalidInput','Invalid Init option') end % % % % % % % ERROR case 'initimg' @@ -189,7 +205,7 @@ error('TIGRE:CGLS:InvalidInput','Invalid image for initialization'); end end - % ========================================================================= + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -200,7 +216,7 @@ error('TIGRE:CGLS:InvalidInput','Invalid quality measurement parameters'); end end - case 'verbose' + case 'verbose' if default verbose=1; else @@ -216,13 +232,19 @@ else gpuids = val; end - case 'redundancy_weighting' + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + case 'restart' if default - redundancy_weights = true; + restart=true; else - redundancy_weights = val; + restart=val; end - otherwise + otherwise error('TIGRE:CGLS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in CGLS()']); end end diff --git a/MATLAB/Algorithms/FISTA.m b/MATLAB/Algorithms/FISTA.m index 33e38e44..0cde93a5 100644 --- a/MATLAB/Algorithms/FISTA.m +++ b/MATLAB/Algorithms/FISTA.m @@ -6,10 +6,10 @@ % Processing Conference, 2018, pp.732-736. Lazy-start FISTA_mod % -% 'hyper': This parameter should approximate the largest -% eigenvalue in the A matrix in the equation Ax-b and Atb. +% 'hyper': This parameter should approximate the largest +% eigenvalue in the A matrix in the equation Ax-b and Atb. % Empirical tests show that for, the headphantom object: -% +% % geo.nVoxel = [64,64,64]' , hyper (approx=) 2.e8 % geo.nVoxel = [512,512,512]' , hyper (approx=) 2.e4 % Default: 2.e8 @@ -17,7 +17,7 @@ % • 'none' : Initializes the image to zeros (default) % • 'FDK' : intializes image to FDK reconstrucition % 'tviter': Number of iterations of Im3ddenoise to use. Default: 20 -% 'lambda': Multiplier for the tvlambda used, which is proportional to +% 'lambda': Multiplier for the tvlambda used, which is proportional to % L (hyper). Default: 0.1 % 'fista_p': P parameter for "faster" FISTA (say 1/50). Default: 1 (standard % FISTA) @@ -29,7 +29,9 @@ % parameters. Input should contain a cell array of desired % quality measurement names. Example: {'CC','RMSE','MSSIM'} % These will be computed in each iteration. - +% 'groundTruth' an image as grounf 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 @@ -46,10 +48,18 @@ % Codes: https://github.com/CERN/TIGRE/ % Coded by: Ander Biguri, Reuben Lindroos %-------------------------------------------------------------------------- -[verbose,res,tviter,hyper,lambda,p,q,QualMeasOpts,gpuids]=parse_inputs(proj,geo,angles,varargin); -%res = zeros(geo.nVoxel','single'); -measurequality=~isempty(QualMeasOpts); +[verbose,res,tviter,hyper,lambda,p,q,QualMeasOpts,gpuids,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); @@ -58,26 +68,27 @@ bm = 1/L; t = 1; r = 1/4; -p = 1; -q = 1; + for ii = 1:niter - res0 = res; + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + res_prev = res; % only store if necesary + end if (ii==1);tic;end % gradient descent step res = res + bm * 2 * Atb(proj - Ax(res,geo,angles, 'Siddon', 'gpuids', gpuids), geo, angles, 'matched', 'gpuids', gpuids); lambdaforTV = 2* bm* lambda; x_recold = x_rec; - x_rec = im3DDenoise(res,'TV',tviter,1/lambdaforTV, 'gpuids', gpuids); + x_rec = im3DDenoise(res,'TV',tviter,1/lambdaforTV, 'gpuids', gpuids); told = t; t = ( p+sqrt(q+r*t*t) ) / 2; res= x_rec + (told-1)/t * (x_rec - x_recold); - + if measurequality - qualMeasOut(:,ii)=Measure_Quality(res0,res,QualMeasOpts); + qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts); end - + if (ii==1)&&(verbose==1) expected_time=toc*niter; disp('FISTA'); @@ -90,8 +101,8 @@ end %% Parse inputs -function [verbose,f0,tviter,hyper,lambda,fista_p,fista_q,QualMeasOpts,gpuids]=parse_inputs(proj,geo,angles,argin) -opts = {'lambda','init','tviter','verbose','hyper','fista_p','fista_q','qualmeas','gpuids'}; +function [verbose,f0,tviter,hyper,lambda,fista_p,fista_q,QualMeasOpts,gpuids,gt]=parse_inputs(proj,geo,angles,argin) +opts = {'lambda','init','tviter','verbose','hyper','fista_p','fista_q','qualmeas','gpuids','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -137,8 +148,8 @@ warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Initial image - % ========================================================================= + % Initial image + % ========================================================================= case 'init' if default || strcmp(val,'none') f0=zeros(geo.nVoxel','single'); @@ -149,7 +160,7 @@ error('TIGRE:FISTA:InvalidInput','Invalid init') end end - % % % % % % % hyperparameter, LAMBDA + % % % % % % % hyperparameter, LAMBDA case 'lambda' if default lambda=0.1; @@ -158,18 +169,18 @@ else lambda=val; end - % hyperparameter - % ========================================================== + % hyperparameter + % ========================================================== case 'hyper' if default hyper = 2.e8; elseif length(val)>1 || ~isnumeric( val) error('TIGRE:FISTA:InvalidInput','Invalid lambda') - else + else hyper = val; end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default tviter = 20; @@ -178,28 +189,28 @@ else tviter = val; end - % FISTA parameter p - % ========================================================================= + % FISTA parameter p + % ========================================================================= case 'fista_p' if default - fista_p = 1; % standard FISTA, 1/50 for "faster" FISTA + fista_p = 1/50; % standard FISTA, 1/50 for "faster" FISTA elseif length(val)>1 || ~isnumeric( val) error('TIGRE:FISTA:InvalidInput','Invalid lambda') else fista_q = val; end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'fista_q' if default - fista_q = 1; % standard FISTA, 1/20 for "faster" FISTA + fista_q = 1/20; % standard FISTA, 1/20 for "faster" FISTA elseif length(val)>1 || ~isnumeric( val) error('TIGRE:FISTA:InvalidInput','Invalid lambda') else fista_q = val; end - % Image Quality Measure - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -210,14 +221,20 @@ error('TIGRE:FISTA:InvalidInput','Invalid quality measurement parameters'); end end - % GPU IDs - % ========================================================================= + % GPU IDs + % ========================================================================= case 'gpuids' if default gpuids = GpuIds(); else gpuids = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:FISTA:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in FISTA()']); end diff --git a/MATLAB/Algorithms/IRN_TV_CGLS.m b/MATLAB/Algorithms/IRN_TV_CGLS.m new file mode 100644 index 00000000..fbe75978 --- /dev/null +++ b/MATLAB/Algorithms/IRN_TV_CGLS.m @@ -0,0 +1,368 @@ +function [x_outer,resL2,qualMeasOut]= IRN_TV_CGLS(proj,geo,angles,niter,varargin) +% IRN_TV_CGLS solves the IRN_TV_CGLS problem using the conjugate gradient least +% squares with Total Variation regularization, using an inner outer scheme +% +% IRN_TV_CGLS(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% IRN_TV_CGLS(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% 'lambda' Value of regularization parameter lambda, default 10. +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% 'groundTruth' an image as grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +% 'restart' true or false. By default the algorithm will restart when +% loss of ortogonality is found. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- + +%% + +% % An Iteratively Reweighted Norm Algorithm for Total Variation Regularization +% % Paul Rodriguez and Brendt Wohlberg +% % 10.1109/ACSSC.2006.354879 + +[verbose,x0,QualMeasOpts,gpuids,lambda,niter_outer,gt,restart]=parse_inputs(proj,geo,angles,varargin); +x=x0; + +niter_break=round(niter/niter_outer); + +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x_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); +remember=0; +iter=0; +% for iii = 1:niter_outer +while iterresL2(iter-1) + % we lost orthogonality, lets restart the algorithm unless the + % user asked us not to. + % undo bad step. + x=x-alpha*p; + % if the restart didn't work. + if remember==iter || ~restart + disp(['Algorithm stoped in iteration ', num2str(iter),' due to loss of ortogonality.']) + return; + end + remember=iter; + iter=iter-1; + if verbose + disp(['Orthogonality lost, restarting at iteration ', num2str(iter) ]) + end + break + end + + % If step is adecuate, then continue withg CGLS + r_aux_1 = r_aux_1-alpha*q_aux_1; + r_aux_2=r_aux_2-alpha*q_aux_2; + + s_aux_1 = Atb(r_aux_1 ,geo,angles,'matched','gpuids',gpuids); + s_aux_2 = sqrt(lambda) * Ltx (W, r_aux_2); + s = s_aux_1 + s_aux_2; + + gamma1=norm(s(:),2)^2; + beta=gamma1/gamma; + gamma=gamma1; + p=s+beta*p; + + % replaces double explicit loops. Breaks the inner loop when its + % time to do a cold restart. + if mod(iter,niter_break)==0 + break + end + if (iter==1 && verbose) + expected_time=toc*niter; + disp('IRN_TV_CGLS'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end + end + x_outer=x; + + +end +end +% % % Non-sense now, just giving output of the right dimensions +% % function out = Lx (sizeD3_1, sizeD3_2, x) +% % % L = speye(3*sizeD3_1,sizeD3_2 ); +% % if prod(size(x)) == sizeD3_2 && 3*prod(size(x(:,:,1:end-1))) == sizeD3_1 +% % out = cat(3,x(:,:,1:end-1),x(:,:,1:end-1),x(:,:,1:end-1)); +% % else +% % error('wrong dimensions') +% % end +% % end +% % +% % function out = Ltx (sizeD3_1, sizeD3_2, x) +% % % L = speye(3*sizeD3_1,sizeD3_2 ); +% % if prod(size(x)) == sizeD3_1 +% % out = x(:,:,1:512); +% % else +% % error('wrong dimensions') +% % end +% % end + +% How to make this efficient? + +function W = build_weights ( x) +% Directional discrete derivatives +% Reflective boundary conditions +Dxx=zeros(size(x),'single'); +Dyx=zeros(size(x),'single'); +Dzx=zeros(size(x),'single'); + +Dxx(1:end-1,:,:)=x(1:end-1,:,:)-x(2:end,:,:); +Dyx(:,1:end-1,:)=x(:,1:end-1,:)-x(:,2:end,:); +Dzx(:,:,1:end-1)=x(:,:,1:end-1)-x(:,:,2:end); + +W = (Dxx.^2+Dyx.^2+Dzx.^2+1e-6).^(-1/4); % Fix this... + +end + +function out = Lx (W, x) +% Directional discrete derivatives +% Reflective boundary conditions +Dxx=zeros(size(x),'single'); +Dyx=zeros(size(x),'single'); +Dzx=zeros(size(x),'single'); + +Dxx(1:end-1,:,:)=x(1:end-1,:,:)-x(2:end,:,:); +Dyx(:,1:end-1,:)=x(:,1:end-1,:)-x(:,2:end,:); +Dzx(:,:,1:end-1)=x(:,:,1:end-1)-x(:,:,2:end); +% Build weights - is it better to find the right rotation and add +% tensors? +out = cat(4,W .* Dxx, W .* Dyx, W .* Dzx); +end + +function out = Ltx (W, x) +Wx_1 = W .* x(:,:,:,1); +Wx_2 = W .* x(:,:,:,2); +Wx_3 = W .* x(:,:,:,3); + +% Left here, but this is how to make a transpose +DxtWx_1=Wx_1; +DytWx_2=Wx_2; +DztWx_3=Wx_3; + +DxtWx_1(2:end-1,:,:)=Wx_1(2:end-1,:,:)-Wx_1(1:end-2,:,:); +DxtWx_1(end,:,:)=-Wx_1(end-1,:,:); + +DytWx_2(:,2:end-1,:)=Wx_2(:,2:end-1,:)-Wx_2(:,1:end-2,:); +DytWx_2(:,end,:)=-Wx_2(:,end-1,:); + +DztWx_3(:,:,2:end-1)=Wx_3(:,:,2:end-1)-Wx_3(:,:,1:end-2); +DztWx_3(:,:,end)=-Wx_3(:,:,end-1); + +out = DxtWx_1 + DytWx_2 + DztWx_3; +end +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids,lambda,niter_outer,gt,restart]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','lambda','niter_outer','groundtruth','restart'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:IRN_TV_CGLS: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:IRN_TV_CGLS: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 extranc 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:IRN_TV_CGLS:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:IRN_TV_CGLS:InvalidInput','Invalid Init option') + end + case 'niter_outer' + if default + niter_outer=10; + else + niter_outer=val; + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:IRN_TV_CGLS:InvalidInput','Invalid image for initialization'); + end + end + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:IRN_TV_CGLS:InvalidInput','Invalid quality measurement parameters'); + end + end + 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 + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'lambda' + if default + lambda=10; + else + lambda=val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + case 'restart' + if default + restart=true; + else + restart=val; + end + otherwise + error('TIGRE:IRN_TV_CGLS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in IRN_TV_CGLS()']); + end +end + + +end \ No newline at end of file diff --git a/MATLAB/Algorithms/LSMR.m b/MATLAB/Algorithms/LSMR.m new file mode 100644 index 00000000..d19e2720 --- /dev/null +++ b/MATLAB/Algorithms/LSMR.m @@ -0,0 +1,343 @@ +function [x,resL2,qualMeasOut]= LSMR(proj,geo,angles,niter,varargin) + +% LSMR solves the CBCT problem using LSMR. +% +% LSMR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ALPHA angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% LSMR(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% 'lambda' Value of parameter lambda, default 0. +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% 'groundTruth' an image as grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +% 'restart' true or false. By default the algorithm will restart when +% loss of ortogonality is found. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- +%% + +[verbose,x,QualMeasOpts,gpuids,lambda,gt,restart]=parse_inputs(proj,geo,angles,varargin); + +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x0=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); + +% David Chin-Lung Fong and Michael Saunders //doi.org/10.1137/10079687X + + +% Enumeration as given in the paper for 'Algorithm LSMR' +iter=0; +remember=0; +while iter= 2, construct and apply \tilde{P} to compute ||r_k|| + rhotilda = sqrt(rhod^2 + thetabar^2); + ctilda = rhod / rhotilda; + stilda = thetabar / rhotilda; + thetatildapre = thetatilda; + thetatilda = stilda * rhobar; + rhod = ctilda * rhobar; + % betatilda = ctilda * betad + stilda * betahat; % msl: in the orinal paper, but not used + betad = -stilda * betad + ctilda * betahat; + + % (10) Update \tilde{t}_k by forward substitution + tautilda = (zetapre - thetatildapre*tautilda) / rhotilda; + taud = (zeta - thetatilda*tautilda) / rhod; + + % (11) Compute ||r_k|| + d = d + betacheck^2; + gamma_var = d + (betad - taud)^2 + betadd^2; + aux = sqrt(gamma_var); % this is the residual teh algorithm follows, but we lose ortogonality, so we compute it explicitly + + % ||A^T r_k || is just |zetabar| + + + + % (6) Test for convergence. + % msl: I still need to implement this. + % msl: There are suggestions on the original paper. Let's talk about it! + + if measurequality + qualMeasOut(:,iter)=Measure_Quality(x0,x,QualMeasOpts); + end + % The following should never happen, but the reallity is that if we use + % the residual from the algorithm, it starts diverging from this explicit residual value. + % This is an interesting fact that I believe may be caused either by + % the mismatch of the backprojection w.r.t the real adjoint, or + % numerical issues related to doing several order of magnitude + % difference operations on single precission numbers. + aux=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); + resL2(iter)=im3Dnorm(aux,'L2'); + if iter>1 && resL2(iter)>resL2(iter-1) + % we lost orthogonality, lets restart the algorithm unless the + % user asked us not to. + + % undo bad step. + x=x-(zeta / (rho*rhobar)) * hbar; + % if the restart didn't work. + if remember==iter || ~restart + disp(['Algorithm stoped in iteration ', num2str(iter),' due to loss of ortogonality.']) + return; + end + remember=iter; + iter=iter-1; + if verbose + disp(['Orthogonality lost, restarting at iteration ', num2str(iter) ]) + end + break + end + + if (iter==1 && verbose) + expected_time=toc*niter; + disp('LSMR'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end + end +end +end + +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids, lambda,gt,restart]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','lambda','groundtruth','restart'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:LSMR: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:LSMR: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 extranc 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:LSMR:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:LSMR:InvalidInput','Invalid Init option') + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:LSMR:InvalidInput','Invalid image for initialization'); + end + end + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:LSMR:InvalidInput','Invalid quality measurement parameters'); + end + end + 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 + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'lambda' + if default + lambda = 0; + else + lambda = val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + case 'restart' + if default + restart=true; + else + restart=val; + end + otherwise + error('TIGRE:LSMR:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in LSMR()']); + end +end + + +end diff --git a/MATLAB/Algorithms/LSQR.m b/MATLAB/Algorithms/LSQR.m new file mode 100644 index 00000000..bb8f5687 --- /dev/null +++ b/MATLAB/Algorithms/LSQR.m @@ -0,0 +1,287 @@ +function [x,resL2,qualMeasOut]= LSQR(proj,geo,angles,niter,varargin) + +% LSQR solves the CBCT problem using LSQR. +% This is mathematically equivalent to CGLS. +% +% LSQR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% LSQR(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% 'groundTruth' an image as grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. +% 'restart' true or false. By default the algorithm will restart when +% loss of ortogonality is found. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- + +%% + +[verbose,x,QualMeasOpts,gpuids,gt,restart]=parse_inputs(proj,geo,angles,varargin); + +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x0=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); + + +% Paige and Saunders //doi.org/10.1145/355984.355989 + +% Enumeration as given in the paper for 'Algorithm LSQR' +iter=0; +remember=0; +while iter1 && resL2(iter)>resL2(iter-1) + % we lost orthogonality, lets restart the algorithm unless the + % user asked us not to. + + % undo bad step. + x=x-(phi / rho) * (v-w)/((theta / rho)); + % if the restart didn't work. + if remember==iter || ~restart + disp(['Algorithm stoped in iteration ', num2str(iter),' due to loss of ortogonality.']) + return; + end + remember=iter; + iter=iter-1; + if verbose + disp(['Orthogonality lost, restarting at iteration ', num2str(iter) ]) + end + break + end + + if (iter==1 && verbose) + expected_time=toc*niter; + disp('LSQR'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end + end + +end + +end +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids,gt,restart]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','groundtruth','restart'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:LSQR: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:LSQR: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 extranc 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:LSQR:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:LSQR:InvalidInput','Invalid Init option') + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:LSQR:InvalidInput','Invalid image for initialization'); + end + end + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:LSQR:InvalidInput','Invalid quality measurement parameters'); + end + end + case 'verbose' + if default + verbose=1; + else + verbose=val; + end + if ~is2014bOrNewer + warning('TIGRE:LSQR:Verbose mode not available for older versions than MATLAB R2014b'); + verbose=false; + end + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + case 'restart' + if default + restart=true; + else + restart=val; + end + otherwise + error('TIGRE:LSQR:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in CGLS()']); + end +end + + +end \ No newline at end of file diff --git a/MATLAB/Algorithms/MLEM.m b/MATLAB/Algorithms/MLEM.m index a14ee4a8..6a803c8e 100644 --- a/MATLAB/Algorithms/MLEM.m +++ b/MATLAB/Algorithms/MLEM.m @@ -16,7 +16,9 @@ % parameters. Input should contain a cell array of desired % quality measurement names. Example: {'CC','RMSE','MSSIM'} % These will be computed in each iteration. - +% 'groundTruth' an image as grounf 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 @@ -33,30 +35,38 @@ % Codes: https://github.com/CERN/TIGRE/ % Coded by: Ander Biguri %-------------------------------------------------------------------------- -[verbose,res,QualMeasOpts,gpuids]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); - -if measurequality - qualMeasOut=zeros(length(QualMeasOpts),niter); +[verbose,res,QualMeasOpts,gpuids,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); + res = max(res,0); -W=Atb(ones(size(proj),'single'),geo,angles,'matched','gpuids',gpuids); -W(W<=0.) = inf; +% Projection weight, W +W=computeW(geo,angles,gpuids); -tic for ii=1:niter - res0 = res; - + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + res_prev = res; % only store if necesary + end + if (ii==1);tic;end + den = Ax(res,geo,angles,'gpuids',gpuids); den(den<=0.)=inf; - auxMLEM=proj./den; - imgupdate = Atb(auxMLEM, geo,angles,'matched','gpuids',gpuids)./W; + imgupdate = Atb(proj./den, geo,angles,'matched','gpuids',gpuids)./W; res = max(res.*imgupdate,0.); if measurequality - qualMeasOut(:,ii)=Measure_Quality(res0,res,QualMeasOpts); + qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts); end if (ii==1)&&(verbose==1) @@ -71,8 +81,8 @@ end %% Parse inputs -function [verbose,f0,QualMeasOpts,gpuids]=parse_inputs(proj,geo,angles,argin) -opts = {'verbose','init','qualmeas','gpuids'}; +function [verbose,f0,QualMeasOpts,gpuids,gt]=parse_inputs(proj,geo,angles,argin) +opts = {'verbose','init','qualmeas','gpuids','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -150,6 +160,12 @@ else gpuids = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:MLEM:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in MLEM()']); end diff --git a/MATLAB/Algorithms/OS_ASD_POCS.m b/MATLAB/Algorithms/OS_ASD_POCS.m index 6c536a80..1332a933 100644 --- a/MATLAB/Algorithms/OS_ASD_POCS.m +++ b/MATLAB/Algorithms/OS_ASD_POCS.m @@ -55,6 +55,9 @@ % '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 grounf 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 @@ -73,8 +76,18 @@ %-------------------------------------------------------------------------- %% parse inputs -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,blocksize,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,blocksize,OrderStrategy,nonneg,QualMeasOpts,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); % first order the projection angles [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); @@ -108,8 +121,6 @@ W = W.*W_r; % include redundancy weighting in W end -qualMeasOut=zeros(length(QualMeasOpts),maxiter); - stop_criteria=0; @@ -120,6 +131,10 @@ 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 necesary + end f0=f; if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -164,7 +179,7 @@ % Save copy of image. fres=f; if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + 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'); @@ -229,9 +244,9 @@ end -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,block_size,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) +function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,block_size,OrderStrategy,nonneg,QualMeasOpts,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','initimg','tviter','verbose','alpha','alpha_red','ratio','maxl2err','blocksize','orderstrategy','blocksize','nonneg','qualmeas','gpuids','redundancy_weighting'}; +opts= {'lambda','lambda_red','init','initimg','tviter','verbose','alpha','alpha_red','ratio','maxl2err','blocksize','orderstrategy','blocksize','nonneg','qualmeas','gpuids','redundancy_weighting','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -278,9 +293,9 @@ warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= - % Its called beta in OS_ASD_POCS + % Lambda + % ========================================================================= + % Its called beta in OS_ASD_POCS case 'lambda' if default beta=1; @@ -290,8 +305,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -301,22 +316,22 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % 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:OS_ASD_POCS:InvalidInput','Invalid init') - + end end - % % % % % % % ERROR + % % % % % % % ERROR case 'initimg' if default continue; @@ -329,48 +344,48 @@ end end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % TV hyperparameter - % ========================================================================= + % TV hyperparameter + % ========================================================================= case 'alpha' if default alpha=0.002; % 0.2 else alpha=val; end - % TV hyperparameter redution - % ========================================================================= + % TV hyperparameter redution + % ========================================================================= case 'alpha_red' if default alpha_red=0.95; else alpha_red=val; end - % Maximum update ratio - % ========================================================================= + % Maximum update ratio + % ========================================================================= case 'ratio' if default rmax=0.95; else rmax=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % Block size for OS-SART - % ========================================================================= + % Block size for OS-SART + % ========================================================================= case 'blocksize' if default block_size=20; @@ -380,24 +395,24 @@ end block_size=val; end - % Order strategy - % ========================================================================= + % Order strategy + % ========================================================================= case 'orderstrategy' if default OrderStrategy='random'; else OrderStrategy=val; end - % Non negative - % ========================================================================= + % Non negative + % ========================================================================= case 'nonneg' if default nonneg=true; else nonneg=val; end - % Image Quality Measure - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -408,8 +423,8 @@ error('TIGRE:OS_ASD_POCS:InvalidInput','Invalid quality measurement parameters'); end end - % GPU Ids - % ========================================================================= + % GPU Ids + % ========================================================================= case 'gpuids' if default gpuids = GpuIds(); @@ -422,6 +437,12 @@ else redundancy_weights = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:OS_ASD_POCS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option']); diff --git a/MATLAB/Algorithms/OS_AwASD_POCS.m b/MATLAB/Algorithms/OS_AwASD_POCS.m index ef88c4ef..b188777e 100644 --- a/MATLAB/Algorithms/OS_AwASD_POCS.m +++ b/MATLAB/Algorithms/OS_AwASD_POCS.m @@ -59,6 +59,9 @@ % '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 grounf 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 @@ -76,8 +79,18 @@ % Coded by: Ander Biguri and Manasavee Lohvithee %% parse inputs -[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,blocksize,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[beta,beta_red,f,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,blocksize,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); % does detector rotation exists? @@ -85,7 +98,6 @@ geo.rotDetector=[0;0;0]; end -%[proj,geo] = doOffsetWang(proj,geo); %% Create weigthing matrices for the SART step % the reason we do this, instead of calling the SART fucntion is not to @@ -125,7 +137,10 @@ DSO=geo.DSO; while ~stop_criteria %POCS - f0=f; + % 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 necesary + end if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -146,7 +161,9 @@ geo.DSO=DSO(jj); end f=f+beta* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,orig_index{jj}).*(proj(:,:,orig_index{jj})-Ax(f,geo,alphablocks{:,jj},'gpuids',gpuids)),geo,alphablocks{:,jj},'gpuids',gpuids)); - f(f<0)=0; + if nonneg + f(f<0)=0; + end end geo.offDetector=offDetector; @@ -155,7 +172,7 @@ geo.DSO=DSO; geo.rotDetector=rotDetector; if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); end % compute L2 error of actual image. Ax-b @@ -201,7 +218,7 @@ %Define c_alpha as in equation 21 in the journal c=dot(dg_vec(:),dp_vec(:))/(dg*dp); %This c is examined to see if it is close to -1.0 - %disp(['Iteration = ' num2str(iter) ', c = ' num2str(c)]); + %disp(['Iteration = ' num2str(iter) ', c = ' num2str(c)]); if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>maxiter if verbose disp(['Stopping criteria met']); @@ -224,8 +241,8 @@ end -function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,block_size,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','alpha','alpha_red','ratio','maxl2err','delta','blocksize','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting'}; +function [beta,beta_red,f0,ng,verbose,alpha,alpha_red,rmax,epsilon,delta,block_size,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','blocksize','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -272,8 +289,8 @@ warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= + % Lambda + % ========================================================================= case 'lambda' if default beta=1; @@ -283,8 +300,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -294,8 +311,8 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % Initial image + % ========================================================================= case 'init' if default || strcmp(val,'none') f0=zeros(geo.nVoxel','single'); @@ -306,57 +323,57 @@ error('TIGRE:OS_AwASD_POCS:InvalidInput','Invalid init') end end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % TV hyperparameter - % ========================================================================= + % TV hyperparameter + % ========================================================================= case 'alpha' if default alpha=0.002; % 0.2 else alpha=val; end - % TV hyperparameter redution - % ========================================================================= + % TV hyperparameter redution + % ========================================================================= case 'alpha_red' if default alpha_red=0.95; else alpha_red=val; end - % Maximum update ratio - % ========================================================================= + % Maximum update ratio + % ========================================================================= case 'ratio' if default rmax=0.95; else rmax=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % ========================================================================= + %Parameter to control the amount of smoothing for pixels at the + %edges + % ========================================================================= case 'delta' if default delta=-0.005; else delta=val; end - % Block size for OS-SART - % ========================================================================= + % Block size for OS-SART + % ========================================================================= case 'blocksize' if default block_size=20; @@ -366,16 +383,16 @@ end block_size=val; end - % Order strategy - % ========================================================================= + % Order strategy + % ========================================================================= case 'orderstrategy' if default OrderStrategy='random'; else OrderStrategy=val; end - % Image Quality Measure - % ========================================================================= + % Image Quality Measure + % ========================================================================= case 'qualmeas' if default QualMeasOpts={}; @@ -404,6 +421,12 @@ else redundancy_weights = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:OS_AwASD_POCS:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in OS_AwASD_POCS()']); diff --git a/MATLAB/Algorithms/OS_AwPCSD.m b/MATLAB/Algorithms/OS_AwPCSD.m index 56190660..3da96e7d 100644 --- a/MATLAB/Algorithms/OS_AwPCSD.m +++ b/MATLAB/Algorithms/OS_AwPCSD.m @@ -1,4 +1,4 @@ -function [ f,errorSART,errorTV,errorL2,qualMeasOut] = OS_AwPCSD(proj,geo,angles,maxiter,varargin) +function [ f,resL2,qualMeasOut] = OS_AwPCSD(proj,geo,angles,maxiter,varargin) %OS_AwPCSD solves the reconstruction problem using adaptive-weighted %projection-controlled steepest descent method % @@ -51,6 +51,9 @@ % '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 grounf 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 @@ -68,20 +71,25 @@ % Coded by: Ander Biguri and Manasavee Lohvithee %-------------------------------------------------------------------------- %% parse inputs -[beta,beta_red,f,ng,verbose,epsilon,delta,blocksize,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); - -measurequality=~isempty(QualMeasOpts); +[beta,beta_red,f,ng,verbose,epsilon,delta,blocksize,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<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 -errorL2=[]; - - -errorSART=[]; -errorTV=[]; [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); % @@ -130,7 +138,10 @@ DSD=geo.DSD; DSO=geo.DSO; while ~stop_criteria %POCS - f0=f; + % 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 necesary + end if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -163,24 +174,19 @@ end %Non-negativity projection on all pixels - f=max(f,0); - + if nonneg + f=max(f,0); + end geo.offDetector=offDetector; geo.offOrigin=offOrigin; if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + 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 errorSART - errorSARTnow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2'); - errorSART=[errorSART errorSARTnow]; - - % Compute change in the image after last SART iteration dp_vec=(f-f0); @@ -207,11 +213,6 @@ delta_p_first=im3Dnorm((Ax(f0,geo,angles,'interpolated','gpuids',gpuids))-proj,'L2'); end - %Compute errorTV - errorTVnow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2'); - errorTV=[errorTV errorTVnow]; - - % Reduce SART step beta=beta*beta_red; @@ -236,15 +237,15 @@ if computeL2 geo.offOrigin=offOrigin; geo.offDetector=offDetector; - errornow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax + resL2(ii)=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 iter~=1 && resL2(ii)>resL2(ii-1) if verbose disp(['Convergence criteria met, exiting on iteration number:', num2str(iter)]); end return; end - errorL2=[errorL2 errornow]; + end @@ -260,8 +261,8 @@ end -function [beta,beta_red,f0,ng,verbose,epsilon,delta,block_size,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,argin) -opts= {'lambda','lambda_red','init','tviter','verbose','maxl2err','delta','blocksize','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting'}; +function [beta,beta_red,f0,ng,verbose,epsilon,delta,block_size,OrderStrategy,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin) +opts= {'lambda','lambda_red','init','tviter','verbose','maxl2err','delta','blocksize','orderstrategy','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -308,8 +309,8 @@ warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % Lambda - % ========================================================================= + % Lambda + % ========================================================================= case 'lambda' if default beta=1; @@ -319,8 +320,8 @@ end beta=val; end - % Lambda reduction - % ========================================================================= + % Lambda reduction + % ========================================================================= case 'lambda_red' if default beta_red=0.99; @@ -330,12 +331,12 @@ end beta_red=val; end - % Initial image - % ========================================================================= + % Initial image + % ========================================================================= case 'init' if default || strcmp(val,'none') f0=zeros(geo.nVoxel','single'); - + else if strcmp(val,'FDK') f0=FDK(proj, geo, angles); @@ -343,35 +344,35 @@ error('TIGRE:MLEM:InvalidInput','Invalid init') end end - % Number of iterations of TV - % ========================================================================= + % Number of iterations of TV + % ========================================================================= case 'tviter' if default ng=20; else ng=val; end - % Maximum L2 error to have a "good image" - % ========================================================================= + % 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 - % ========================================================================= + % ========================================================================= + % Parameter to control the amount of smoothing for pixels at the + % edges + % ========================================================================= case 'delta' if default delta=-0.005; else delta=val; end - % ========================================================================= - % Block size for OS-SART - % ========================================================================= + % ========================================================================= + % Block size for OS-SART + % ========================================================================= case 'blocksize' if default block_size=20; @@ -381,18 +382,18 @@ end block_size=val; end - % Order strategy - % ========================================================================= + % Order strategy + % ========================================================================= case 'orderstrategy' if default OrderStrategy='random'; else OrderStrategy=val; end - % ========================================================================= - % Image Quality Measure - % ========================================================================= - case 'qualmeas' + % ========================================================================= + % Image Quality Measure + % ========================================================================= + case 'qualmeas' if default QualMeasOpts={}; else @@ -402,16 +403,16 @@ error('TIGRE:OS_AwPCSD:InvalidInput','Invalid quality measurement parameters'); end end - % Non negative - % ========================================================================= + % Non negative + % ========================================================================= case 'nonneg' if default nonneg=true; else nonneg=val; end - % GPU Ids - % ========================================================================= + % GPU Ids + % ========================================================================= case 'gpuids' if default gpuids = GpuIds(); @@ -424,6 +425,12 @@ else redundancy_weights = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:OS_AwPCSD:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in OS_AwPCSD()']); diff --git a/MATLAB/Algorithms/OS_SART.m b/MATLAB/Algorithms/OS_SART.m index 64c4acf9..90c203c3 100644 --- a/MATLAB/Algorithms/OS_SART.m +++ b/MATLAB/Algorithms/OS_SART.m @@ -1,4 +1,4 @@ -function [res,errorL2,qualMeasOut]=OS_SART(proj,geo,angles,niter,varargin) +function [res,resL2,qualMeasOut]=OS_SART(proj,geo,angles,niter,varargin) % OS_SART solves Cone Beam CT image reconstruction using Oriented Subsets % Simultaneous Algebraic Reconstruction Technique algorithm % @@ -43,7 +43,9 @@ % '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 grounf truth, to be used if quality measures +% are requested, to plot their change w.r.t. this known +% data. % OUTPUTS: % % [img] will output the reconstructed image @@ -71,11 +73,23 @@ %% Deal with input parameters -[blocksize,lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[blocksize,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 @@ -120,7 +134,6 @@ ynesterov_prev=ynesterov; end %% Iterate -errorL2=[]; offOrigin=geo.offOrigin; offDetector=geo.offDetector; rotDetector=geo.rotDetector; @@ -135,8 +148,8 @@ 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 - res_prev=res; + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + res_prev = res; % only store if necesary end @@ -182,15 +195,6 @@ end - % If quality is being measured - if measurequality - - % Can save quality measure for every iteration here - % See if some image quality measure should be used for every - % iteration? - qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts); - end - % reduce hyperparameter if nesterov gamma=(1-lambda); @@ -199,23 +203,27 @@ else lambda=lambda*lambdared; end - + + if measurequality + qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts); + end + if computeL2 || nesterov % Compute error norm2 of b-Ax geo.offOrigin=offOrigin; geo.offDetector=offDetector; geo.DSD=DSD; geo.rotDetector=rotDetector; - errornow=im3Dnorm(proj-Ax(res,geo,angles,'Siddon','gpuids',gpuids),'L2'); + resL2(ii)=im3Dnorm(proj-Ax(res,geo,angles,'Siddon','gpuids',gpuids),'L2'); % If the error is not minimized - if ii~=1 && errornow>errorL2(end) % This 1.1 is for multigrid, we need to focus to only that case + if ii~=1 && resL2(ii)>resL2(ii-1) if verbose disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]); end return end % Store Error - errorL2=[errorL2 errornow]; + end % If timing was asked for if ii==1 && verbose==1 @@ -267,8 +275,8 @@ end %% Parse inputs -function [block_size,lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,alpha,argin) -opts={'blocksize','lambda','init','initimg','verbose','lambda_red','qualmeas','orderstrategy','nonneg','gpuids','redundancy_weighting'}; +function [block_size,lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin) +opts={'blocksize','lambda','init','initimg','verbose','lambda_red','qualmeas','orderstrategy','nonneg','gpuids','redundancy_weighting','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -314,7 +322,7 @@ warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b'); verbose=false; end - % % % % % % % hyperparameter, LAMBDA + % % % % % % % hyperparameter, LAMBDA case 'lambda' if default lambda=1; @@ -411,6 +419,12 @@ else redundancy_weights = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:OS_SART:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in OS_SART_CBCT()']); end diff --git a/MATLAB/Algorithms/PCSD.m b/MATLAB/Algorithms/PCSD.m index 8e6d243f..5407b787 100644 --- a/MATLAB/Algorithms/PCSD.m +++ b/MATLAB/Algorithms/PCSD.m @@ -31,6 +31,9 @@ % '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 grounf 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 @@ -49,13 +52,19 @@ %-------------------------------------------------------------------------- %% parse inputs -[beta,beta_red,f,ng,verbose,epsilon,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); +[beta,beta_red,f,ng,verbose,epsilon,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); - -if measurequality - qualMeasOut=zeros(length(QualMeasOpts),maxiter); +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') @@ -104,7 +113,10 @@ DSO=geo.DSO; %% while ~stop_criteria %POCS - f0=f; + % 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 necesary + end if (iter==0 && verbose==1);tic;end iter=iter+1; @@ -136,7 +148,9 @@ end %Non-negativity projection on all pixels - f=max(f,0); + if nonneg + f=max(f,0); + end geo.offDetector=offDetector; geo.offOrigin=offOrigin; @@ -144,7 +158,7 @@ geo.DSO=DSO; geo.rotDetector=rotDetector; if measurequality - qualMeasOut(:,iter)=Measure_Quality(f0,f,QualMeasOpts); + qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts); end % Compute L2 error of actual image. Ax-b diff --git a/MATLAB/Algorithms/README.md b/MATLAB/Algorithms/README.md new file mode 100644 index 00000000..c8605beb --- /dev/null +++ b/MATLAB/Algorithms/README.md @@ -0,0 +1,65 @@ +Algorithm Classification +=== + +There are many algorithms in TIGRE and it can be hard to chose were to start. +So let me introduce a short classification of the algorithms in TIGRE and a small suggestion on where to start. +In any case, always check the demos, as all features and algorithms of TIGRE should be showcased there. + + +There are two main ways we can classify the algorithms. First is by the way they solve the data minimization problem, +i.e. the way they make sure the image reocnstructed matches the measured data. +The second is by their regularization, i.e. by the additional constraints that we add to the problem. A common one, and the most +used one in TIGRE is Total Variation (TV), that tries to make the images "piecewise flat". + + +# Gradient descend + +Tradititonal iterative algorithms use gradient descend-like algorithm (kaczman method, ART), to solve the image. + +Non regularized: + - SIRT + - OS-SART + - SART +Not based on the maths of gradient descend, but similar (solves the proximal of the problem, using gradient descend) + - FISTA + +Regularized with TV (check demo Algorithms04) + - ASD-POCS + - B-ASD-POCS-beta + - OS-ASD-POCS + - AwASDS-POCS + - PCSD + - AwPCSD + - SART-TV (this is practically FISTA-TV) + +**Were do I start?**: For any problem, the first iterative algorithm to try is possibly OS-SART. Good convergence rate with relatively good speed per iteartion. +For TV, start with ASD-POCS, but caution, the hyperparameters have massive influence on the quality of the reconstruction. + +# Maximum likelihood + +This assumes the noise and the data follows Poisson distribution, so mostly useful for very low dose scans. + + - MLEM + +# Krylov Subspace algorithms + +These are fast coverging iterative algorithms. The have some issues with regards of semiconvergence and loss of ortogonality, so in some cases they may not produce best results, +but the main advantage is that few iterations should produce a good image. + +Non regularized: + - CGLS + - LSQR + - LSMR + - BA-GMRES + - AB-GMRES +Tikhonov regularization + - hybrid LSQR + - LSMR with non-zero lambda +TV regularization + - IRN-CGLS-TV + - hybrid-fLSQR-TV + +**Were do I start?**: LSQR. If LSQR doesn't provide stable solutions, BA-GMRES and AB-GMRES are supposed to fix that, but they require many copies (one per iteration) of the iamge, so they are very memory consuming. For regularized solutions, IRN-CGLS-TV, as the hybrid-fLSQR-TV only works for very small images, +it requires large computational and memory resources. + + \ No newline at end of file diff --git a/MATLAB/Algorithms/SART.m b/MATLAB/Algorithms/SART.m index 5d3fc71f..478cb10a 100644 --- a/MATLAB/Algorithms/SART.m +++ b/MATLAB/Algorithms/SART.m @@ -1,4 +1,4 @@ -function [res,errorL2,qualMeasOut]=SART(proj,geo,angles,niter,varargin) +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 % @@ -40,6 +40,9 @@ % '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 grounf 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 @@ -59,15 +62,26 @@ %% Deal with input parameters blocksize=1; -[lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); +[lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +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 -errorL2=[]; [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); index_angles=cell2mat(orig_index); @@ -120,8 +134,8 @@ 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 - res_prev=res; + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + res_prev = res; % only store if necesary end % reorder angles @@ -157,7 +171,7 @@ 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)); + 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); @@ -166,7 +180,6 @@ % If quality is being measured if measurequality - % HERE GOES qualMeasOut(:,ii)=Measure_Quality(res,res_prev,QualMeasOpts); end @@ -182,15 +195,14 @@ geo.offDetector=offDetector; geo.DSD=DSD; geo.rotDetector=rotDetector; - errornow=im3Dnorm(proj-Ax(res,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax + 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 && errornow>errorL2(end) + if ii~=1 && resL2(ii)>resL2(ii-1) if verbose disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]); end return end - errorL2=[errorL2 errornow]; end if (ii==1 && verbose==1) @@ -241,8 +253,8 @@ end -function [lambda,res,lambdared,verbose,QualMeasOpts,OrderStrategy,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,alpha,argin) -opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','orderstrategy','nonneg','gpuids','redundancy_weighting'}; +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); @@ -375,6 +387,12 @@ 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 diff --git a/MATLAB/Algorithms/SART_TV.m b/MATLAB/Algorithms/SART_TV.m index 50ae10bf..6b6b7bc7 100644 --- a/MATLAB/Algorithms/SART_TV.m +++ b/MATLAB/Algorithms/SART_TV.m @@ -1,4 +1,4 @@ -function [res,errorL2,qualMeasOut]=SART_TV(proj,geo,angles,niter,varargin) +function [res,resL2,qualMeasOut]=SART_TV(proj,geo,angles,niter,varargin) % SART_TV solves Cone Beam CT image reconstruction using Oriented Subsets % Simultaneous Algebraic Reconstruction Technique algorithm % @@ -47,6 +47,9 @@ % '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 grounf 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 @@ -65,16 +68,25 @@ %-------------------------------------------------------------------------- %% Deal with input parameters -[lambda,res,lamdbared,verbose,QualMeasOpts,TViter,TVlambda,OrderStrategy,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[lambda,res,lamdbared,verbose,QualMeasOpts,TViter,TVlambda,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 -errorL2=[]; blocksize=1; [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy); @@ -116,8 +128,8 @@ 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 - res_prev=res; + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + res_prev = res; % only store if necesary end @@ -150,7 +162,6 @@ % If quality is being measured if measurequality - % HERE GOES qualMeasOut(:,ii)=Measure_Quality(res,res_prev,QualMeasOpts); end @@ -164,15 +175,14 @@ geo.offDetector=offDetector; geo.DSD=DSD; geo.rotDetector=rotDetector; - errornow=im3Dnorm(proj(:,:,index_angles)-Ax(res,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax + resL2(ii)=im3Dnorm(proj(:,:,index_angles)-Ax(res,geo,angles,'gpuids',gpuids),'L2'); % Compute error norm2 of b-Ax % If the error is not minimized. - if ii~=1 && errornow>errorL2(end) + if ii~=1 && resL2(ii)>resL2(ii-1) if verbose disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]); end return end - errorL2=[errorL2 errornow]; end if (ii==1 && verbose==1) @@ -223,8 +233,8 @@ end -function [lambda,res,lamdbared,verbose,QualMeasOpts,TViter,TVlambda,OrderStrategy,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,alpha,argin) -opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','tviter','tvlambda','orderstrategy','nonneg','gpuids','redundancy_weighting'}; +function [lambda,res,lamdbared,verbose,QualMeasOpts,TViter,TVlambda,OrderStrategy,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin) +opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','tviter','tvlambda','orderstrategy','nonneg','gpuids','redundancy_weighting','groundtruth'}; defaults=ones(length(opts),1); % Check inputs nVarargs = length(argin); @@ -351,10 +361,10 @@ OrderStrategy=val; end - case 'nonneg' + case 'nonneg' if default nonneg=true; - else + else nonneg=val; end case 'gpuids' @@ -369,6 +379,12 @@ else redundancy_weights = val; end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end otherwise error('TIGRE:SART_TV:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in SART()']); end diff --git a/MATLAB/Algorithms/SIRT.m b/MATLAB/Algorithms/SIRT.m index ad62f448..53232713 100644 --- a/MATLAB/Algorithms/SIRT.m +++ b/MATLAB/Algorithms/SIRT.m @@ -1,4 +1,4 @@ -function [res,errorL2,qualMeasOut]=SIRT(proj,geo,angles,niter,varargin) +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 % @@ -35,6 +35,9 @@ % '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 grounf 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 @@ -54,10 +57,21 @@ %% Deal with input parameters -[lambda,res,lambdared,verbose,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,angles,varargin); -measurequality=~isempty(QualMeasOpts); +[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 @@ -100,15 +114,13 @@ end %% Iterate -errorL2=[]; - % 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 - res_prev=res; + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + res_prev = res; % only store if necesary end % --------- Memory expensive----------- % proj_err=proj-Ax(res,geo,angles); % (b-Ax) @@ -146,15 +158,15 @@ end if computeL2 || nesterov - errornow=im3Dnorm(proj-Ax(res,geo,angles,'gpuids',gpuids),'L2','gpuids',gpuids); % Compute error norm2 of b-Ax + 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 && errornow>errorL2(end) + if ii~=1 && resL2(ii)>resL2(ii-1) if verbose disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]); end return end - errorL2=[errorL2 errornow]; + end @@ -209,8 +221,8 @@ end -function [lambda,res,lambdared,verbose,QualMeasOpts,nonneg,gpuids,redundancy_weights]=parse_inputs(proj,geo,alpha,argin) -opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','nonneg','gpuids','redundancy_weighting'}; +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); @@ -337,6 +349,12 @@ 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 diff --git a/MATLAB/Algorithms/hybrid_LSQR.m b/MATLAB/Algorithms/hybrid_LSQR.m new file mode 100644 index 00000000..b40a3e10 --- /dev/null +++ b/MATLAB/Algorithms/hybrid_LSQR.m @@ -0,0 +1,357 @@ +function [x,resL2, qualMeasOut,lambda_vec]= hybrid_LSQR(proj,geo,angles,niter,varargin) + +% hybrid_LSQR solves the CBCT problem using LSQR. +% +% hybrid_LSQR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% hybrid_LSQR(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% 'lambda' Value of parameter lambda, default autocomputed. +% 'Noiselevel' the expected nosie level, in %, replaces lambda. +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +% 'groundTruth' an image as grounf 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- + +%% + +[verbose,x0,QualMeasOpts,gpuids, lambda, NoiseLevel,gt]=parse_inputs(proj,geo,angles,varargin); + +%%% PARAMETER CHOICE HIERARCHY: given lambda, DP, GCV + +if isnan(lambda) + if isnan(NoiseLevel) + RegParam = 'gcv'; + % Malena: if this is not good, we can use alternative formulations + else + RegParam = 'discrepit'; + % Malena: if this is slow, we can use an adaptive method + end +else + RegParam = 'given_lambda'; +end +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x_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); +% Paige and Saunders //doi.org/10.1145/355984.355989 + +% This should be one to avoid loosing orthogonality, but can be switched +% off for testing. +reorth = 1; + +% Initialise matrices +U = zeros(numel(proj), niter+1,'single'); +V = zeros(prod(geo.nVoxel), niter,'single'); % Malena: Check if prod(geo.nVoxel) is correct +B = zeros(niter+1,niter,'single'); % Projected matrix +proj_rhs = zeros(niter+1,1,'single'); % Projected right hand side +lambda_vec = zeros(niter,1); + +% Enumeration as given in the paper for 'Algorithm LSQR' +% (1) Initialize +u = proj-Ax(x0,geo,angles,'Siddon','gpuids',gpuids); +normr = norm(u(:),2); +u = u/normr; +U(:,1) = u(:); + +beta = normr; +proj_rhs(1) = normr; + +% (2) Start iterations +for ii=1:niter + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + if ii==1 + x_prev=x0; + else + x_prev = x0 + reshape(V(:,1:ii)*y,size(x0)); % only store if necesary + end + end + if (ii==1 && verbose);tic;end + + % Update V_ii + v = Atb(u,geo,angles,'matched','gpuids',gpuids); + + if ii>1 + v(:) = v(:) - beta*V(:,ii-1); + end + if reorth % Maybe change to matrix operations! + for jj = 1:ii-1 + v(:) = v(:) - (V(:,jj)'*v(:))*V(:,jj); + end + end + alpha = norm(v(:),2); % msl: do we want to check if it is 0? + v = v/alpha; + V(:,ii) = v(:); + + % Update U_{ii+1} + u = Ax(v,geo,angles,'Siddon','gpuids',gpuids) - alpha*u; + if reorth % Maybe change to matrix operations! + for jj = 1:ii-1 + u(:) = u(:) - (U(:,jj)'*u(:))*U(:,jj); + end + end + beta = norm(u(:),2); + u = u / beta; + U(:,ii+1) = u(:); + + % Update projected matrix + B(ii,ii) = alpha; + B(ii+1,ii) = beta; + % Malena. Proposed update: we should check algorithms breaks; + % 'if abs(alpha) <= eps || abs(beta) <= eps' - end and save + + % Solve the projected problem + % (using the SVD of the small projected matrix) + Bk = B(1:ii+1,1:ii); + [Uk, Sk, Vk] = svd(Bk); + if ii==1 + Sk = Sk(1,1); + else + Sk = diag(Sk); + end + rhsk = proj_rhs(1:ii+1); + rhskhat = Uk'*rhsk; + lsqr_res = abs(rhskhat(ii+1))/normr; + + + if strcmp(RegParam,'discrepit') + eta = 1; + if lsqr_res > eta*NoiseLevel + lambda = 0; + else + lambda = fzero(@(l)discrepancy(l, Bk, rhsk, eta*NoiseLevel), [0, 1e10]); + end + lambda_vec(ii) = lambda; % We should output this, maybe? + elseif strcmp(RegParam,'gcv') + lambda = fminbnd(@(l)gcv(l, rhskhat, Sk), 0, double(Sk(1))); + % Adapt from IRtools + lambda_vec(ii) = lambda; % We should output this, maybe? + elseif strcmp(RegParam,'given_lambda') + lambda_vec(ii) = lambda; + end + + Dk = Sk.^2 + lambda^2; + rhskhat = Sk .* rhskhat(1:ii); + yhat = rhskhat(1:ii)./Dk; + y = Vk * yhat; + +% resL2(ii)=norm(rhsk - Bk*y); % residual norm + x = x0 + reshape(V(:,1:ii)*y,size(x0)); + if measurequality + qualMeasOut(:,ii)=Measure_Quality(x_prev,x,QualMeasOpts); + end + aux=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); + resL2(ii)=im3Dnorm(aux,'L2'); + if ii>1 && resL2(ii)>resL2(ii-1) + disp(['Algorithm stoped in iteration ', num2str(ii),' due to loss of ortogonality.']) + return; + end + + if (ii==1 && verbose) + expected_time=toc*niter; + disp('hybrid LSQR'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end +end +x = x0 + reshape(V(:,1:ii)*y,size(x0)); + + +end + +%%% Regularization parameter choices + +function out = discrepancy(l, A, b, nnoise) +n = size(A,2); +if n == 1 + out = 0; +else + xl = [A; l*eye(n)]\[b; zeros(n,1)]; + out = (norm(A*xl -b)/norm(b))^2 - nnoise^2; +end +end + +function out = gcv(lambda, bhat, s) +% GCV for the projected problem - no weights +% If Bk is the projected matrix and Bk=Uk*Sk*Vk^T +% lambda is the regularisation parameter +% bhat is Uk'*bk +% s=diag(Sk) + +m = length(bhat); +n = length(s); + +t0 = sum(abs(bhat(n+1:m)).^2); + +s2 = abs(s) .^ 2; +lambda2 = lambda^2; + +t1 = lambda2 ./ (s2 + lambda2); +t2 = abs(bhat(1:n) .* t1) .^2; + +out = (sum(t2) + t0) / ((sum(t1)+m-n)^2); + +end + + +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids,lambda,NoiseLevel,gt]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','lambda','noiselevel','groundtruth'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:LSQR: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:LSQR: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 extranc 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:LSQR:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:LSQR:InvalidInput','Invalid Init option') + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:LSQR:InvalidInput','Invalid image for initialization'); + end + end + case 'lambda' + if default + lambda=NaN; + else + lambda=val; + end + case 'noiselevel' + if default + NoiseLevel=NaN; + else + NoiseLevel=val; + end + + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:LSQR:InvalidInput','Invalid quality measurement parameters'); + end + end + case 'verbose' + if default + verbose=1; + else + verbose=val; + end + if ~is2014bOrNewer + warning('TIGRE:LSQR:Verbose mode not available for older versions than MATLAB R2014b'); + verbose=false; + end + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + otherwise + error('TIGRE:LSQR:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in CGLS()']); + end +end + + +end \ No newline at end of file diff --git a/MATLAB/Algorithms/hybrid_fLSQR_TV.m b/MATLAB/Algorithms/hybrid_fLSQR_TV.m new file mode 100644 index 00000000..e0f1b3ed --- /dev/null +++ b/MATLAB/Algorithms/hybrid_fLSQR_TV.m @@ -0,0 +1,513 @@ +function [x,errorL2, qualMeasOut,lambda_vec]= hybrid_fLSQR_TV(proj,geo,angles,niter,varargin) + +% LSQR solves the CBCT problem using LSQR. +% This is mathematically equivalent to CGLS. +% +% LSQR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem +% using the projection data PROJ taken over ANGLES angles, corresponding +% to the geometry descrived in GEO, using NITER iterations. +% +% LSQR(PROJ,GEO,ANGLES,NITER,OPT,VAL,...) uses options and values for solving. The +% possible options in OPT are: +% +% +% 'Init' Describes diferent initialization techniques. +% * 'none' : Initializes the image to zeros (default) +% * 'FDK' : intializes image to FDK reconstrucition +% * '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 recomended unless you really +% know what you are doing. +% 'InitImg' an image for the 'image' initialization. Avoid. +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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: Malena Sabate Landman, Ander Biguri +%-------------------------------------------------------------------------- + +%% + +[verbose,x0,QualMeasOpts,gpuids, lambda, NoiseLevel,gt]=parse_inputs(proj,geo,angles,varargin); + +%%% PARAMETER CHOICE HIERARCHY: given lambda, DP, GCV + +if isnan(lambda) + if isnan(NoiseLevel) + RegParam = 'gcv'; + % Malena: if this is not good, we can use alternative formulations + else + RegParam = 'discrepit'; + % Malena: if this is slow, we can use an adaptive method + end +else + RegParam = 'given_lambda'; +end + +measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:))); +if ~any(isnan(gt(:))) + QualMeasOpts{end+1}='error_norm'; + x_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); + +% Paige and Saunders //doi.org/10.1145/355984.355989 + +% Initialise matrices +U = single(zeros(numel(proj), niter+1)); +V = single(zeros(prod(geo.nVoxel), niter)); +Z = single(zeros(prod(geo.nVoxel), niter)); % Flexible basis +M = zeros(niter+1,niter); % Projected matrix 1 +T = zeros(niter+1); % Projected matrix 2 +proj_rhs = zeros(niter+1,1); % Projected right hand side +lambda_vec = zeros(niter,1); + +z = zeros(geo.nVoxel'); + +% Enumeration as given in the paper for 'Algorithm LSQR' +% (1) Initialize +% norm_proj = norm (proj(:),2); +u = proj-Ax(x0,geo,angles,'Siddon','gpuids',gpuids); + +% Build xA0 +null_D = single((prod(geo.nVoxel)^(-1/2))*ones(geo.nVoxel')); + +k_aux = Ax(null_D,geo,angles,'Siddon','gpuids',gpuids); +xA0 = (prod(geo.nVoxel)^(-1/2))*norm(k_aux(:),2)^(-2)* (k_aux(:)'*u(:)); % scalar +xA0 = single(ones(geo.nVoxel')) * xA0; + +% Other way of computing the same, what is better? +% vecAN = Ax(null_D,geo,angles,'Siddon','gpuids',gpuids); +% [Q0, R0] = qr(vecAN(:), 0); +% x0 = Q0'*u(:); x0 = R0\x0; x0 = null_D(:)*x0; +% Ax0 = A_times_vec(A, x0) + +k_aux = (prod(geo.nVoxel)^(-1/2))*norm(k_aux(:),2)^(-2)*Atb(k_aux,geo,angles,'matched','gpuids',gpuids); + +u = u-Ax(xA0,geo,angles,'Siddon','gpuids',gpuids); + +normr = norm(u(:),2); +u = u/normr; +U(:,1) = u(:); + +if max(x0(:)) == 0 + W = ones(size(x0),'single'); +else + W = build_weights (x0); +end + +proj_rhs(1) = normr; +errorL2 = zeros(1,niter); + + +% (2) Start iterations +for ii=1:niter + if measurequality && ~strcmp(QualMeasOpts,'error_norm') + x_prev = x; % only store if necesary + end + if (ii==1 && verbose);tic;end + + % Update V, Z, and projected matrix T + v = Atb(u,geo,angles,'matched','gpuids',gpuids); + for jj = 1:ii-1 + T(jj,ii) = V(:,jj)'*v(:); + v(:) = v(:) - T(jj,ii) *V(:,jj); + end + T(ii,ii) = norm(v(:),2); + v = v/T(ii,ii); + % Choose if fixing a tolerance or a number of iterations of both + % This should maybe be done in single precisions... + z(:) = mvpE(k_aux, v(:), 'transp'); + aux_z = lsqr(@(x,tflag) Ltx(W, x, tflag), z(:), [], 50); + z(:) = lsqr(@(x,tflag) Lx(W, x, tflag), aux_z, [], 50); + z(:) = mvpE(k_aux, z(:), 'notransp'); + z = single(z); + + V(:,ii) = v(:); + Z(:,ii) = z(:); + + % Update U and projected matrix M + u = Ax(z,geo,angles,'Siddon','gpuids',gpuids); + for jj = 1:ii + M(jj,ii) = U(:,jj)'*u(:); + u(:) = u(:) - M(jj,ii) *U(:,jj); + end + M(ii+1,ii) = norm(u(:),2); + u = u / M(ii+1,ii); + U(:,ii+1) = u(:); + + % Malena. Proposed update: we should check algorithms breaks; + % 'if abs(alpha) <= eps || abs(beta) <= eps' - end and save + + %%% Solve the regularised projected problem + % (using the SVD of the small projected matrix) + Mk = M(1:ii+1,1:ii); + + % Prepare the projected regularization term + WZ = zeros(3*prod(geo.nVoxel),ii); + for jj=1:ii + % This can be done more efficiently... + % DZ can be saved and updates at each iteration + out = Lx (W, Z(:,jj), 'notransp'); + WZ(:,jj) = out; %[out{1}(:);out{2}(:);out{3}(:)]; + end + [~, ZRk] = qr(WZ(:,1:ii), 0); + ZRksq = ZRk(1:ii,1:ii); + rhsk = proj_rhs(1:ii+1); + + if strcmp(RegParam,'discrepit') + eta = 1.01; + if discrepancy_Tik(0, Mk, ZRksq, rhsk, eta*NoiseLevel) > 0 + lambda = 0; + else + lambda = fzero(@(l)discrepancy_Tik(l, Mk, ZRksq, rhsk, eta*NoiseLevel), [0, 1e10]); + end + lambda_vec(ii) = lambda; % We should output this, maybe? + elseif strcmp(RegParam,'gcv') + [Uk, ~, ~, Ck, Sk] = gsvd(Mk, ZRksq); + rhskhat = Uk'*rhsk; + if ii==1 + gammak = Ck(1)/Sk(1); + else + gammak = sqrt(diag(Ck'*Ck)./diag(Sk'*Sk)); + end + lambda = fminbnd(@(l)gcv(l, rhskhat, gammak), 0, double(gammak(ii))); + lambda_vec(ii) = lambda; % We should output this, maybe? + + elseif strcmp(RegParam,'given_lambda') + lambda_vec(ii) = lambda; + end + + MZk = [Mk; lambda*ZRksq]; + rhsZk = [rhsk; zeros(ii,1)]; + y = MZk\rhsZk; + +% errorL2(ii)=norm(rhsk - Mk*y); + + d = Z(:,1:ii)*y; + x = x0 + reshape(d,size(x0)) + xA0; + + W = build_weights (x); + + % Test for convergence. + % msl: I still need to implement this. + % msl: There are suggestions on the original paper. Let's talk about it! + + if measurequality + qualMeasOut(:,ii)=Measure_Quality(x_prev,x,QualMeasOpts); + end + aux=proj-Ax(x,geo,angles,'Siddon','gpuids',gpuids); + errorL2(ii)=im3Dnorm(aux,'L2'); + if ii>1 && resL2(ii)>resL2(ii-1) + disp(['Algorithm stoped in iteration ', num2str(ii),' due to loss of ortogonality.']) + return; + end + if (ii==1 && verbose) + expected_time=toc*niter; + disp('hybrid fLSQR TV'); + disp(['Expected duration : ',secs2hms(expected_time)]); + disp(['Expected finish time: ',datestr(datetime('now')+seconds(expected_time))]); + disp(''); + end +end + +end + +%%% Regularization parameter choices + +function out = discrepancy_Tik(lambda, A, ZRksq, b, nnoise) +n = size(A,2); +if n == 1 + out = 0; +else + xl = [A; lambda*ZRksq]\[b; zeros(size(ZRksq,1),1)]; + out = (norm(A*xl - b)/norm(b))^2 - nnoise^2; +end +end + +function out = gcv(lambda, bhat, s) +% GCV for the projected problem - no weights +% If Bk is the projected matrix and Bk=Uk*Sk*Vk^T +% lambda is the regularisation parameter +% bhat is Uk'*bk +% s=diag(Sk) + +m = length(bhat); +n = length(s); + +t0 = sum(abs(bhat(n+1:m)).^2); + +s2 = abs(s) .^ 2; +lambda2 = lambda^2; + +t1 = lambda2 ./ (s2 + lambda2); +t2 = abs(bhat(1:n) .* t1) .^2; + +out = (sum(t2) + t0) / ((sum(t1)+m-n)^2); + +end + +%%%%%% + +function W = build_weights ( x) + % Directional discrete derivatives + % Reflective boundary conditions + Dxx=zeros(size(x),'single'); + Dyx=zeros(size(x),'single'); + Dzx=zeros(size(x),'single'); + + Dxx(1:end-1,:,:)=x(1:end-1,:,:)-x(2:end,:,:); + Dyx(:,1:end-1,:)=x(:,1:end-1,:)-x(:,2:end,:); + Dzx(:,:,1:end-1)=x(:,:,1:end-1)-x(:,:,2:end); + + W = (Dxx.^2+Dyx.^2+Dzx.^2+1e-6).^(-1/4); % Fix this... + +end + +function out = Lx (W, x, transp_flag) +if strcmp(transp_flag,'transp') + x = reshape(x,[size(W),3]); + Wx_1 = W .* x(:,:,:,1); + Wx_2 = W .* x(:,:,:,2); + Wx_3 = W .* x(:,:,:,3); + + % Left here, but this is how to make a transpose + DxtWx_1=Wx_1; + DytWx_2=Wx_2; + DztWx_3=Wx_3; + + DxtWx_1(2:end-1,:,:)=Wx_1(2:end-1,:,:)-Wx_1(1:end-2,:,:); + DxtWx_1(end,:,:)=-Wx_1(end-1,:,:); + + DytWx_2(:,2:end-1,:)=Wx_2(:,2:end-1,:)-Wx_2(:,1:end-2,:); + DytWx_2(:,end,:)=-Wx_2(:,end-1,:); + + DztWx_3(:,:,2:end-1)=Wx_3(:,:,2:end-1)-Wx_3(:,:,1:end-2); + DztWx_3(:,:,end)=-Wx_3(:,:,end-1); + + out = DxtWx_1 + DytWx_2 + DztWx_3; + out = out(:); + +elseif strcmp(transp_flag,'notransp') + + x = reshape(x,size(W)); + + % Directional discrete derivatives + % Reflective boundary conditions + Dxx=zeros(size(x),'single'); + Dyx=zeros(size(x),'single'); + Dzx=zeros(size(x),'single'); + + Dxx(1:end-1,:,:)=x(1:end-1,:,:)-x(2:end,:,:); + Dyx(:,1:end-1,:)=x(:,1:end-1,:)-x(:,2:end,:); + Dzx(:,:,1:end-1)=x(:,:,1:end-1)-x(:,:,2:end); + % Build weights - is it better to find the right rotation and add + % tensors? + out = cat(4,W .* Dxx, W .* Dyx, W .* Dzx); + out = out(:); +end +end + +function out = Ltx (W, x, transp_flag) +if strcmp(transp_flag,'notransp') + x = reshape(x,[size(W),3]); + Wx_1 = W .* x(:,:,:,1); + Wx_2 = W .* x(:,:,:,2); + Wx_3 = W .* x(:,:,:,3); + + % Left here, but this is how to make a transpose + DxtWx_1=Wx_1; + DytWx_2=Wx_2; + DztWx_3=Wx_3; + + DxtWx_1(2:end-1,:,:)=Wx_1(2:end-1,:,:)-Wx_1(1:end-2,:,:); + DxtWx_1(end,:,:)=-Wx_1(end-1,:,:); + + DytWx_2(:,2:end-1,:)=Wx_2(:,2:end-1,:)-Wx_2(:,1:end-2,:); + DytWx_2(:,end,:)=-Wx_2(:,end-1,:); + + DztWx_3(:,:,2:end-1)=Wx_3(:,:,2:end-1)-Wx_3(:,:,1:end-2); + DztWx_3(:,:,end)=-Wx_3(:,:,end-1); + + out = DxtWx_1 + DytWx_2 + DztWx_3; + out = out(:); + +elseif strcmp(transp_flag,'transp') + + x = reshape(x,size(W)); + + % Directional discrete derivatives + % Reflective boundary conditions + Dxx=zeros(size(x),'single'); + Dyx=zeros(size(x),'single'); + Dzx=zeros(size(x),'single'); + + Dxx(1:end-1,:,:)=x(1:end-1,:,:)-x(2:end,:,:); + Dyx(:,1:end-1,:)=x(:,1:end-1,:)-x(:,2:end,:); + Dzx(:,:,1:end-1)=x(:,:,1:end-1)-x(:,:,2:end); + % Build weights - is it better to find the right rotation and add + % tensors? + out = cat(4,W .* Dxx, W .* Dyx, W .* Dzx); + out = out(:); +end +end + +function out = mvpE(k_aux, x , transp_flag) +if strcmp(transp_flag,'transp') + out = x(:) - k_aux(:)*sum(x(:)); +elseif strcmp(transp_flag,'notransp') + out = x(:) - (k_aux(:)'*x(:)); +end +end + +% function out = mvpEt(k_aux, x , transp_flag) +% if strcmp(transp_flag,'notransp') +% out = x(:) - k_aux(:)*(ones(size(x(:)))'*x(:)); +% elseif strcmp(transp_flag,'transp') +% out = x(:) - ones(size(x(:)))*(k_aux(:)'*x(:)); +% end +% end + + +%% parse inputs' +function [verbose,x,QualMeasOpts,gpuids, lambda, NoiseLevel,gt]=parse_inputs(proj,geo,angles,argin) +opts= {'init','initimg','verbose','qualmeas','gpuids','lambda','noiselevel','groundtruth'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:LSQR: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:LSQR: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 extranc 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:LSQR:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'init' + x=[]; + if default || strcmp(val,'none') + x=zeros(geo.nVoxel','single'); + continue; + end + if strcmp(val,'FDK') + x=FDK(proj,geo,angles); + continue; + end + if strcmp(val,'multigrid') + x=init_multigrid(proj,geo,angles); + continue; + end + if strcmp(val,'image') + initwithimage=1; + continue; + end + if isempty(x) + error('TIGRE:LSQR:InvalidInput','Invalid Init option') + end + % % % % % % % ERROR + case 'initimg' + if default + continue; + end + if exist('initwithimage','var') + if isequal(size(val),geo.nVoxel') + x=single(val); + else + error('TIGRE:LSQR:InvalidInput','Invalid image for initialization'); + end + end + case 'lambda' + if default + lambda=NaN; + else + lambda=val; + end + case 'noiselevel' + if default + NoiseLevel=NaN; + else + NoiseLevel=val; + end + + % ========================================================================= + case 'qualmeas' + if default + QualMeasOpts={}; + else + if iscellstr(val) + QualMeasOpts=val; + else + error('TIGRE:LSQR:InvalidInput','Invalid quality measurement parameters'); + end + end + case 'verbose' + if default + verbose=1; + else + verbose=val; + end + if ~is2014bOrNewer + warning('TIGRE:LSQR:Verbose mode not available for older versions than MATLAB R2014b'); + verbose=false; + end + case 'gpuids' + if default + gpuids = GpuIds(); + else + gpuids = val; + end + case 'groundtruth' + if default + gt=nan; + else + gt=val; + end + otherwise + error('TIGRE:LSQR:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in CGLS()']); + end +end + + +end \ No newline at end of file diff --git a/MATLAB/Demos/d08_Algorithms03.m b/MATLAB/Demos/d08_Algorithms03.m index 027a5cb1..8997187c 100644 --- a/MATLAB/Demos/d08_Algorithms03.m +++ b/MATLAB/Demos/d08_Algorithms03.m @@ -1,81 +1,173 @@ -%% DEMO 8: Algorithms 03. Krylov subspace -% -% -% In this demo the usage of the the Krylov subspace family is explained. -% This family of algorithms iterates trhough the eigenvectors of the -% residual (Ax-b) of the problem in descending order, achieving increased -% convergence rates comparing to SART family. -% -% In cases where the data is good quality, SART type families tend to reach -% to a better image, but when the data gets very big, or has bad quality, -% CGLS is a good and fast algorithm -%-------------------------------------------------------------------------- -%-------------------------------------------------------------------------- -% 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 -%-------------------------------------------------------------------------- -%% Initialize -clear; -close all; - -%% Define Geometry -geo=defaultGeometry('nVoxel',[512,512,512]','nDetector',[512,512]); - -%% Load data and generate projections -% see previous demo for explanation -angles=linspace(0,2*pi,100); -head=headPhantom(geo.nVoxel); -projections=Ax(head,geo,angles,'interpolated'); -noise_projections=addCTnoise(projections); - -%% Usage CGLS -% -% -% CGLS has the common 4 inputs for iterative algorithms in TIGRE: -% -% Projections, geometry, angles, and number of iterations -% -% Additionally it contains optional initialization tehcniques, but we -% reccomend not using them. CGLS is already quite fast and using them may -% lead to divergence. -% The options are: -% 'Init' Describes diferent initialization techniques. -% � 'none' : Initializes the image to zeros (default) -% � 'FDK' : intializes image to FDK reconstrucition -% � '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 recomended unless you really -% know what you are doing. -% 'InitImg' an image for the 'image' initialization. Avoid. - -% use CGLS -[imgCGLS, errL2CGLS]=CGLS(noise_projections,geo,angles,60); -% SIRT for comparison. -[imgSIRT,errL2SIRT]=SIRT(noise_projections,geo,angles,60); - -%% plot results -% -% We can see that CGLS gets to the same L2 error in less amount of -% iterations. - -plot([errL2SIRT;[errL2CGLS nan(1,length(errL2SIRT)-length(errL2CGLS))]]'); -title('L2 error') -legend('SIRT','CGLS') - -% plot images -plotImg([imgCGLS imgSIRT],'Dim','Z','Step',2) -%plot errors -plotImg(abs([head-imgCGLS head-imgSIRT]),'Dim','Z','Slice',64) +%% DEMO 8: Algorithms 03. Krylov subspace +% +% +% In this demo the usage of the the Krylov subspace family is explained. +% This family of algorithms iterates trhough the eigenvectors of the +% residual (Ax-b) of the problem in descending order, achieving increased +% convergence rates comparing to SART family. +% +% In cases where the data is good quality, SART type families tend to reach +% to a better image, but when the data gets very big, or has bad quality, +% CGLS is a good and fast algorithm +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% 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 +%-------------------------------------------------------------------------- +%% Initialize +clear; +close all; + +%% Define Geometry +geo=defaultGeometry('nVoxel',[128,128,128]','nDetector',[128,128]); + +%% Load data and generate projections +% see previous demo for explanation +angles=linspace(0,2*pi,100); +head=headPhantom(geo.nVoxel); +projections=Ax(head,geo,angles,'interpolated'); +projections=addCTnoise(projections); + +%% Algorithhms +% these algorithms are showcased and referenced in the paper: +% "On Krylov Methods for Large Scale CBCT Reconstruction", M Sabate Landman +% et al. To be published. + +niter=30; + +% We also introduce a new feature in this demo too: the posibility to add the +% ground Truth image for error norms, if you have it. This is only for +% algorithm study, as in real situations you don't have this (otherwise no +% need for recon, right?) but its good to show performance and +% semiconvergence. + +% SIRT, for comparison +[imgsirt,ressirt,errorsirt]=SIRT(projections,geo,angles,niter,'groundTruth',head); +% CGLS, traditional and mostly used Krylov algorithm +[imgcgls,rescgls,errcgls]=CGLS(projections,geo,angles,niter,'groundTruth',head); +% LSQR, stable version of CGLS +[imglsqr,reslsqr,errlsqr]=LSQR(projections,geo,angles,niter,'groundTruth',head); +% LSMR, a different krylov algorithm. It has Tikhonov regularization for +% stability, controled with the parameter 'lambda' +[imglsmr,reslsmr,errlsmr]=LSMR(projections,geo,angles,niter,'groundTruth',head,'lambda',0); +[imglsmr2,reslsmr2,errlsmr2]=LSMR(projections,geo,angles,niter,'groundTruth',head,'lambda',30); +% hybrid LSQR, a reorthogonalizing version of LSQR that can lead to better +% stability. +[imghlsqr,reshlsqr,errhlsqr]=hybrid_LSQR(projections,geo,angles,niter,'groundTruth',head); +% ABBA-GMRES are a set of algorithms that accept unmatched backprojectors, +% which TIGRE (and most CT libraries) have. +% read more at: "GMRES methods for tomographic reconstruction with an unmatched back projector", +% Per Christian Hansen, Ken Hayami, Keiichi Morikuni +[imgabgmres,resabgmres,errabgmres]=AB_GMRES(projections,geo,angles,niter,'groundTruth',head); +[imgbagmres,resbagmres,errbagmres]=BA_GMRES(projections,geo,angles,niter,'groundTruth',head); +% these can have other backprojectors, not just the standard ones. You can +% backproject with FDK for ultra-fast convergence, for example. +[imgabgmresfdk,resabgmresfdk,errabgmresfdk]=AB_GMRES(projections,geo,angles,niter,'groundTruth',head,'backprojector','FDK'); +[imgbagmresfdk,resbagmresfdk,errbagmresfdk]=BA_GMRES(projections,geo,angles,niter,'groundTruth',head,'backprojector','FDK'); + +%% plot results +% +% Notice semiconvergence (this also happens with the other algorithms in +% other demos). i.e. residual goes down, but error goes up. +% Notice how LSMR with regularization doesn't suffer from this. +% Notice fast covnergence of ABBA with FDK backprojector, but not as good +% as the other algorithms after long iterations +% Notice that they are sometimes stopped early, due to divergence +% (unmatched backprojector+numerical issues) + +% clean up data: +errorsirt(errorsirt==0)=nan; +errcgls(errcgls==0)=nan; +errlsqr(errlsqr==0)=nan; +errlsmr(errlsmr==0)=nan; +errlsmr2(errlsmr2==0)=nan; +errhlsqr(errhlsqr==0)=nan; +errabgmres(errabgmres==0)=nan; +errbagmres(errbagmres==0)=nan; +errabgmresfdk(errabgmresfdk==0)=nan; +errbagmresfdk(errbagmresfdk==0)=nan; +%% plots N1: +set(0,'defaultTextInterpreter','latex'); +set(0,'DefaultTextFontName','Helvetica') +figure(1) +subplot(121) +semilogy(errorsirt/norm(head(:)),'linewidth',2); +hold on; +semilogy(errcgls/norm(head(:)),'linewidth',2); +semilogy(errlsqr/norm(head(:)),'linewidth',2); +semilogy(errlsmr/norm(head(:)),'linewidth',2); +semilogy(errlsmr2/norm(head(:)),'linewidth',2); +semilogy(errhlsqr/norm(head(:)),'linewidth',2); +semilogy(errabgmres/norm(head(:)),'linewidth',2); +semilogy(errbagmres/norm(head(:)),'linewidth',2); +semilogy(errabgmresfdk/norm(head(:)),'linewidth',2); +semilogy(errbagmresfdk/norm(head(:)),'linewidth',2); +xlim([1,niter]); +ylim([7*10^-2,10^0]) +xlabel("Iteration number") +ylabel("$\|x-x_{gt}\|/\|x_{gt}\|$",'interpreter','latex') +title("Error norm") +legend(["SIRT","CGLS","LSQR","LSMR $\lambda=0$","LSMR $\lambda=30$","hybrid LSQR","AB-GMRES","BA-GMRES","AB-GRMES ($B_{FDK}$)","BA-GRMES ($B_{FDK}$)"],'NumColumns',2,'interpreter','latex') + +% set(gca,'fontsize',16) +set(gcf, 'Color', 'w'); +grid on + +subplot(122) + +semilogy(ressirt/norm(projections(:)),'linewidth',2); +hold on; +semilogy(rescgls/norm(projections(:)),'linewidth',2); +semilogy(reslsqr/norm(projections(:)),'linewidth',2); +semilogy(reslsmr/norm(projections(:)),'linewidth',2); +semilogy(reslsmr2/norm(projections(:)),'linewidth',2); +semilogy(reshlsqr/norm(projections(:)),'linewidth',2); + +semilogy(resabgmres/norm(projections(:)),'linewidth',2); +semilogy(resbagmres/norm(projections(:)),'linewidth',2); +semilogy(resabgmresfdk/norm(projections(:)),'linewidth',2); +semilogy(resbagmresfdk/norm(projections(:)),'linewidth',2); +xlim([1,niter]); +ylim([9*10^-3,10^0]) +xlabel("Iteration number") +ylabel("$\|Ax-b\|/\|b\|$",'interpreter','latex') +title("Residual norm") +legend(["SIRT","CGLS","LSQR","LSMR $\lambda=0$","LSMR $\lambda=30$","hybrid LSQR","AB-GMRES","BA-GMRES","AB-GRMES ($B_{FDK}$)","BA-GRMES ($B_{FDK}$)"],'NumColumns',2,'interpreter','latex') +grid on +set(gcf, 'Color', 'w'); +set(gcf, 'Units', 'Inches', 'Position', [1, 1, 15/1.2, 7/1.2], 'PaperUnits', 'Inches', 'PaperSize', [9/1.2 7/1.2]) +% set(gca,'fontsize',16) +set(0,'defaultTextInterpreter','tex'); + +%% Hybrid methods with different regularisation parameter choices + +% you can explicitly defined the parameter in the mathematical terms +[imghLSQR_l10, residual_hLSQR_l10]=hybrid_LSQR(projections,geo,angles,30, 'lambda', 10); +% You can give it a "noise level" (in %) instead, and it will chose the lamba inside +[imghLSQR_nl002, residual_hLSQR_nl002,~, lambda_vec_nl002]=hybrid_LSQR(projections,geo,angles,30, 'NoiseLevel', 0.02); +% if you don't give it any, it will use Generalized Cross Validation to approximate a good lambda +[imghLSQR_gcv, residual_hLSQR_gcv,~, lambda_vec_gcv]=hybrid_LSQR(projections,geo,angles,30); + +%% plot images +plotImg([imgcgls imghLSQR_l10 imghLSQR_nl002 imghLSQR_gcv],'Dim','Z','Step',2) + +% Look at the parameters +figure +plot(lambda_vec_nl002) +hold on +plot(lambda_vec_gcv) +title("lambda vs iteratios") +legend(["Noise Level", "Generalized Cross Validation"]) + +%% diff --git a/MATLAB/Demos/d09_Algorithms04.m b/MATLAB/Demos/d09_Algorithms04.m index a02e0da9..838c92a8 100644 --- a/MATLAB/Demos/d09_Algorithms04.m +++ b/MATLAB/Demos/d09_Algorithms04.m @@ -191,23 +191,52 @@ imgSARTTV=SART_TV(noise_projections,geo,angles,50,'TViter',100,'TVlambda',50); + +% IRN_TV_CGLS +%========================================================================== +%========================================================================== +% CGLS with TV regularization in an inner/outer iteration scheme +% +% 'lambda' hyperparameter in TV norm. It gives the ratio of +% importance of the image vs the minimum total variation. +% default is 15. Lower means less TV denoising. +% +% 'niter_outer' Number of outer iterations. Each outer iteration will +% perform niter number of inner iterations, in the example +% below, 20.Albeit this seems that it does many more +% iterations than the other algorithms, this is an inherently +% faster algorithm, both in convergence and time. + +imgIRN_TV_CGLS=IRN_TV_CGLS(noise_projections,geo,angles,20,'lambda',5,'niter_outer',2); +% hybrid_fLSQR_TV +%========================================================================== +%========================================================================== +% flexyble hybrod LSQR. Has TV regularization with reorthogonalization. +% +% 'lambda' hyperparameter in TV norm. It gives the ratio of +% importance of the image vs the minimum total variation. +% default is 15. Lower means less TV denoising. +% +imgflsqr=hybrid_fLSQR_TV(projections,geo,angles,20,'lambda',5); + %% Lets visualize the results % Notice the smoother images due to TV regularization. % -% thorax OS-SART ASD-POCS +% head OS-SART ASD-POCS % -% OSC-TV B-ASD-POCS-beta SART-TV +% OS-ASD-POCS B-ASD-POCS-beta SART-TV +% +% IRN-TV-CGLS hybrid-fLSQR heads -plotImg([ imgOSASDPOCS imgBASDPOCSbeta imgSARTTV; head imgOSSART imgASDPOCS ] ,'Dim','Z','Step',2,'clims',[0 1]) +plotImg([imgIRN_TV_CGLS imgflsqr head; + imgOSASDPOCS imgBASDPOCSbeta imgSARTTV; + head, imgOSSART, imgASDPOCS ] ,'Dim','Z','Step',2,'clims',[0 1]) % error - -plotImg(abs([ head-imgOSASDPOCS head-imgBASDPOCSbeta head-imgSARTTV;head-head head-imgOSSART head-imgASDPOCS ]) ,'Dim','Z','Slice',64) - - - +plotImg(abs([head-imgIRN_TV_CGLS head-imgflsqr head-head; + head-imgOSASDPOCS head-imgBASDPOCSbeta head-imgSARTTV; + head-head, head-imgOSSART, head-imgASDPOCS ]) ,'Dim','Z','Slice',64,'clims',[0 0.1]) %% - % Obligatory XKCD reference: https://xkcd.com/833/ diff --git a/MATLAB/Utilities/Quality_measures/Measure_Quality.m b/MATLAB/Utilities/Quality_measures/Measure_Quality.m index 53044d13..33ce3233 100644 --- a/MATLAB/Utilities/Quality_measures/Measure_Quality.m +++ b/MATLAB/Utilities/Quality_measures/Measure_Quality.m @@ -21,21 +21,19 @@ opt=QualMeasOpts{ii}; switch opt - %%%%%%%%RMSE case 'RMSE' q=RMSE(res_prev,res); - %%%%%%CC case 'CC' q=CC(res_prev,res); - %%%%%%MSSIM case 'MSSIM' q=MSSIM(res_prev,res); case 'UQI' q=UQI(res_prev,res); - + case 'error_norm' + q=im3Dnorm(res_prev-res,'L2'); end diff --git a/Python/demos/d08_Algorithms03.py b/Python/demos/d08_Algorithms03.py index 9eb37078..6222585f 100644 --- a/Python/demos/d08_Algorithms03.py +++ b/Python/demos/d08_Algorithms03.py @@ -68,10 +68,26 @@ # know what you are doing. # 'InitImg' an image for the 'image' initialization. Avoid. -# use CGLS -imgCGLS, errL2CGLS = algs.cgls(noise_projections, geo, angles, 60, computel2=True) +# # use CGLS +imgCGLS, normL2CGLS = algs.cgls(noise_projections, geo, angles, 30, computel2=True) +# use LSQR +imgLSQR, normL2LSQR = algs.lsqr(noise_projections, geo, angles, 30, computel2=True) +# use LSMR +imgLSMR, normL2LSMR = algs.lsmr(noise_projections, geo, angles, 30, computel2=True,lmbda=0) +imgLSMR2, normL2LSMR2 = algs.lsmr(noise_projections, geo, angles, 30, computel2=True,lmbda=30) +# use LSQR +imghLSQR, normhL2LSQR = algs.hybrid_lsqr(noise_projections, geo, angles, 30, computel2=True) + +# AB/BA-GMRES +imgabgmres, normhabgmres = algs.ab_gmres(noise_projections, geo, angles, 30, computel2=True) +imgbagmres, normhbagmres = algs.ba_gmres(noise_projections, geo, angles, 30, computel2=True) +# # AB/BA-GMRES with FDK backprojector +imgabgmresfdk, normhabgmresfdk = algs.ab_gmres(noise_projections, geo, angles, 30, computel2=True,backprojector="FDK") +imgbagmresfdk, normhbagmresfdk = algs.ba_gmres(noise_projections, geo, angles, 30, computel2=True,backprojector="FDK") + + # SIRT for comparison. -imgSIRT, errL2SIRT = algs.sirt(noise_projections, geo, angles, 60, computel2=True) +imgSIRT, normL2SIRT = algs.sirt(noise_projections, geo, angles, 60, computel2=True) #%% plot results # @@ -79,13 +95,14 @@ # iterations. -plt.plot(np.vstack((errL2CGLS[0, :], errL2SIRT[0, :])).T) + +plt.plot(np.vstack((normL2CGLS[0, :], normL2SIRT[0, 0:30],normL2LSMR[0, :],normL2LSMR2[0, :],normhL2LSQR[0, :],normhabgmres[0,:],normhbagmres[0,:],normhabgmresfdk[0,:],normhbagmresfdk[0,:])).T) plt.title("L2 error") plt.xlabel("Iteration") -plt.ylabel("$ log_{10}(|Ax-b|) $") -plt.gca().legend(("CGLS", "SIRT")) +plt.ylabel("$ |Ax-b| $") +plt.gca().legend(("CGLS", "SIRT","LSMR lambda=0", "LSMR lambda=30","hybrid LSQR","AB-GMRES","BA-GMRES","AB-GMRES FDK","BA-GMRES FDK")) plt.show() # plot images -tigre.plotimg(np.concatenate([imgCGLS, imgSIRT], axis=1), dim="z", step=2) +tigre.plotimg(np.concatenate([np.concatenate([imgCGLS, imgSIRT, imgLSQR,imgabgmres,imgabgmresfdk],axis=1),np.concatenate([imgLSMR, imgLSMR2, imghLSQR,imgbagmres,imgbagmresfdk], axis=1)], axis=2), dim="z", step=2,clims=[0, 2]) # plot errors -tigre.plotimg(np.abs(np.concatenate([head - imgCGLS, head - imgSIRT], axis=1)), dim="z", slice=32) +tigre.plotimg(np.concatenate([np.concatenate([head-imgCGLS, head-imgSIRT, head-imgLSQR, head-imgabgmres, head-imgabgmresfdk],axis=1),np.concatenate([head-imgLSMR, head-imgLSMR2, head-imghLSQR,head-imgbagmres,head-imgbagmresfdk], axis=1)], axis=2), dim="z", slice=32) diff --git a/Python/demos/d09_Algorithms04.py b/Python/demos/d09_Algorithms04.py index 3957b54d..7a0918b5 100644 --- a/Python/demos/d09_Algorithms04.py +++ b/Python/demos/d09_Algorithms04.py @@ -170,8 +170,8 @@ # ========================================================================== # ========================================================================== # -# This is a more edge preserving algorithms than ASD_POCS, in theory. -# delta is the cuttof vlaue of anromalized edge exponential weight.... +# This is a more edge preserving algorithms than ASD_POCS. +# delta is the cut-off value of a normalized edge exponential weight.... # not super clear, but it cotnrols at which point you accept something as real vs noise edge. imgAWASDPOCS = algs.awasd_pocs( @@ -189,9 +189,19 @@ delta=np.array([-0.005]), ) + + +# IRN-CGLS-TV: CGLS with TV regularization +# ========================================================================== +# ========================================================================== +# +# Should converge fairly fast to a TV regularized solution +imgCGLSTV=algs.irn_tv_cgls(noise_projections,geo,angles,10,lmbda=5,niter_outer=2) + + #%% plot results # plot images tigre.plotimg( - np.concatenate([imgAWASDPOCS, imgOSASDPOCS, imgASDPOCS, imgOSSART], axis=1), dim="z", step=2 + np.concatenate([imgAWASDPOCS, imgOSASDPOCS, imgASDPOCS, imgOSSART, imgCGLSTV], axis=1), dim="z", step=2 ) diff --git a/Python/tigre/algorithms/__init__.py b/Python/tigre/algorithms/__init__.py index 1d53f0de..0fe068ca 100644 --- a/Python/tigre/algorithms/__init__.py +++ b/Python/tigre/algorithms/__init__.py @@ -11,6 +11,13 @@ from .ista_algorithms import ista from .iterative_recon_alg import iterativereconalg from .krylov_subspace_algorithms import cgls +from .krylov_subspace_algorithms import lsqr +from .krylov_subspace_algorithms import hybrid_lsqr +from .krylov_subspace_algorithms import lsmr +from .krylov_subspace_algorithms import irn_tv_cgls +from .krylov_subspace_algorithms import hybrid_flsqr_tv +from .krylov_subspace_algorithms import ab_gmres +from .krylov_subspace_algorithms import ba_gmres from .pocs_algorithms import asd_pocs from .pocs_algorithms import os_asd_pocs from .pocs_algorithms import awasd_pocs @@ -41,6 +48,8 @@ "os_aw_pcsd", "fbp", "cgls", + "lsqr", + "lsmr", "fista", "ista", "mlem", diff --git a/Python/tigre/algorithms/iterative_recon_alg.py b/Python/tigre/algorithms/iterative_recon_alg.py index ad288db6..9a6a2d85 100644 --- a/Python/tigre/algorithms/iterative_recon_alg.py +++ b/Python/tigre/algorithms/iterative_recon_alg.py @@ -150,6 +150,8 @@ def __init__(self, proj, geo, angles, niter, **kwargs): name="Iterative Reconstruction", sup_kw_warning=False, gpuids=None, + niter_outer=4, + restart=True ) allowed_keywords = [ "V", @@ -167,7 +169,8 @@ def __init__(self, proj, geo, angles, niter, **kwargs): "tvlambda", "hyper", "fista_p", - "fista_q" + "fista_q", + "niter_outer" ] self.__dict__.update(options) self.__dict__.update(**kwargs) diff --git a/Python/tigre/algorithms/krylov_subspace_algorithms.py b/Python/tigre/algorithms/krylov_subspace_algorithms.py index fd78c90c..f493e74d 100644 --- a/Python/tigre/algorithms/krylov_subspace_algorithms.py +++ b/Python/tigre/algorithms/krylov_subspace_algorithms.py @@ -3,11 +3,14 @@ import time import copy import numpy as np +import scipy as sp import tigre from tigre.algorithms.iterative_recon_alg import IterativeReconAlg from tigre.algorithms.iterative_recon_alg import decorator from tigre.utilities.Atb import Atb from tigre.utilities.Ax import Ax +import tigre.algorithms as algs +import scipy.sparse.linalg if hasattr(time, "perf_counter"): @@ -30,26 +33,10 @@ def __init__(self, proj, geo, angles, niter, **kwargs): # Don't precompute V and W. kwargs.update(dict(W=None, V=None)) kwargs.update(dict(blocksize=angles.shape[0])) - self.log_parameters = False self.re_init_at_iteration = 0 IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) - if self.log_parameters: - parameter_history = {} - iterations = self.niter - parameter_history["alpha"] = np.zeros([iterations], dtype=np.float32) - parameter_history["beta"] = np.zeros([iterations], dtype=np.float32) - parameter_history["gamma"] = np.zeros([iterations], dtype=np.float32) - parameter_history["q_norm"] = np.zeros([iterations], dtype=np.float32) - parameter_history["s_norm"] = np.zeros([iterations], dtype=np.float32) - self.parameter_history = parameter_history - - self.__r__ = self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) - self.__p__ = Atb(self.__r__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) - p_norm = np.linalg.norm(self.__p__.ravel(), 2) - self.__gamma__ = p_norm * p_norm - - def reinitialise_cgls(self): + def initialize_algo(self): self.__r__ = self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) self.__p__ = Atb(self.__r__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) p_norm = np.linalg.norm(self.__p__.ravel(), 2) @@ -57,8 +44,11 @@ def reinitialise_cgls(self): # Overide def run_main_iter(self): + self.l2l = np.zeros((1, self.niter), dtype=np.float32) avgtime = [] + + self.initialize_algo() for i in range(self.niter): if self.verbose: self._estimate_time_until_completion(i) @@ -72,27 +62,20 @@ def run_main_iter(self): self.res += alpha * self.__p__ avgtoc = default_timer() avgtime.append(abs(avgtic - avgtoc)) - for item in self.__dict__: - if ( - isinstance(getattr(self, item), np.ndarray) - and np.isnan(getattr(self, item)).any() - ): - raise ValueError("nan found for " + item + " at iteraton " + str(i)) - - aux = self.proj - tigre.Ax( - self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids - ) - self.l2l[0, i] = np.linalg.norm(aux) + + self.l2l[0, i] = np.linalg.norm(self.proj - tigre.Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids)) if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + self.res -= alpha * self.__p__ + if self.verbose: print("re-initilization of CGLS called at iteration:" + str(i)) - if self.re_init_at_iteration + 1 == i: - if self.verbose: - print("Algorithm exited with two consecutive reinitializations.") + if self.re_init_at_iteration + 1 == i or not self.restart: + print("CGLS exited due to divergence.") return self.res - self.res -= alpha * self.__p__ - self.reinitialise_cgls() - self.re_init_at_iteration = i + self.re_init_at_iteration=i + i=i-1 + self.initialize_algo() + break self.__r__ -= alpha * q s = tigre.Atb(self.__r__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) @@ -100,12 +83,6 @@ def run_main_iter(self): gamma1 = s_norm * s_norm beta = gamma1 / self.__gamma__ - if self.log_parameters: - self.parameter_history["alpha"][i] = alpha - self.parameter_history["beta"][i] = beta - self.parameter_history["gamma"][i] = self.__gamma__ - self.parameter_history["q_norm"][i] = q_norm - self.parameter_history["s_norm"][i] = s_norm self.__gamma__ = gamma1 self.__p__ = s + beta * self.__p__ @@ -119,5 +96,765 @@ def run_main_iter(self): + "(s)" ) - cgls = decorator(CGLS, name="cgls") + +class LSQR(IterativeReconAlg): + __doc__ = ( + " LSQR solves the CBCT problem using the least squares\n" + " LSQR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + def initialize_algo(self): + # Paige and Saunders //doi.org/10.1145/355984.355989 + + # Enumeration as given in the paper for 'Algorithm LSQR' + # (1) Initialize + self.__u__=self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + + normr = np.linalg.norm(self.__u__.ravel(), 2) + self.__u__ = self.__u__/normr + + self.__beta__ = normr + self.__phibar__ = normr + self.__v__ = Atb(self.__u__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) + + self.__alpha__ =np.linalg.norm(self.__v__.ravel(), 2) + self.__v__ = self.__v__/self.__alpha__ + self.__rhobar__ = self.__alpha__ + self.__w__ = np.copy(self.__v__) + + def run_main_iter(self): + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + avgtime = [] + self.initialize_algo() + for i in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(i) + if self.Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + avgtic = default_timer() + + #% (3)(a) + self.__u__ = tigre.Ax(self.__v__, self.geo, self.angles, "Siddon", gpuids=self.gpuids) - self.__alpha__*self.__u__ + self.__beta__ = np.linalg.norm(self.__u__.ravel(),2) + self.__u__ = self.__u__ / self.__beta__ + + #% (3)(b) + self.__v__ = tigre.Atb(self.__u__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) - self.__beta__*self.__v__ + self.__alpha__ = np.linalg.norm(self.__v__.ravel(),2) + self.__v__ = self.__v__ / self.__alpha__ + + #% (4)(a-g) + rho = np.sqrt(self.__rhobar__**2 + self.__beta__**2) + c = self.__rhobar__ / rho + s = self.__beta__ / rho + theta = s * self.__alpha__ + self.__rhobar__ = - c * self.__alpha__ + phi = c * self.__phibar__ + self.__phibar__ = s * self.__phibar__ + + #% (5) Update x, w + self.res = self.res + (phi / rho) * self.__w__ + self.__w__ = self.__v__ - (theta / rho) * self.__w__ + + avgtoc = default_timer() + avgtime.append(abs(avgtic - avgtoc)) + + if self.Quameasopts is not None: + self.error_measurement(res_prev, i) + + self.l2l[0, i] = np.linalg.norm(self.proj - tigre.Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids)) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + self.res -= (phi / rho) * (self.__v__-self.__w__)/((theta / rho)) + if self.verbose: + print("re-initilization of LSQR called at iteration:" + str(i)) + if self.re_init_at_iteration + 1 == i or not self.restart: + print("LSQR exited due to divergence.") + return self.res + self.re_init_at_iteration=i + i=i-1 + self.initialize_algo() + break + +lsqr = decorator(LSQR, name="lsqr") + +class hybrid_LSQR(IterativeReconAlg): + __doc__ = ( + " hybrid_LSQR solves the CBCT problem using the hybrid_LSQR\n" + " hybrid_LSQR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + self.__U__ = np.zeros((self.niter+1,np.prod(self.geo.nDetector)*len(self.angles)),dtype=np.float32) + self.__V__ = np.zeros((self.niter,(np.prod(self.geo.nVoxel))),dtype=np.float32) + self.__B__ = np.zeros((self.niter,self.niter+1),dtype=np.float32) #% Projected matrix + self.__proj_rhs__ = np.zeros((self.niter+1,1),dtype=np.float32) #% Projected right hand side + + def initialize_algo(self): + # Paige and Saunders //doi.org/10.1145/355984.355989 + # Enumeration as given in the paper for 'Algorithm LSQR' + # % Initialise matrices + + # (1) Initialize + self.__u__=self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + + normr = np.linalg.norm(self.__u__.ravel(), 2) + self.__u__ = self.__u__/normr + self.__U__[0]=self.__u__.ravel() + + self.__beta__ = normr + self.__proj_rhs__[0]=normr + + + def run_main_iter(self): + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + avgtime = [] + self.initialize_algo() + for i in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(i) + + avgtic = default_timer() + + v = Atb(self.__u__,self.geo,self.angles,backprojection_type="matched",gpuids=self.gpuids) + + if i>0: + v = np.reshape(v.ravel() - self.__beta__*self.__V__[i-1],v.shape) + + + for j in range(i-1): + v=np.reshape(v.ravel()-(self.__V__[j]*v.ravel())*self.__V__[j],v.shape) + + + alpha = np.linalg.norm(v.ravel(), 2) + v = v/alpha + self.__V__[i] = v.ravel() + + #% Update U_{ii+1} + self.__u__ = tigre.Ax(v, self.geo, self.angles, "Siddon", gpuids=self.gpuids) - alpha*self.__u__ + + for j in range(i-1): + self.__u__=np.reshape(self.__u__.ravel()-(self.__U__[j]*self.__u__.ravel())*self.__U__[j],self.__u__.shape) + + + self.__beta__ = np.linalg.norm(self.__u__.ravel(), 2) + self.__u__ = self.__u__ / self.__beta__ + self.__U__[i+1] = self.__u__.ravel() + + #% Update projected matrix + self.__B__[i,i] = alpha + self.__B__[i,i+1] = self.__beta__ + #% Malena. Proposed update: we should check algorithms breaks; + #% 'if abs(alpha) <= eps || abs(beta) <= eps' - end and save + + #% Solve the projected problem + #% (using the SVD of the small projected matrix) + Bk = self.__B__[0:i+1,0:i+2] + Uk, Sk, Vk = np.linalg.svd(np.transpose(Bk)) + + if i==0: + Sk = Sk[0] + + rhsk = self.__proj_rhs__[0:i+2] + rhskhat = np.matmul(np.transpose(Uk),rhsk) # + Dk = Sk**2 + self.lmbda**2 + + rhskhat = Sk * rhskhat[0:i+1,0] + yhat = rhskhat[0:i+1]/Dk + y = np.matmul(np.transpose(Vk), yhat) + + + + + self.l2l[0, i] = np.linalg.norm(self.proj - tigre.Ax(self.res + np.reshape(np.matmul(np.transpose(self.__V__[0:i+1]),y),self.res.shape), self.geo, self.angles, "Siddon", gpuids=self.gpuids)) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + if self.re_init_at_iteration + 1 == i or not self.restart: + print("hybrid LSQR exited due to divergence at iteration "+str(i)) + return self.res + np.reshape(np.matmul(np.transpose(self.__V__[0:i+1]),y),self.res.shape) + + #% Test for convergence. + #% msl: I still need to implement this. + #% msl: There are suggestions on the original paper. Let's talk about it! + + self.res = self.res + np.reshape(np.matmul(np.transpose(self.__V__),y),self.res.shape) + return self.res +hybrid_lsqr = decorator(hybrid_LSQR, name="hybrid_lsqr") + +class LSMR(IterativeReconAlg): + __doc__ = ( + " LSMR solves the CBCT problem using LSMR\n" + " LSMR(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + def initialize_algo(self): + #% David Chin-Lung Fong and Michael Saunders //doi.org/10.1137/10079687X + #% Enumeration as given in the paper for 'Algorithm LSMR' + #% (1) Initialize + self.__u__=self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + normr = np.linalg.norm(self.__u__.ravel(), 2) + self.__beta__ = normr + self.__u__ = self.__u__/normr + + self.__v__ = Atb(self.__u__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) + self.__alpha__ =np.linalg.norm(self.__v__.ravel(), 2) + self.__v__ = self.__v__/self.__alpha__ + + self.__alphabar__ = self.__alpha__ + self.__zetabar__ = self.__alpha__ * self.__beta__ + self.__rho__ = 1 + self.__rhobar__ = 1 + self.__cbar__ = 1 + self.__sbar__ = 0 + self.__h__ = self.__v__ + self.__hbar__ = 0 + + #% Compute the residual norm ||r_k|| + self.__betadd__ = self.__beta__ + self.__betad__ = 0 + self.__rhod__ = 1 + self.__tautilda__ = 0 + self.__thetatilda__ = 0 + self.__zeta__ = 0 + self.__d__ = 0 + + def run_main_iter(self): + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + avgtime = [] + self.initialize_algo() + for i in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(i) + if self.Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + avgtic = default_timer() + + #% (3) Continue the bidiagonalization + self.__u__ = tigre.Ax(self.__v__, self.geo, self.angles, "Siddon", gpuids=self.gpuids) - self.__alpha__*self.__u__ + self.__beta__ = np.linalg.norm(self.__u__.ravel(),2) + self.__u__ = self.__u__ / self.__beta__ + + self.__v__ = tigre.Atb(self.__u__, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) - self.__beta__*self.__v__ + self.__alpha__ = np.linalg.norm(self.__v__.ravel(),2) + self.__v__ = self.__v__ / self.__alpha__ + + #% (4) Construct and apply rotation \hat{P}_k + alphahat = np.sqrt(self.__alphabar__**2 + self.lmbda**2) + chat = self.__alphabar__/alphahat + shat = self.lmbda/alphahat + + #% (5) Construct and apply rotation P_k + rhopre = self.__rho__; + self.__rho__ = np.sqrt(alphahat**2 + self.__beta__**2) + c = alphahat / self.__rho__ + s = self.__beta__ / self.__rho__ + theta = s * self.__alpha__ + self.__alphabar__ = c * self.__alpha__ + + #% (6) Construct and apply rotation \bar{P}_k + thetabar = self.__sbar__ * self.__rho__ + rhobarpre = self.__rhobar__ + self.__rhobar__ = np.sqrt((self.__cbar__ *self.__rho__)**2 + theta**2) + self.__cbar__ = self.__cbar__ * self.__rho__ / self.__rhobar__ + self.__sbar__ = theta / self.__rhobar__ + zetapre = self.__zeta__ + self.__zeta__ = self.__cbar__ * self.__zetabar__ + self.__zetabar__ = -self.__sbar__ * self.__zetabar__ + + #% (7) Update \bar{h}, x, h + self.__hbar__ = self.__h__ - (thetabar*self.__rho__/(rhopre*rhobarpre))*self.__hbar__ + self.res = self.res + (self.__zeta__ / (self.__rho__*self.__rhobar__)) * self.__hbar__ + self.__h__ = self.__v__ - (theta / self.__rho__) * self.__h__ + + #% (8) Apply rotation \hat{P}_k, P_k + betaacute = chat* self.__betadd__ + betacheck = - shat* self.__betadd__ + + #% Computing ||r_k|| + + betahat = c * betaacute + betadd = -s * betaacute + + #% Update estimated quantities of interest. + #% (9) If k >= 2, construct and apply \tilde{P} to compute ||r_k|| + rhotilda = np.sqrt(self.__rhod__**2 + thetabar**2) + ctilda = self.__rhod__ / rhotilda + stilda = thetabar / rhotilda + thetatildapre = self.__thetatilda__ + self.__thetatilda__ = stilda * self.__rhobar__ + self.__rhod__ = ctilda * self.__rhobar__ + #% betatilda = ctilda * betad + stilda * betahat; % msl: in the orinal paper, but not used + self.__betad__ = -stilda * self.__betad__ + ctilda * betahat + + #% (10) Update \tilde{t}_k by forward substitution + self.__tautilda__ = (zetapre - thetatildapre* self.__tautilda__) / rhotilda + taud = (self.__zeta__ - self.__thetatilda__*self.__tautilda__) / self.__rhod__ + + #% (11) Compute ||r_k|| + self.__d__ = self.__d__ + betacheck**2 + gamma_var = self.__d__ + (self.__betad__ - taud)**2 + betadd**2 + + avgtoc = default_timer() + avgtime.append(abs(avgtic - avgtoc)) + + if self.Quameasopts is not None: + self.error_measurement(res_prev, i) + + self.l2l[0, i] = np.linalg.norm(self.proj - tigre.Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids)) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + self.res -= (self.__zeta__ / (self.__rho__*self.__rhobar__)) * self.__hbar__ + if self.verbose: + print("re-initilization of LSMR called at iteration:" + str(i)) + if self.re_init_at_iteration + 1 == i or not self.restart: + print("LSMR exited due to divergence.") + return self.res + self.re_init_at_iteration=iter + iter=iter-1 + self.initialize_algo() + break +lsmr = decorator(LSMR, name="lsmr") + +class IRN_TV_CGLS(IterativeReconAlg): + __doc__ = ( + " IRN_TV_CGLS solves the CBCT problem using CGLS with TV constraints\n" + " IRN_TV_CGLS(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + + def __build_weights__(self): + Dxx=np.copy(self.res) + Dyx=np.copy(self.res) + Dzx=np.copy(self.res) + + Dxx[0:-2,:,:]=self.res[0:-2,:,:]-self.res[1:-1,:,:] + Dyx[:,0:-2,:]=self.res[:,0:-2,:]-self.res[:,1:-1,:] + Dzx[:,:,0:-2]=self.res[:,:,0:-2]-self.res[:,:,1:-1] + + return (Dxx**2+Dyx**2+Dzx**2+1e-6)**(-1/4) + + def Lx(self,W,img): + Dxx=np.copy(img) + Dyx=np.copy(img) + Dzx=np.copy(img) + + Dxx[0:-2,:,:]=img[0:-2,:,:]-img[1:-1,:,:] + Dyx[:,0:-2,:]=img[:,0:-2,:]-img[:,1:-1,:] + Dzx[:,:,0:-2]=img[:,:,0:-2]-img[:,:,1:-1] + + return np.stack((W*Dxx,W*Dyx,W*Dzx),axis=0) + + def Ltx(self,W,img3): + Wx_1 = W * img3[0,:,:,:] + Wx_2 = W * img3[1,:,:,:] + Wx_3 = W * img3[2,:,:,:] + + DxtWx_1=Wx_1 + DytWx_2=Wx_2 + DztWx_3=Wx_3 + + DxtWx_1[1:-2,:,:]=Wx_1[1:-2,:,:]-Wx_1[0:-3,:,:] + DxtWx_1[-1,:,:]=-Wx_1[-2,:,:] + + DytWx_2[:,1:-2,:]=Wx_2[:,1:-2,:]-Wx_2[:,0:-3,:] + DytWx_2[:,-1,:]=-Wx_2[:,-2,:] + + DztWx_3[:,:,1:-2]=Wx_3[:,:,1:-2]-Wx_3[:,:,0:-3] + DztWx_3[:,:,-1]=-Wx_3[:,:,-2] + + return DxtWx_1 + DytWx_2 + DztWx_3 + + def run_main_iter(self): + self.l2l = np.zeros((1, self.niter*self.niter_outer), dtype=np.float32) + avgtime = [] + + res0=self.res + + for outer in range(self.niter_outer): + if self.verbose: + niter=self.niter + self.niter=self.niter_outer + self._estimate_time_until_completion(outer) + self.niter=niter + if self.Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + avgtic = default_timer() + + + W=self.__build_weights__() + self.res=res0 + + prox_aux_1 =Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + prox_aux_2 = self.Lx(W,self.res)*np.sqrt(self.lmbda) + + r_aux_1 = self.proj - prox_aux_1 + r_aux_2 = -prox_aux_2 + #% Malena: changed the format, r_aux_2 is 3 + #% r = cat(3,r_aux_1, r_aux_2); % Malena: size guide, erase later, N x N x (100 + N-1) + p_aux_1 = tigre.Atb(r_aux_1, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) + p_aux_2 = np.sqrt(self.lmbda)*self.Ltx(W, r_aux_2) + p = p_aux_1 + p_aux_2 + + gamma=np.linalg.norm(p.ravel(),2)**2 + + for i in range(self.niter): + res0=self.res + + q_aux_1 = tigre.Ax(p, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + q_aux_2 = self.Lx(W,p)*np.sqrt(self.lmbda) + + #% q = cat(3, q_aux_1, q_aux_2{1},q_aux_2{2},q_aux_2{3}); % Probably never need to actually do this + #% alpha=gamma/norm(q(:),2)^2; + alpha=gamma/(np.linalg.norm(q_aux_1.ravel(),2)**2 + np.linalg.norm(q_aux_2[0].ravel(),2)**2 + np.linalg.norm(q_aux_2[1].ravel(),2)**2+np.linalg.norm(q_aux_2[2].ravel(),2)**2) + self.res=self.res+alpha*p + aux=self.proj-tigre.Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + #% residual norm or the original least squares (not Tikhonov). + #% Think if that is what we want of the NE residual + self.l2l[0, outer*self.niter+i] = np.linalg.norm(aux.ravel(),2) + + + #% If step is adecuate, then continue withg CGLS + r_aux_1 = r_aux_1-alpha*q_aux_1 + r_aux_2=r_aux_2-alpha*q_aux_2 + + s_aux_1 = tigre.Atb(r_aux_1, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) + s_aux_2 = np.sqrt(self.lmbda) * self.Ltx(W, r_aux_2) + s = s_aux_1 + s_aux_2 + + gamma1=np.linalg.norm(s.ravel(),2)**2 + beta=gamma1/gamma + gamma=gamma1 + p=s+beta*p + + avgtoc = default_timer() + avgtime.append(abs(avgtic - avgtoc)) + + if self.Quameasopts is not None: + self.error_measurement(res_prev, outer) + +irn_tv_cgls = decorator(IRN_TV_CGLS, name="irn_tv_cgls") + + +class AB_GMRES(IterativeReconAlg): + __doc__ = ( + " AB_GMRES solves the CBCT problem using preconditioned GMRES\n" + " AB_GMRES(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + if "backprojector" in kwargs: + backproject=kwargs.pop("backprojector") + else: + backproject="matched" + if backproject == "matched": + self.backproject=Atb + elif backproject == "FDK": + self.backproject=algs.fdk + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + + def __compute_res__(self,x,w,y): + y=y.astype(np.float32) + for i in range(w.shape[0]): + x=x+self.backproject(np.reshape(w[i],self.proj.shape),self.geo,self.angles,gpuids=self.gpuids)*y[i] + return x + + def run_main_iter(self): + + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + w = np.zeros((self.niter+1,np.prod(self.geo.nDetector)*len(self.angles)),dtype=np.float32) + r=self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + w[0] = r.ravel()/np.linalg.norm(r.ravel(), 2) + h=np.zeros((self.niter,self.niter+1),dtype=np.float32) + for k in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(k) + + qk=Ax(self.backproject(np.reshape(w[k],self.proj.shape),self.geo,self.angles,gpuids=self.gpuids),self.geo, self.angles, "Siddon", gpuids=self.gpuids) + e1=np.zeros(k+2) + e1[0]=1 + for i in range(k+1): + h[k,i]=sum(qk.ravel()*w[i]) + qk=qk.ravel()-h[k,i]*w[i] + + h[k,k+1]=np.linalg.norm(qk.ravel(),2) + w[k+1]=qk.ravel()/h[k,k+1] + y=np.linalg.lstsq(np.transpose(h[0:k+1,0:k+2]),e1*np.linalg.norm(r.ravel(),2),rcond=None) + y=y[0] + self.l2l[0, i] = np.linalg.norm((self.proj - tigre.Ax(self.__compute_res__(self.res,w[0:k+1],y),self.geo,self.angles, "Siddon",gpuids=self.gpuids)).ravel(),2) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + if self.re_init_at_iteration + 1 == i or not self.restart: + print("AB-GMRES exited due to divergence at iteration "+str(i)) + return self.__compute_res__(self.res,w[0:k+1],y) + + self.res=self.__compute_res__(self.res,w[0:-1],y) + return self.res +ab_gmres = decorator(AB_GMRES, name="ab_gmres") + + + +class BA_GMRES(IterativeReconAlg): + __doc__ = ( + " BA_GMRES solves the CBCT problem using preconditioned GMRES\n" + " AB_GMRES(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + if "backprojector" in kwargs: + backproject=kwargs.pop("backprojector") + else: + backproject="matched" + if backproject == "matched": + self.backproject=Atb + elif backproject == "FDK": + self.backproject=algs.fdk + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + def __compute_res__(self,x,w,y): + y=y.astype(np.float32) + for i in range(w.shape[0]): + x=x+np.reshape(w[i],self.res.shape)*y[i] + return x + + def run_main_iter(self): + + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + w = np.zeros((self.niter+1,(np.prod(self.geo.nVoxel))),dtype=np.float32) + r=self.backproject(self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids), self.geo, self.angles, gpuids=self.gpuids) + w[0] = r.ravel()/np.linalg.norm(r.ravel(), 2) + h=np.zeros((self.niter,self.niter+1),dtype=np.float32) + + for k in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(k) + + qk=self.backproject(Ax(np.reshape(w[k],self.res.shape),self.geo,self.angles, "Siddon",gpuids=self.gpuids),self.geo, self.angles, gpuids=self.gpuids) + e1=np.zeros(k+2) + e1[0]=1 + for i in range(k+1): + h[k,i]=sum(qk.ravel()*w[i]) + qk=qk.ravel()-h[k,i]*w[i] + + h[k,k+1]=np.linalg.norm(qk.ravel(),2) + w[k+1]=qk.ravel()/h[k,k+1] + y=np.linalg.lstsq(np.transpose(h[0:k+1,0:k+2]),e1*np.linalg.norm(r.ravel(),2),rcond=None) + y=y[0] + + self.l2l[0, i] = np.linalg.norm((self.proj - tigre.Ax(self.__compute_res__(self.res,w[0:k+1],y), self.geo, self.angles, "Siddon", gpuids=self.gpuids)).ravel(),2) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + if self.re_init_at_iteration + 1 == i or not self.restart: + print("BA-GMRES exited due to divergence at iteration "+str(i)) + return self.__compute_res__(self.res,w[0:k+1],y) + + self.res=self.__compute_res__(self.res,w[0:-1],y) + return self.res +ba_gmres = decorator(BA_GMRES, name="ba_gmres") + + +class hybrid_fLSQR_TV(IterativeReconAlg): + __doc__ = ( + " hybrid_fLSQR_TV solves the CBCT problem using preconditioned hybrid flexuble LSQR with TV regularization\n" + " AB_GMRES(PROJ,GEO,ANGLES,NITER) solves the reconstruction problem\n" + " using the projection data PROJ taken over ALPHA angles, corresponding\n" + " to the geometry descrived in GEO, using NITER iterations." + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + # Don't precompute V and W. + kwargs.update(dict(W=None, V=None)) + kwargs.update(dict(blocksize=angles.shape[0])) + self.re_init_at_iteration = 0 + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + self.__U__ = np.zeros((self.niter+1,np.prod(self.geo.nDetector)*len(self.angles)),dtype=np.float32) + self.__V__ = np.zeros((self.niter,(np.prod(self.geo.nVoxel))),dtype=np.float32) + self.__Z__ = np.zeros((self.niter,(np.prod(self.geo.nVoxel))),dtype=np.float32) + + self.__M__ = np.zeros((self.niter,self.niter+1),dtype=np.float32) #% Projected matrix + self.__T__ = np.zeros((self.niter,self.niter),dtype=np.float32) #% Projected matrix + self.__proj_rhs__ = np.zeros((self.niter+1,1),dtype=np.float32) #% Projected right hand side + + + + def __build_weights__(self): + Dxx=np.copy(self.res) + Dyx=np.copy(self.res) + Dzx=np.copy(self.res) + + Dxx[0:-2,:,:]=self.res[0:-2,:,:]-self.res[1:-1,:,:] + Dyx[:,0:-2,:]=self.res[:,0:-2,:]-self.res[:,1:-1,:] + Dzx[:,:,0:-2]=self.res[:,:,0:-2]-self.res[:,:,1:-1] + + return (Dxx**2+Dyx**2+Dzx**2+1e-6)**(-1/4) + + def Lx(self,W,img): + img=np.reshape(img,self.res.shape) + Dxx=np.copy(img) + Dyx=np.copy(img) + Dzx=np.copy(img) + + Dxx[0:-2,:,:]=img[0:-2,:,:]-img[1:-1,:,:] + Dyx[:,0:-2,:]=img[:,0:-2,:]-img[:,1:-1,:] + Dzx[:,:,0:-2]=img[:,:,0:-2]-img[:,:,1:-1] + + return np.stack((W*Dxx,W*Dyx,W*Dzx),axis=0) + + def Ltx(self,W,img3): + img3 =np.reshape(img3,(3,*self.res.shape)) + Wx_1 = W * img3[0,:,:,:] + Wx_2 = W * img3[1,:,:,:] + Wx_3 = W * img3[2,:,:,:] + + DxtWx_1=Wx_1 + DytWx_2=Wx_2 + DztWx_3=Wx_3 + + DxtWx_1[1:-2,:,:]=Wx_1[1:-2,:,:]-Wx_1[0:-3,:,:] + DxtWx_1[-1,:,:]=-Wx_1[-2,:,:] + + DytWx_2[:,1:-2,:]=Wx_2[:,1:-2,:]-Wx_2[:,0:-3,:] + DytWx_2[:,-1,:]=-Wx_2[:,-2,:] + + DztWx_3[:,:,1:-2]=Wx_3[:,:,1:-2]-Wx_3[:,:,0:-3] + DztWx_3[:,:,-1]=-Wx_3[:,:,-2] + + return DxtWx_1 + DytWx_2 + DztWx_3 + + def mvT(self,k_aux,x ): + return x.ravel() - k_aux.ravel()*np.sum(x.ravel()) + def mv(self,k_aux,x ): + return x.ravel() - np.dot(k_aux.ravel(),x.ravel()) + + def run_main_iter(self): + # % Initialise matrices + + # (1) Initialize + u=self.proj - Ax(self.res, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + + k_aux = Ax(np.ones(self.res.shape,dtype=np.float32)/np.sqrt(np.prod(self.geo.nVoxel)), self.geo, self.angles, "Siddon", gpuids=self.gpuids) + xA0 = 1/(np.sqrt(np.prod(self.geo.nVoxel))*np.linalg.norm(k_aux.ravel(),2)**2)*np.dot(k_aux.ravel(),u.ravel()) + xA0 =np.ones(self.res.shape,dtype=np.float32)*xA0 + + k_aux = 1/(np.sqrt(np.prod(self.geo.nVoxel))*np.linalg.norm(k_aux.ravel(),2)**2)*Atb(k_aux, self.geo, self.angles, backprojection_type="matched", gpuids=self.gpuids) + + u = u - Ax(xA0, self.geo, self.angles, "Siddon", gpuids=self.gpuids) + + normr=np.linalg.norm(u.ravel(),2) + u=u/normr + self.__U__[0]=u.ravel() + + self.__proj_rhs__[0]=normr + + if np.max(self.res)==0: + W=np.ones(self.res.shape,dtype=np.float32) + else: + W=self.__build_weights__() + + self.l2l = np.zeros((1, self.niter), dtype=np.float32) + + for i in range(self.niter): + if self.verbose: + self._estimate_time_until_completion(i) + v=Atb(u,self.geo,self.angles,backprojection_type="matched", gpuids=self.gpuids) + for j in range(i+1): + self.__T__[i,j] = np.dot(self.__V__[j],v.ravel()) + v = v.ravel() - self.__T__[i,j]*self.__V__[j] + v=np.reshape(v,self.res.shape) + self.__T__[i,i] = np.linalg.norm(v.ravel(),2) + v=v/self.__T__[i,i] + + z=self.mvT(k_aux,v) + + L = sp.sparse.linalg.LinearOperator((np.prod(self.res.shape),np.prod(self.res.shape)*3), matvec=lambda x: self.Ltx(W,x).ravel(),rmatvec=lambda x: self.Lx(W,x).ravel()) + aux_z=sp.sparse.linalg.lsqr(L,z.ravel(),iter_lim=50) + L = sp.sparse.linalg.LinearOperator((np.prod(self.res.shape)*3,np.prod(self.res.shape)), matvec=lambda x: self.Lx(W,x).ravel(),rmatvec=lambda x: self.Ltx(W,x).ravel()) + z=sp.sparse.linalg.lsqr(L,aux_z[0].ravel(),iter_lim=50) + z=z[0].astype(np.float32) + + z=self.mv(k_aux,z) + + self.__V__[i]=v.ravel() + self.__Z__[i]=z.ravel() + + # Update U and projected matrix M + u = Ax(np.reshape(z,self.res.shape), self.geo, self.angles, "Siddon", gpuids=self.gpuids) + for j in range(i+1): + self.__M__[i,j] = np.dot(self.__U__[j],u.ravel()) + u = u.ravel() - np.dot(self.__M__[i,j],self.__U__[j]) + self.__M__[i,i+1]=np.linalg.norm(u.ravel(),2) + u=u/self.__M__[i,i+1] + u=np.reshape(u,self.proj.shape) + self.__U__[i+1]=u.ravel() + + ### Solve the regularized projected problem + # (using the DVF of the small projected matrix) + + Mk=self.__M__[0:i+1,0:i+2] + # Prepare the projected regularization term + WZ=np.zeros((i+1,(3*np.prod(self.geo.nVoxel))),dtype=np.float32) + for j in range(i+1): + # This can be more efficiently... + # DZ can be saved and updated at each iteration + WZ[j]=self.Lx(W,self.__Z__[j]).ravel() + __,ZRk =sp.linalg.qr(np.transpose(WZ),mode='economic') + ZRksq=ZRk[0:i+1,0:i+1] + rhsk=self.__proj_rhs__[0:i+2] + MZk=np.concatenate((np.transpose(Mk),self.lmbda*ZRksq)) + rhsZk=np.concatenate((rhsk,np.zeros((i+1,1),dtype=np.float32))) + y = np.linalg.lstsq(MZk, rhsZk,rcond=None) + print(y[0]) + d=np.matmul(np.transpose(self.__Z__[0:i+1]),y[0]) + + x = self.res + np.reshape(d,self.res.shape) + xA0 + self.l2l[0, i] = np.linalg.norm((self.proj - tigre.Ax(x, self.geo, self.angles, "Siddon", gpuids=self.gpuids)).ravel(),2) + if i > 0 and self.l2l[0, i] > self.l2l[0, i - 1]: + if self.re_init_at_iteration + 1 == i or not self.restart: + print("BA-GMRES exited due to divergence at iteration "+str(i)) + return x + + return x + +hybrid_flsqr_tv = decorator(hybrid_fLSQR_TV, name="hybrid_flsqr_tv") diff --git a/README.md b/README.md index b5bb454e..5fd56725 100644 --- a/README.md +++ b/README.md @@ -51,13 +51,13 @@ TIGRE is a GPU-based CT reconstruction software repository that contains a wide - **Iterative algorithms** - - Gradient-based algorithms (SART, OS-SART, SIRT) with multiple tuning parameters (Nesterov acceleration, initialization, parameter reduction, ...) + - Gradient-based algorithms (SART, OS-SART, SIRT, ASD-POCS, OS-ASD-POCS, B-ASD-POCS-β, PCSD, AwPCSD, Aw-ASD-POCS) with multiple tuning parameters (Nesterov acceleration, initialization, parameter reduction, ...) - - Krylov subspace algorithms (CGLS) + - Krylov subspace algorithms (CGLS, LSQR, hybrid LSQR, LSMR, IRN-TV-CGLS, hybrid-fLSQR-TV, AB/BA-GMRES) - Statistical reconstruction (MLEM) - - Total variation regularization based algorithms: proximal-based (FISTA, SART-TV) and POCS-based (ASD-POCS, OS-ASD-POCS, B-ASD-POCS-β, PCSD, AwPCSD, Aw-ASD-POCS) + - Variational methods (FISTA, SART-TV) - TV denoising for 3D images.