Added objective function, improved output, and added references.

demo_polyquant_cbct.m

@@ -19,7 +19,7 @@
 % o Adjust regularisation and convergence parameters an compare to fanbeam.
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 % Created:      26/04/2019
-% Last edit:    29/05/2019
+% Last edit:    31/05/2019
 % Jonathan Hugh Mason
 %
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -45,11 +45,14 @@
 
 %% Polyquant setup and reconstruction
 mode = [];
+mode.useConst = true;
 mode.verbose = 2;
 mode.tau = 2;
 mode.nSplit = 32;
 mode.maxIter = 300;
 mode.scatFun = scatEst;
 lambda = 2;  % can be optimised for better results
 mode.proxFun = @(z,t) prox_tv3d_nn(z,t*lambda);
-out = polyquant(mode,head.specData,head.proj,i0,A,head.eden);
+mode.regFun = @(z) norm_tv3d(z);
+out = polyquant(mode,head.specData,head.proj,i0,A,head.eden);
+fprintf('Reconstructed with PSNR = %.2f dB\n',20*log10(max(head.eden(:))./out.rmse(end)));

demo_polyquant_fanbeam.m

@@ -35,7 +35,7 @@
 %   equivalent to a linearised model (monoenergetic).
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 % Created:      26/04/2019
-% Last edit:    29/05/2019
+% Last edit:    31/05/2019
 % Jonathan Hugh Mason
 %
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -57,15 +57,18 @@
 %% Polyquant setup 
 sample = 'chest';  % 'brain', 'head', 'abdomen', 'pelvis', 'implant'
 mode = [];
+mode.useConst = true;  % just to offset objective function to better range
 mode.contrast = [0.5,1.4];  % display contrast
-mode.verbose = 2;
+mode.verbose = 2;  % 0 = no output; 1 = text output; 2 = text+image output
 mode.tau = 5;  % conservative choice
 mode.nSplit = 24;
 mode.maxIter = 500;  % more iterations recommended for implant
 lambda = 0.5;  % can be optimised for better results
 mode.numLinFit = 2;  % can be set to 2 for tissue, but use 3 for implant
 mode.proxFun = @(z,t) prox_tv_nn(z,t*lambda);
+mode.regFun = @(z) norm_tv(z);
 xTrue = mat_to_den(attenuationDb,single(fandata.(sample).mat));
 
 %% Perform the reconstruction
-out = polyquant(mode,fandata.specData,fandata.(sample).proj,fandata.i0,A,xTrue);
+out = polyquant(mode,fandata.specData,fandata.(sample).proj,fandata.i0,A,xTrue);
+fprintf('Reconstructed with PSNR = %.2f dB\n',20*log10(max(xTrue(:))./out.rmse(end)));

demo_polyquant_scatter.m

@@ -25,7 +25,7 @@
 % o Adjust regularisation and convergence parameters an compare to fanbeam.
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 % Created:      26/04/2019
-% Last edit:    29/05/2019
+% Last edit:    31/05/2019
 % Jonathan Hugh Mason
 %
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -51,6 +51,7 @@
 A = Gcone(cg, ig, 'type', 'sf2', 'class', 'Fatrix');
 %% Polyquant setup and reconstruction
 mode = [];
+mode.useConst = true;
 mode.verbose = 2;
 mode.tau = 5;  % a little more aggressive
 mode.offset = true;
@@ -60,6 +61,7 @@
 mode.numLinFit = 2;  % since there's no metal implants
 lambda = 0.1;  % can be optimised for better results
 mode.proxFun = @(z,t) prox_tv3d_nn(z,t*lambda);
+mode.regFun = @(z) norm_tv3d(z);
 angArray = ((-90:2.25:270-2.25));
 %% Scatter estimation function
 % The factor 1.5 on the incident intensity is to compensate for the 
@@ -69,4 +71,5 @@
                  poly_sks(1.5*i0,projA,projB,projC,rho,angArray(subSet),...
                  pelvis.specData,scatParam,32,cg,ig,[0.3,15]);
 %% Perfrm Polyquant reconstrution
-out = polyquant(mode,pelvis.specData,pelvis.proj,i0,A,pelvis.eden);
+out = polyquant(mode,pelvis.specData,pelvis.proj,i0,A,pelvis.eden);
+fprintf('Reconstructed with PSNR = %.2f dB\n',20*log10(max(pelvis.eden(:))./out.rmse(end)));

polyquant.m

