support for extra-detector offset / sample-shifting in Absorption Reconconstruction

Signed-off-by: Wolfgang Ludwig <wolfgang.ludwig@esrf.fr>

5_reconstruction/gtAstraPrepareAbsorptionStack.m

@@ -25,7 +25,8 @@ conf = struct( ...
     'rot_axes', [0 0 1], ...
     'psf', [], ...
     'save', true, ...
-    'filename', 'absorption_stack' );
+    'filename', 'absorption_stack', ...
+    'extra_detector_offset', []);
 conf = parse_pv_pairs(conf, varargin);
 
 acq = parameters.acq;
@@ -80,15 +81,21 @@ bbpos = repmat(acq.bb, [numel(Omega_deg), 1]);
 % bbpos = repmat(bbpos, [numel(Omega_deg), 1]);
 if isfield(parameters.acq, 'correct_sample_shifts')
     if parameters.acq.correct_sample_shifts
-        [~, shifts_sam] = gtMatchGetSampleShifts(parameters, om);
+        [shifts_lab, shifts_sam] = gtMatchGetSampleShifts(parameters, Omega_deg);
     else
         shifts_sam = zeros(numel(Omega_deg), 3);
+        shifts_lab = zeros(numel(Omega_deg), 3);
     end
 else
     shifts_sam = zeros(numel(Omega_deg), 3);
+    shifts_lab = zeros(numel(Omega_deg), 3);
+end
+%figure; scatter(shifts_lab(:, 1), shifts_lab(:, 2));
+if isempty(conf.extra_detector_offset)
+    extra_detector_offset = (detgeo.detrefu - detgeo.detrefpos(2) / detgeo.pixelsizeu) - acq.rotu  %sub-pixel shift for 2D-FBP algorithm
+else
+    extra_detector_offset = conf.extra_detector_offset;
 end
-
-extra_detector_offset = (detgeo.detrefu - detgeo.detrefpos(2) / detgeo.pixelsizeu) - acq.rotu;  %sub-pixel shift for 2D-FBP algorithm
 geom = gtGeoProjForReconstruction([], Omega_deg', [], bbpos, [], ...
     detgeo, labgeo, samgeo, recgeo, 'ASTRA_absorption', shifts_sam);
 if (~isempty(conf.rot_angles))
@@ -104,8 +111,6 @@ if (~isempty(conf.rot_angles))
     geom = gtMathsMatrixProduct(geom, rot_tensor_tot);
     geom = reshape(geom, [], 12);
 end
-% This was ust a trick to correct for tilts, but it was probably a bad idea
-% geom(:, 10:12) = round(geom(:, 10:12));
 
 fprintf('\b\b: Done.\nBuilding sinogram: \n');
 c = tic();
@@ -132,11 +137,6 @@ stack = -log(max(0.001, stack));
 
 fprintf('Done in %f seconds.\n', toc(c));
 
-% deal with the sense of the rotation
-if all(labgeo.rotdir == [0 0 1])
-    % rotation axis is anti-clockwise
-    Omega_rad = -Omega_rad;
-end
 
 fprintf('...done!\n')
 
@@ -163,7 +163,7 @@ proj = struct( ...
     'vol_size_x', new_stack_size(1), ...
     'vol_size_y', new_stack_size(1), ...
     'vol_size_z', new_stack_size(3) );
-
+    'num_iter', parameters.rec.absorption.num_iter;
 if (conf.save)
     name = fullfile(acq.dir, [conf.filename '.mat']);
     disp(['Writing sinogram in ' name]);

5_reconstruction/gtAstra_FBP2D.m

@@ -43,11 +43,20 @@ function [absvol, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_off
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % projection geometry (type, #rows, #columns, vectors)
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    proj_geom = astra_create_proj_geom('parallel', 1.0, hsize, proj.Omega);
+    % deal with the sense of the rotation,
+    if all(round(parameters.labgeo.rotdir) == [0 0 1])
+        % rotation axis is anti-clockwise: -Omega_rad will be used for
+        % fbp - astra_create_proj_geom, which does not know about rotdir...
+        Omega_fbp = -proj.Omega;
+    else
+        Omega_fbp = proj.Omega;
+    end
+
+    proj_geom = astra_create_proj_geom('parallel', 1.0, hsize, Omega_fbp);
     omegas = abs(proj.Omega) / pi * 180;
     if isfield(parameters.acq, 'correct_sample_shifts')
         if parameters.acq.correct_sample_shifts
-            [shifts_lab, shifts_sam] = gtMatchGetSampleShifts(parameters, om);
+            [shifts_lab, shifts_sam] = gtMatchGetSampleShifts(parameters, omegas);
         else
             shifts_sam = zeros(numel(omegas), 3);
             shifts_lab = zeros(numel(omegas), 3);
@@ -61,7 +70,7 @@ function [absvol, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_off
         disp('shifting projections');
         proj_geom.option.ExtraDetectorOffset = proj_geom.option.ExtraDetectorOffset + ones(1, nproj) * proj.extra_detector_offset;
     end
-
+    figure; scatter(omegas, proj_geom.option.ExtraDetectorOffset)
     vol_geom = astra_create_vol_geom(hsize, hsize);
     proj_id = astra_create_projector('linear', proj_geom, vol_geom);
     sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, 0);