From 074c9e3e393cb2aa3394b88e282c2232ce163cfd Mon Sep 17 00:00:00 2001 From: twhitbread Date: Fri, 2 Dec 2022 10:51:24 +0000 Subject: [PATCH 1/5] Use original geometry for FDK redundancy weighting --- MATLAB/Algorithms/FDK.m | 2 +- MATLAB/Utilities/redundancy_weighting.m | 11 +++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/MATLAB/Algorithms/FDK.m b/MATLAB/Algorithms/FDK.m index d0832e6f..01670ed1 100644 --- a/MATLAB/Algorithms/FDK.m +++ b/MATLAB/Algorithms/FDK.m @@ -48,7 +48,7 @@ % for all cases, not just wang? [zproj, zgeo] = zeropadding(proj, geo); % Preweighting using Wang function - proj=zproj.*redundancy_weighting(zgeo); + proj=zproj.*redundancy_weighting(zgeo, geo); % Replace original proj and geo % proj = proj_w; geo = zgeo; diff --git a/MATLAB/Utilities/redundancy_weighting.m b/MATLAB/Utilities/redundancy_weighting.m index 81223faf..125565c4 100644 --- a/MATLAB/Utilities/redundancy_weighting.m +++ b/MATLAB/Utilities/redundancy_weighting.m @@ -1,4 +1,4 @@ -function w = redundancy_weighting(geo) +function w = redundancy_weighting(geo, originalGeo) % Preweighting using Wang function % Ref: % Wang, Ge. X-ray micro-CT with a displaced detector array. Medical Physics, 2002,29(7):1634-1636. @@ -6,17 +6,20 @@ if ~isfield(geo,'COR') geo.COR=0; end +if nargin == 1 + originalGeo = geo; +end offset = geo.offDetector(1); offset = offset + (geo.DSD(1) / geo.DSO(1)) * geo.COR(1); % added correction us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) + abs(offset); us = us * geo.DSO(1)/geo.DSD(1); -theta = (geo.sDetector(1)/2 - abs(offset))... - * sign(offset); +theta = (originalGeo.sDetector(1)/2 - abs(originalGeo.offDetector(1)))... + * sign(originalGeo.offDetector(1)); abstheta = abs(theta * geo.DSO(1)/geo.DSD(1)); w = ones([geo.nDetector(2),geo.nDetector(1)]); -if apply_wang_weights(geo) +if apply_wang_weights(originalGeo) for ii = 1:geo.nDetector t = us(ii); if(abs(t) <= abstheta) From fad215e2f01a7d50ad9d3a31f8c948f7e01efd4b Mon Sep 17 00:00:00 2001 From: twhitbread Date: Fri, 2 Dec 2022 13:44:04 +0000 Subject: [PATCH 2/5] Fix error message --- MATLAB/Utilities/Atb.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MATLAB/Utilities/Atb.m b/MATLAB/Utilities/Atb.m index b21c6785..daa21396 100644 --- a/MATLAB/Utilities/Atb.m +++ b/MATLAB/Utilities/Atb.m @@ -54,7 +54,7 @@ end %% geometry geo=checkGeo(geo,angles); -assert(isequal([size(projections,2) size(projections,1)],geo.nDetector.'),'TIGRE:checkGeo:BadGeometry','nVoxel does not match with provided image size'); +assert(isequal([size(projections,2) size(projections,1)],geo.nDetector.'),'TIGRE:checkGeo:BadGeometry','nDetector does not match with provided image size'); %% Thats it, lets call the mex fucntion From 2b5774d6f215ce5db0dc6f3821999a36b4c449af Mon Sep 17 00:00:00 2001 From: twhitbread Date: Wed, 7 Dec 2022 17:11:02 +0000 Subject: [PATCH 3/5] Revert bug and correct typos --- MATLAB/Algorithms/MLEM.m | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/MATLAB/Algorithms/MLEM.m b/MATLAB/Algorithms/MLEM.m index 6a803c8e..4c3d8542 100644 --- a/MATLAB/Algorithms/MLEM.m +++ b/MATLAB/Algorithms/MLEM.m @@ -1,22 +1,23 @@ function [res,qualMeasOut]=MLEM(proj,geo,angles,niter,varargin) -%MLEM solves the tomographic problem by using Maximum Likelihood Expection -% Maximitation algorithm. +%MLEM solves the tomographic problem by using Maximum Likelihood Expectation +% Maximisation algorithm. % % MLEM(PROJ,GEO,ALPHA,NITER,opt) solves the reconstruction problem % using the projection data PROJ taken over ALPHA angles, corresponding -% to the geometry descrived in GEO, using NITER iterations. +% to the geometry described in GEO, using NITER iterations. % -% 'verbose': get feedback or not. Default: 1 +% 'verbose': Get feedback or not. Default: 1 % -% 'init': Describes diferent initialization techniques. +% 'init': Describes different initialization techniques. % • 'none' : Initializes the image to ones (default) -% • 'FDK' : intializes image to FDK reconstrucition +% • 'FDK' : Initializes image to FDK reconstruction % -% 'QualMeas' Asks the algorithm for a set of quality measurement +% 'QualMeas': Asks the algorithm for a set of quality measurement % parameters. Input should contain a cell array of desired % quality measurement names. Example: {'CC','RMSE','MSSIM'} % These will be computed in each iteration. -% 'groundTruth' an image as grounf truth, to be used if quality measures +% +% 'groundTruth': An image as ground truth, to be used if quality measures % are requested, to plot their change w.r.t. this known % data. %-------------------------------------------------------------------------- @@ -43,26 +44,27 @@ 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") + warning("Image metrics requested but none caught as output. Call the algorithm with 3 outputs to store them") measurequality=false; end qualMeasOut=zeros(length(QualMeasOpts),niter); res = max(res,0); -% Projection weight, W -W=computeW(geo,angles,gpuids); +% Back-projection weight, V +V = Atb(ones(size(proj),'single'),geo,angles,'matched','gpuids',gpuids); +V(V<=0.) = inf; for ii=1:niter if measurequality && ~strcmp(QualMeasOpts,'error_norm') - res_prev = res; % only store if necesary + res_prev = res; % only store if necessary end if (ii==1);tic;end den = Ax(res,geo,angles,'gpuids',gpuids); den(den<=0.)=inf; - imgupdate = Atb(proj./den, geo,angles,'matched','gpuids',gpuids)./W; + imgupdate = Atb(proj./den, geo,angles,'matched','gpuids',gpuids)./V; res = max(res.*imgupdate,0.); if measurequality @@ -87,7 +89,7 @@ % Check inputs nVarargs = length(argin); if mod(nVarargs,2) - error('TIGRE:FISTA:InvalidInput','Invalid number of inputs') + error('TIGRE:MLEM:InvalidInput','Invalid number of inputs') end % check if option has been passed as input @@ -103,7 +105,7 @@ for ii=1:length(opts) opt=opts{ii}; default=defaults(ii); - % if one option isnot default, then extranc value from input + % if one option is not default, then extract value from input if default==0 ind=double.empty(0,1);jj=1; while isempty(ind) From 019788548c2bbfe232db5f0eb7e17aa0c9d7de1e Mon Sep 17 00:00:00 2001 From: Ander Biguri Date: Wed, 14 Dec 2022 15:53:03 +0000 Subject: [PATCH 4/5] Removed optional parameter to redundacy_weight --- MATLAB/Algorithms/FDK.m | 8 +++----- MATLAB/Utilities/redundancy_weighting.m | 12 +++++------- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/MATLAB/Algorithms/FDK.m b/MATLAB/Algorithms/FDK.m index 01670ed1..67e418e5 100644 --- a/MATLAB/Algorithms/FDK.m +++ b/MATLAB/Algorithms/FDK.m @@ -46,12 +46,10 @@ if dowang % Zero-padding to avoid FFT-induced aliasing %TODO: should't this be % for all cases, not just wang? - [zproj, zgeo] = zeropadding(proj, geo); % Preweighting using Wang function - proj=zproj.*redundancy_weighting(zgeo, geo); - % Replace original proj and geo - % proj = proj_w; - geo = zgeo; + proj=proj.*redundancy_weighting(geo); + [proj, geo] = zeropadding(proj, geo); + end if size(geo.offDetector,2)==1 diff --git a/MATLAB/Utilities/redundancy_weighting.m b/MATLAB/Utilities/redundancy_weighting.m index 125565c4..4ae34467 100644 --- a/MATLAB/Utilities/redundancy_weighting.m +++ b/MATLAB/Utilities/redundancy_weighting.m @@ -1,4 +1,4 @@ -function w = redundancy_weighting(geo, originalGeo) +function w = redundancy_weighting(geo) % Preweighting using Wang function % Ref: % Wang, Ge. X-ray micro-CT with a displaced detector array. Medical Physics, 2002,29(7):1634-1636. @@ -6,20 +6,18 @@ if ~isfield(geo,'COR') geo.COR=0; end -if nargin == 1 - originalGeo = geo; -end + offset = geo.offDetector(1); offset = offset + (geo.DSD(1) / geo.DSO(1)) * geo.COR(1); % added correction us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) + abs(offset); us = us * geo.DSO(1)/geo.DSD(1); -theta = (originalGeo.sDetector(1)/2 - abs(originalGeo.offDetector(1)))... - * sign(originalGeo.offDetector(1)); +theta = (end.sDetector(1)/2 - abs(end.offDetector(1)))... + * sign(end.offDetector(1)); abstheta = abs(theta * geo.DSO(1)/geo.DSD(1)); w = ones([geo.nDetector(2),geo.nDetector(1)]); -if apply_wang_weights(originalGeo) +if apply_wang_weights(geo) for ii = 1:geo.nDetector t = us(ii); if(abs(t) <= abstheta) From f58fb831748ff9f8c028d5247d771980763eb205 Mon Sep 17 00:00:00 2001 From: Ander Biguri Date: Wed, 14 Dec 2022 15:57:43 +0000 Subject: [PATCH 5/5] Fix typo --- MATLAB/Utilities/redundancy_weighting.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/MATLAB/Utilities/redundancy_weighting.m b/MATLAB/Utilities/redundancy_weighting.m index 4ae34467..bb53d53d 100644 --- a/MATLAB/Utilities/redundancy_weighting.m +++ b/MATLAB/Utilities/redundancy_weighting.m @@ -12,8 +12,8 @@ us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) + abs(offset); us = us * geo.DSO(1)/geo.DSD(1); -theta = (end.sDetector(1)/2 - abs(end.offDetector(1)))... - * sign(end.offDetector(1)); +theta = (geo.sDetector(1)/2 - abs(geo.offDetector(1)))... + * sign(geo.offDetector(1)); abstheta = abs(theta * geo.DSO(1)/geo.DSD(1)); w = ones([geo.nDetector(2),geo.nDetector(1)]);