@@ -11,7 +11,7 @@
 % (xTrue)       -- ground truth image (optional)
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 % Created:      07/03/2018
-% Last edit:    26/04/2019
+% Last edit:    31/05/2019
 % Jonathan Hugh Mason
 %
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -49,8 +49,14 @@
 if ~isfield(mode,'L')
     mode.L = lipscitz_estimate(specData,I0,mode.scat,y,Ab*x0,Af);
 end
-%
+
 alpha = mode.nSplit*mode.tau/mode.L;  % the step-size
+if mode.useConst
+    const = y-y.*log(y+eps);
+    const = sum(const(:));  % a constant offset for objective function
+else
+    const = 0;
+end
 
 x1 = x0;
 timeTot = tic;
@@ -59,7 +65,7 @@
    t = 1;
 end
 
-out.res(1) = rms(x1(:)-xTrue(:));
+out.rmse(1) = rms(x1(:)-xTrue(:));
 if mode.verbose == 2
     if ndims(xTrue) == 3
         subplot(2,3,1),imshow(imrotate(xTrue(:,:,20),-90),mode.contrast);
@@ -72,8 +78,11 @@
     drawnow;
 end
 grAx = @(x1,is,ys,ind,subSet) polyquant_grad(specData,A,At,is,x1,ys,ind,mode.scatFun,subSet,w);
-
+objFac = zeros(size(y)); out.scat = zeros(size(y));
 %% The main iterative loop
+if mode.verbose > 0
+    fprintf('Starting Polyquant reconstruction:\n');
+end
 for k = 1:mode.maxIter
     ind = mod(k,mode.nSplit)+1;
     if mode.bitRev
@@ -90,17 +99,23 @@
     end
 
     gradAx = grAx(x1,is,ys,ind,subSet);
-    out.scat(:,:,subSet) = gradAx.s;
+    if ndims(x0) == 3
+        out.scat(:,:,subSet) = gradAx.s;
+        objFac(:,:,subSet) = gradAx.objFac;
+    else
+        out.scat(:,subSet) = gradAx.s;
+        objFac(:,subSet) = gradAx.objFac;
+    end
     xNew = mode.proxFun(x1-alpha*gradAx.grad,alpha);
 
     t1 = 0.5*(1+sqrt(1+4*t^2));
     x1 = xNew+(t-1)/t1*(xNew-x0);
     x0 = xNew;
     t = t1;
-
-    out.res(k+1) = rms(x1(:)-xTrue(:));
+    out.rmse(k+1) = rms(x1(:)-xTrue(:));
+    out.obj(k+1) = sum(double(objFac(:)+out.scat(:)-y(:).*log(objFac(:)+out.scat(:)+eps)))-const+mode.regFun(x1);
     if mode.verbose > 0
-      fprintf('\r  Iter = %d; subset = %d; res = %d        ',k,ind,out.res(k+1));
+      fprintf('\rIter = %i;\t RMSE = %.4e;\t obj = %.4e;\t subset = %i    ',k,out.rmse(k+1),out.obj(k+1),ind);
     end
     if mode.verbose == 2
         str = ['polyquant at iteration: ',num2str(k)];
@@ -117,9 +132,9 @@
 end
 time = toc(timeTot);
 out.time = time;
-out.rec = xNew;
+out.recon = xNew;
 if mode.verbose > 0
-    fprintf('\n Finished in %d seconds\n',time);
+    fprintf('\n Finished in %.2e seconds\n',time);
 end
 
 end
@@ -172,6 +187,7 @@
 end
 
 strOut.grad = out;
+strOut.objFac = mainFac;
 strOut.s = s;
 end
 
@@ -209,6 +225,7 @@
 if ~isfield(mode,'regFun'),     mode.regFun = @(z) 0; end
 if ~isfield(mode,'proxFun'),    mode.proxFun = @(z,t) prox_nz(z); end
 if ~isfield(mode,'contrast'),   mode.contrast = [0,2]; end
+if ~isfield(mode,'useConst'),   mode.useConst = false; end
 if ~isfield(mode,'scatFun')
     mode.scat = 0;
     mode.scatFun = @(z,~,~,~,~,~,~) 0; 
@@ -221,6 +238,7 @@
 end
 
 function out = offset_weight(proj,cg)
+    % Offset weighting for half-fan case from [G. Wang, Med Phys. 2002]
     out = proj;
     us = ((cg.ns/2-0.5):-1:(-cg.ns/2+0.5))*cg.ds - cg.offset_s*cg.ds;
     overlap = max(us);

utilities/bit_rev.m

@@ -1,4 +1,5 @@
 function out = bit_rev(in,sz)
+% Bit-reversal ordering as in [G. Herman and L. Meyer, IEEE TMI 1993]
 facs = factor(sz);
 unit = zeros(size(facs));
 unit(1) = 1;