Fix issue with offset between Absorption and Grain reconstructions

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

5_reconstruction/gtAstraPrepareAbsorptionStack.m

@@ -73,11 +73,13 @@ imgs_nums = 0 : rec_abs.interval : totproj-1;
 Omega_rad = angle_rad(imgs_nums + 1);
 Omega_deg = Omega_rad * 180 / pi;
 
+%bb = [acq.bb(1), detgeo.detrefv - acq.bb(4) / 2 + 0.5, acq.bb(3), acq.bb(4)];
+
 bbpos = repmat(acq.bb, [numel(Omega_deg), 1]);
 % bbpos = [gtGeoGrainCenterSam2Det(samgeo.orig, 0, parameters) - (acq.bb(3:4) - 1) / 2, acq.bb(3:4)];
 % bbpos = repmat(bbpos, [numel(Omega_deg), 1]);
 [~, gr_shift] = gtMatchGetSampleShifts(parameters, Omega_deg);
-
+extra_detector_offset = acq.rotu - (detgeo.detrefu - detgeo.detrefpos(2) / detgeo.pixelsizeu);  %sub-pixel shift for 2D-FBP algorithm
 geom = gtGeoProjForReconstruction([], Omega_deg', [], bbpos, [], ...
     detgeo, labgeo, samgeo, recgeo, 'ASTRA_absorption', gr_shift);
 if (~isempty(conf.rot_angles))
@@ -117,7 +119,7 @@ for jj = 1:numel(imgs_nums)
 
     fprintf(repmat('\b', 1, num_chars))
 end
-stack = -log(max(0.01, stack));
+stack = -log(max(0.001, stack));
 
 fprintf('Done in %f seconds.\n', toc(c));
 
@@ -132,7 +134,7 @@ fprintf('...done!\n')
 % Dealing with the padding
 new_stack_size = stack_size + [2, 0, 2] * rec_abs.padding;
 new_stack_shifts = [1, 0, 1] * rec_abs.padding;
-stack(stack < 0) = 0;
+%stack(stack < 0) = 0;
 new_stack = gtPlaceSubVolume(zeros(new_stack_size, 'single'), stack, new_stack_shifts);
 new_stack(1:rec_abs.padding, :) = new_stack((rec_abs.padding + 1) * ones(rec_abs.padding, 1), :);
 new_stack((end - rec_abs.padding + 1):end, :) = new_stack((end - rec_abs.padding) * ones(rec_abs.padding, 1), :);
