Merge pull request #416 from twhitbread/FixFDKOffsetBug

Use original geometry for FDK redundancy weighting

Author: AnderBiguri

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);
-    % 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

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)

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
 

MATLAB/Utilities/redundancy_weighting.m

@@ -6,13 +6,14 @@
 if ~isfield(geo,'COR')
     geo.COR=0;
 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 = (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)]);