From 59fbde0754719f20e62fe7e85c890fce1056b75a Mon Sep 17 00:00:00 2001
From: Wolfgang Ludwig <wolfgang.ludwig@esrf.fr>
Date: Tue, 5 Jun 2018 18:36:32 +0200
Subject: [PATCH] Fix issue with offset between Absorption and Grain
reconstructions Signed-off-by: Wolfgang Ludwig <wolfgang.ludwig@esrf.fr>
---
.../gtAstraPrepareAbsorptionStack.m | 9 ++++---
5_reconstruction/gtAstra_FBP2D.m | 22 ++++++++--------
zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m | 16 +++++++++---
zUtil_Geo/gtGeoDetDefaultParameters.m | 12 +++++----
zUtil_Geo/gtGeoProjForReconstruction.m | 25 +++++++++++--------
zUtil_Geo/gtGeoRecDefaultParameters.m | 11 +++++---
zUtil_Parameters/build_list_v2.m | 3 ++-
7 files changed, 60 insertions(+), 38 deletions(-)
diff --git a/5_reconstruction/gtAstraPrepareAbsorptionStack.m b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
index a60e0160..1a1b2464 100644
--- a/5_reconstruction/gtAstraPrepareAbsorptionStack.m
+++ b/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, ...
diff --git a/5_reconstruction/gtAstra_FBP2D.m b/5_reconstruction/gtAstra_FBP2D.m
index bb0aab30..cb5a930a 100644
--- a/5_reconstruction/gtAstra_FBP2D.m
+++ b/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)
diff --git a/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m b/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
index f6011046..d529be66 100644
--- a/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
+++ b/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
diff --git a/zUtil_Geo/gtGeoDetDefaultParameters.m b/zUtil_Geo/gtGeoDetDefaultParameters.m
index 205fe856..ed267cb3 100644
--- a/zUtil_Geo/gtGeoDetDefaultParameters.m
+++ b/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
diff --git a/zUtil_Geo/gtGeoProjForReconstruction.m b/zUtil_Geo/gtGeoProjForReconstruction.m
index c5020451..26e290c3 100644
--- a/zUtil_Geo/gtGeoProjForReconstruction.m
+++ b/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];
diff --git a/zUtil_Geo/gtGeoRecDefaultParameters.m b/zUtil_Geo/gtGeoRecDefaultParameters.m
index 7bfbde78..045cdd2a 100644
--- a/zUtil_Geo/gtGeoRecDefaultParameters.m
+++ b/zUtil_Geo/gtGeoRecDefaultParameters.m
@@ -1,12 +1,13 @@
-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
diff --git a/zUtil_Parameters/build_list_v2.m b/zUtil_Parameters/build_list_v2.m
index 24cd2370..3649440c 100644
--- a/zUtil_Parameters/build_list_v2.m
+++ b/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
--
GitLab