@@ -142,6 +144,7 @@ proj = struct( ...
     'psf', rec_abs.psf, ...
     'rot_angles', conf.rot_angles, ... % Additional rotation
     'rot_axes', conf.rot_axes, ...   % Additional rotation
+    'extra_detector_offset', extra_detector_offset, ... % Correction for sub-pixel shift of rotation axis in 2D-FBP algorithm
     'stack', new_stack, ...
     'Omega', Omega_rad, ...
     'geom', geom, ...

5_reconstruction/gtAstra_FBP2D.m

@@ -34,10 +34,7 @@ function [absvol, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_off
     hsize = size(proj.stack, 1);  % number of pixels in projection image
     nproj = size(proj.stack, 2);  % number of projections
 
-    hsize_zp = ceil(hsize * sqrt(2));   % zeropad projection images
-    offset = ceil((hsize_zp - hsize) / 2);
-
-    sino = zeros(nproj, hsize_zp);
+    sino = zeros(nproj, hsize);
 
     xvol = hsize;
     yvol = hsize;
@@ -46,13 +43,13 @@ 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_zp, proj.Omega);
+    proj_geom = astra_create_proj_geom('parallel', 1.0, hsize, proj.Omega);
     omegas = abs(proj.Omega) / pi * 180;
     [shift_lab, shift_sam] = gtMatchGetSampleShifts(parameters, omegas);
     proj_geom.option.ExtraDetectorOffset = shift_lab(:,2)' ./ parameters.detgeo.pixelsizeu;
-    if exist('extra_detector_offset', 'var')
+    if isfield(proj, 'extra_detector_offset')
         disp('shifting projections');
-        proj_geom.option.ExtraDetectorOffset = proj_geom.option.ExtraDetectorOffset + ones(1, nproj) * extra_detector_offset;
+        proj_geom.option.ExtraDetectorOffset = proj_geom.option.ExtraDetectorOffset + ones(1, nproj) * proj.extra_detector_offset;
     end
 
     vol_geom = astra_create_vol_geom(hsize, hsize);
@@ -77,11 +74,12 @@ function [absvol, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_off
         % Not sure if we really need to re-create the algorithm each time..
         alg_id = astra_mex_algorithm('create', cfg);
 
-        % 'zeropadding' projections with the value
-        sino = gtPlaceSubImage(squeeze(proj_data(:, :, ii))', sino, offset, 1);
-
-        sino(:, 1:offset-1) = repmat(sino(:, offset), 1, offset-1);
-        sino(:, offset+hsize:end) = repmat(sino(:, offset+hsize-1), 1, hsize_zp-offset-hsize+1);
+        % 'zeropadding' projections is already dealt with in
+        % gtAstraPrepareAbsorptionStack
+        %sino = gtPlaceSubImage(squeeze(proj_data(:, :, ii))', sino, offset, 1);
+        sino = proj_data(:, :, ii)';
+        %sino(:, 1:offset-1) = repmat(sino(:, offset), 1, offset-1);
+        %sino(:, offset+hsize:end) = repmat(sino(:, offset+hsize-1), 1, hsize_zp-offset-hsize+1);
 
         astra_mex_data2d('set', sinogram_id, sino)
 

zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m

@@ -35,7 +35,17 @@ function [proj, verts] = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, om
     % Calculate the projection geometry for the spots:
     %  - Vector used for ROI grain reconstruction in ASTRA (grain shifted
     %    to center of roi volume)
-    [gr_shift_lab, gr_shift_sam] = gtMatchGetSampleShifts(parameters, omegas);
+    if isfield(parameters.acq, 'correct_sample_shifts')
+        if parameters.acq.correct_sample_shifts
+            [gr_shift_lab, gr_shift_sam] = gtMatchGetSampleShifts(parameters, omegas);
+        else
+            gr_shift_lab = [];
+            gr_shift_sam = [];
+        end
+    else
+        gr_shift_lab = [];
+        gr_shift_sam = [];
+    end
     proj_geom = gtGeoProjForReconstruction(diff_beam_dirs, ...
         omegas, gr_center, bbpos_det_grain, [], ...
         detgeo, labgeo, samgeo, recgeo, 'ASTRA_grain', gr_shift_sam);
@@ -65,12 +75,12 @@ function [proj, verts] = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, om
     proj.centerpix = gtGeoSam2Sam(gr_center, samgeo, recgeo, true);
 
     use_polyhedron = ~cloning & (numel(find(selected)) >= 3);
-    if (use_polyhedron)
+    if 0 %(use_polyhedron)
         % This should behave better with vertical detector
         % (if no bug is introduced with it :D)
         verts = gtFwdSimComputeCircumscribingPolyhedron(...
             gr_center, diff_beam_dirs(selected, :), omegas(selected), ...
-            bb(selected, :), parameters, det_ind);
+            bb(selected, :), gr_shift_lab(selected, :), parameters, det_ind, 'convhull', false);
 
         vol_size = 2 * max(abs(max(verts, [], 1)), abs(min(verts, [], 1)));
     else

zUtil_Geo/gtGeoDetDefaultParameters.m

@@ -89,10 +89,12 @@ function d_yz = sfDetrefpos(acq_bb, par_labgeo)
 % account for an offset of the sample bounding box position:
 %   Y lateral: It assumes that the projection of the rotation axis is at the
 %   center of the sample bounding box.
-%   Z vertical: set in a way that the vertical center of the sample bounding
-%   box will have coordinate Z=0.
-% It assumes rotpos = [0 0 0] and view into the beam.
+%   Z vertical: By definition we set the detector center as 0 position
+%
+%   NB: Obsolete definition:  Vertical center of the sample bounding (beam center) corresponds to Z=0.
+%   It assumes rotpos = [0 0 0] and view into the beam.
 
-    d_yz(1) = -((acq_bb(1) + acq_bb(3)/2) - par_labgeo.detrefu) * par_labgeo.pixelsizeu;
-    d_yz(2) =  ((acq_bb(2) + acq_bb(4)/2) - par_labgeo.detrefv) * par_labgeo.pixelsizev;
+    d_yz(1) = -((acq_bb(1) + acq_bb(3)/2) - par_labgeo.detrefu - 0.5) * par_labgeo.pixelsizeu;
+    %d_yz(2) =  ((acq_bb(2) + acq_bb(4)/2) - par_labgeo.detrefv - 0.5) * par_labgeo.pixelsizev;
+    d_yz(2) = 0;  % By definition
 end
\ No newline at end of file

zUtil_Geo/gtGeoProjForReconstruction.m

@@ -105,18 +105,21 @@ function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,.
     end
     projvec_rec = gtMathsNormalizeVectorsList(projvec_rec);
 
-    if (strcmpi(rectype, 'ASTRA_grain') || strcmpi(rectype, 'ASTRA_absorption'))
-        % Volume center in REC reference
-        if (strcmpi(rectype, 'ASTRA_absorption'))
-            vol_center_rec = recgeo.orig;
-        else
+
+    switch rectype
+        case 'ASTRA_absorption'
+            % Volume center in REC reference
+            samorig_rec = gtGeoSam2Sam(samgeo.orig, samgeo, recgeo, false, false);
+            proj_shift  = samorig_rec(ones_omega, :) - vol_center_shift_rec;
+            bbcent_rec  = bbcent_rec + proj_shift;
+
+        case 'ASTRA_grain'
             vol_center_rec = gtGeoSam2Sam(vol_center, samgeo, recgeo, false, false);
-        end
-        vol_center_rec = vol_center_rec(ones_omega, :) + vol_center_shift_rec;
-        bbcent_rec = bbcent_rec - vol_center_rec;
-    else
-        % case {'ASTRA_full'}
-        bbcent_rec = bbcent_rec - vol_center_shift_rec;
+            proj_shift     = -vol_center_rec(ones_omega, :) - vol_center_shift_rec;
+            bbcent_rec     = bbcent_rec + proj_shift;
+
+        case 'ASTRA_full'
+            bbcent_rec = bbcent_rec - vol_center_shift_rec;
     end
 
     proj_geom = [projvec_rec, bbcent_rec, detdiru_rec, detdirv_rec];

zUtil_Geo/gtGeoRecDefaultParameters.m

-function par_recgeo = gtGeoRecDefaultParameters(detgeo, samgeo)
+function par_recgeo = gtGeoRecDefaultParameters(detgeo, samgeo, acq)
 % GTGEOSAMDEFAULTPARAMETERS  Stores and returns the default Sample geomtery
 %                            parameters
 %
-%     par_recgeo = gtGeoRecDefaultParameters([detgeo [, samgeo]])
+%     par_recgeo = gtGeoRecDefaultParameters([detgeo [, samgeo, acq]])
 %     ---------------------------------------
 %     The default Reconstruction reference is set as the Lab but rotating with the
 %     sample. Voxel size is determined from detector pixelsize.
-%
+%     If acq is provided, the vertical position is offset to match the illuminated region defined by acq.bb
+%     otherwise it is set to 0 (center of detector)
 
     if (exist('samgeo', 'var') && ~isempty(samgeo))
         par_recgeo.orig = samgeo.orig;
@@ -26,5 +27,9 @@ function par_recgeo = gtGeoRecDefaultParameters(detgeo, samgeo)
         det_pixel_size = 1e-3;
     end
 
+    if (exist('acq', 'var') && ~isempty(acq))
+        par_recgeo.orig(3) = -(acq.bb(2) + acq.bb(4) / 2 - detgeo.detrefv - 0.5) * detgeo.pixelsizev;
+    end
+
     par_recgeo.voxsize = det_pixel_size([1 1 1]);
 end

zUtil_Parameters/build_list_v2.m

@@ -610,7 +610,8 @@ list.fsim(end+1, :) = {'sel_int_stdev', ...
     'Selects spots for reconsruction that are in the top 4 sigma of the intensity rankings', 'double', 1};
 list.fsim(end+1, :) = {'sel_avint_stdev', ...
     'Selects spots for reconsruction that are in the top 3.5 sigma of the average intensity per pixel rankings', 'double', 1};
-
+list.fsim(end+1, :) = {'mode', ...
+    'Mode of operation {''indexter''}/''farfield''/''cloning''', 'char', 1};
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % reconstruction