From 4d0791d83b93f55a536263cc644686e28a781988 Mon Sep 17 00:00:00 2001
From: Wolfgang Ludwig <wolfgang.ludwig@esrf.fr>
Date: Wed, 13 Dec 2017 08:22:54 +0100
Subject: [PATCH] Correct for sample displacements (e.g. wrong ftomo expert
setting...)
Signed-off-by: Wolfgang Ludwig <wolfgang.ludwig@esrf.fr>
---
3_pairmatchingGUI/gtMatchFindPairCandidates.m | 12 ++++++++----
3_pairmatchingGUI/gtMatchFitParameters.m | 14 ++++++++------
3_pairmatchingGUI/gtMatchIdentifyTheta.m | 4 ++--
3_pairmatchingGUI/gtMatchInitialize.m | 8 ++++++++
3_pairmatchingGUI/gtMatchRankCandidates.m | 18 ++++++++++++++----
3_pairmatchingGUI/gtMatchSaveResults.m | 4 +++-
3_pairmatchingGUI/gtMatchUpdatePairs.m | 2 +-
4_grains/gtCalculateGrain.m | 6 ++++--
.../gtAstraPrepareAbsorptionStack.m | 7 ++++---
5_reconstruction/gtAstra_FBP2D.m | 16 +++++++++-------
zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m | 15 ++++++++++-----
zUtil_Geo/gtGeoDiffVecInLabPair.m | 13 ++++++++-----
zUtil_Geo/gtGeoDiffVecInSample.m | 2 +-
zUtil_Geo/gtGeoProjForReconstruction.m | 17 ++++++++++++++---
zUtil_Indexter/gtINDEXCreateInputData.m | 2 +-
15 files changed, 95 insertions(+), 45 deletions(-)
diff --git a/3_pairmatchingGUI/gtMatchFindPairCandidates.m b/3_pairmatchingGUI/gtMatchFindPairCandidates.m
index 89dda606..0b272873 100644
--- a/3_pairmatchingGUI/gtMatchFindPairCandidates.m
+++ b/3_pairmatchingGUI/gtMatchFindPairCandidates.m
@@ -19,6 +19,7 @@ spB = handles.spB;
nof_spA = length(handles.randindA);
iw = handles.LastRandomIndex;
state = 1;
+shifts = handles.shifts;
% If matching have already been completed
if (iw >= nof_spA)
@@ -63,7 +64,7 @@ while (ishandle(handles.figure1) && (state == 1))
% Select next index for A from random list where they were shuffled
iw = iw + 1;
indA = handles.randindA(iw);
-
+ shiftA = shifts(round(allsp.centim(indA) + 1), :);
% Check if spot A bbox size is large enough
if ((allsp.bbsizeu(indA) >= mp.minsizeu) && (allsp.bbsizev(indA) >= mp.minsizev))
@@ -116,7 +117,6 @@ while (ishandle(handles.figure1) && (state == 1))
if (~isempty(spC))
nofCs = length(spC.ind);
-
goodcands = false(1,nofCs);
meanerrs = zeros(1,nofCs);
thetas = zeros(1,nofCs);
@@ -128,11 +128,11 @@ while (ishandle(handles.figure1) && (state == 1))
% Loop through all its candidates ...
for ic = 1:nofCs
-
+ shiftB(ic, :) = shifts(round(allsp.centim(spC.ind(ic))) + 1, :);
% Check theta angle and plane family for each candidate
[theta,theta_diff,usedthetaind,theta_OK,diffvec] = gtMatchIdentifyTheta(...
[allsp.centu(indA),allsp.centv(indA)], [spC.centu(ic),spC.centv(ic)],...
- handles.labgeo, handles.UsedThetas, mp.thr_theta_used, mp.thr_theta_scale);
+ handles.labgeo, handles.UsedThetas, mp.thr_theta_used, mp.thr_theta_scale, shiftA, shiftB(ic, :));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If candidate is good
@@ -206,6 +206,10 @@ while (ishandle(handles.figure1) && (state == 1))
pairscand.diffvecX(indA,1:nofgoods) = diffvecsX(goodcands);
pairscand.diffvecY(indA,1:nofgoods) = diffvecsY(goodcands);
pairscand.diffvecZ(indA,1:nofgoods) = diffvecsZ(goodcands);
+ pairscand.shiftA(indA,:) = shiftA;
+ pairscand.shiftBX(indA,1:nofgoods) = shiftB(goodcands, 1);
+ pairscand.shiftBY(indA,1:nofgoods) = shiftB(goodcands, 2);
+ pairscand.shiftBZ(indA,1:nofgoods) = shiftB(goodcands, 3);
end
end
diff --git a/3_pairmatchingGUI/gtMatchFitParameters.m b/3_pairmatchingGUI/gtMatchFitParameters.m
index ec759db7..2ad44d87 100644
--- a/3_pairmatchingGUI/gtMatchFitParameters.m
+++ b/3_pairmatchingGUI/gtMatchFitParameters.m
@@ -92,11 +92,13 @@ if handles.fitting.equalize_hkls
end
-pairs.centAU = handles.pairsfilt.centAU(usepair);
-pairs.centAV = handles.pairsfilt.centAV(usepair);
-pairs.centBU = handles.pairsfilt.centBU(usepair);
-pairs.centBV = handles.pairsfilt.centBV(usepair);
-pairs.hkl = handles.pairsfilt.hkl(usepair,:);
+pairs.centAU = handles.pairsfilt.centAU(usepair);
+pairs.centAV = handles.pairsfilt.centAV(usepair);
+pairs.centBU = handles.pairsfilt.centBU(usepair);
+pairs.centBV = handles.pairsfilt.centBV(usepair);
+pairs.shiftA = handles.pairsfilt.shiftA(usepair,:);
+pairs.shiftB = handles.pairsfilt.shiftB(usepair,:);
+pairs.hkl = handles.pairsfilt.hkl(usepair,:);
@@ -314,7 +316,7 @@ function devs = sfDeviation(newvars,newvarsind,oldvars,labgeobase,cryst,...
% Bragg angles theta via the diffraction vector
[~,~,theta] = gtGeoDiffVecInLabPair([pairs.centAU,pairs.centAV], ...
- [pairs.centBU,pairs.centBV], labgeo);
+ [pairs.centBU,pairs.centBV], labgeo, pairs.shiftA, pairs.shiftB);
% 'B matrix': HKL -> Cartesian transformation
Bmat = gtCrystHKL2CartesianMatrix(latticepar);
diff --git a/3_pairmatchingGUI/gtMatchIdentifyTheta.m b/3_pairmatchingGUI/gtMatchIdentifyTheta.m
index 914a4bcc..d21b2cd7 100644
--- a/3_pairmatchingGUI/gtMatchIdentifyTheta.m
+++ b/3_pairmatchingGUI/gtMatchIdentifyTheta.m
@@ -7,11 +7,11 @@
%
function [theta,theta_diff,thetaind,theta_OK,diffvec] = gtMatchIdentifyTheta(...
- centAUV,centBUV,labgeo,compthetas,thr_theta,thr_theta_scale)
+ centAUV,centBUV,labgeo,compthetas,thr_theta,thr_theta_scale, shiftA, shiftB)
% Diffraction vector and theta
-[diffvec,~,theta] = gtGeoDiffVecInLabPair(centAUV,centBUV,labgeo);
+[diffvec,~,theta] = gtGeoDiffVecInLabPair(centAUV,centBUV,labgeo, shiftA, shiftB);
% Discrepancy in theta
diff = theta - compthetas;
diff --git a/3_pairmatchingGUI/gtMatchInitialize.m b/3_pairmatchingGUI/gtMatchInitialize.m
index 83154b5b..53acfee6 100644
--- a/3_pairmatchingGUI/gtMatchInitialize.m
+++ b/3_pairmatchingGUI/gtMatchInitialize.m
@@ -80,6 +80,14 @@ handles.matchpars = orderfields(handles.matchpars, perm);
parnames = fieldnames(handles.matchpars);
handles.matchparnames = parnames;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% Sample Shifts (correction for wrong setting of FTOMO_PAR
+%%% "update_ref_pos"... leading to sample displacements at reference
+%%% positions
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+handles.shifts = gtMatchGetSampleShifts(parameters);
+
+
% Get parameters description
parlist = build_list_v2();
diff --git a/3_pairmatchingGUI/gtMatchRankCandidates.m b/3_pairmatchingGUI/gtMatchRankCandidates.m
index 7a2fbf8d..d5554bb2 100644
--- a/3_pairmatchingGUI/gtMatchRankCandidates.m
+++ b/3_pairmatchingGUI/gtMatchRankCandidates.m
@@ -24,7 +24,7 @@ drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Put list of candidates in one array; keep already computed info
-% [spotA index, spotB index, mean error, theta, diff vector, thetadiff, usedthetaind]
+% [spotA index, spotB index, mean error, theta, diff vector, thetadiff, usedthetaind, shiftA, shiftB]
candslist = zeros(nofcands,3);
lastind = 0;
@@ -56,8 +56,15 @@ for ii = 1:size(pc.ind,1)
candslist(lastind+1:lastind+pc.num(ii),8) = pc.thetadiff(ii,1:pc.num(ii));
% Copy usedthetaind
- candslist(lastind+1:lastind+pc.num(ii),9) = pc.usedthetaind(ii,1:pc.num(ii));
-
+ candslist(lastind+1:lastind+pc.num(ii),9) = pc.usedthetaind(ii,1:pc.num(ii));
+
+ % Copy shifts
+ candslist(lastind+1:lastind+pc.num(ii),10) = pc.shiftA(ii, 1);
+ candslist(lastind+1:lastind+pc.num(ii),11) = pc.shiftA(ii, 2);
+ candslist(lastind+1:lastind+pc.num(ii),12) = pc.shiftA(ii, 3);
+ candslist(lastind+1:lastind+pc.num(ii),13) = pc.shiftBX(ii,1:pc.num(ii));
+ candslist(lastind+1:lastind+pc.num(ii),14) = pc.shiftBY(ii,1:pc.num(ii));
+ candslist(lastind+1:lastind+pc.num(ii),15) = pc.shiftBZ(ii,1:pc.num(ii));
% Update last index in candslist
lastind = lastind + pc.num(ii);
@@ -91,7 +98,8 @@ handles.pairs.theta = candsranked(pairsind,4);
handles.pairs.diffveclab = candsranked(pairsind,[5 6 7]);
handles.pairs.thetadiff = candsranked(pairsind,8);
handles.pairs.usedthetaind = candsranked(pairsind,9);
-
+handles.pairs.shiftA = candsranked(pairsind,[10 11 12]);
+handles.pairs.shiftB = candsranked(pairsind,[13 14 15]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Add other pair info
@@ -119,6 +127,8 @@ handles.pairs.centimB = handles.allsp.centim(handles.pairs.indB);
handles.pairs.omegaA = (handles.pairs.centimA + handles.pairs.centimB)*90/...
handles.parameters.acq.nproj - 90;
handles.pairs.omegaB = handles.pairs.omegaA + 180;
+handles.pairs.shiftA = handles.shifts(round(handles.pairs.centimA) + 1, :);
+handles.pairs.shiftB = handles.shifts(round(handles.pairs.centimB) + 1, :);
handles.pairs.eta = gtGeoEtaFromDiffVec(handles.pairs.diffveclab,handles.labgeo);
diff --git a/3_pairmatchingGUI/gtMatchSaveResults.m b/3_pairmatchingGUI/gtMatchSaveResults.m
index e63152b6..2799ff01 100644
--- a/3_pairmatchingGUI/gtMatchSaveResults.m
+++ b/3_pairmatchingGUI/gtMatchSaveResults.m
@@ -117,8 +117,10 @@ handles.pairsfilt.diffvecsam = gtGeoLab2Sam(handles.pairsfilt.diffveclab,...
% Difspot A in Sample reference
+% Attention - spotA positions are shifted in order to correct for sample
+% movements
centAlab = gtGeoDet2Lab([handles.pairsfilt.centAU,handles.pairsfilt.centAV],...
- handles.labgeo, 0);
+ handles.labgeo, 0) - handles.pairsfilt.shiftA;
handles.pairsfilt.centsamA = gtGeoLab2Sam(centAlab, handles.pairsfilt.omegaA,...
handles.labgeo, samgeo, 0, 0);
diff --git a/3_pairmatchingGUI/gtMatchUpdatePairs.m b/3_pairmatchingGUI/gtMatchUpdatePairs.m
index b1c2677f..dce08a42 100644
--- a/3_pairmatchingGUI/gtMatchUpdatePairs.m
+++ b/3_pairmatchingGUI/gtMatchUpdatePairs.m
@@ -11,7 +11,7 @@ end
% Diffracton vectors and theta
[handles.pairs.diffveclab,~,handles.pairs.theta] = gtGeoDiffVecInLabPair(...
[handles.pairs.centAU,handles.pairs.centAV], ...
- [handles.pairs.centBU,handles.pairs.centBV], handles.labgeo);
+ [handles.pairs.centBU,handles.pairs.centBV], handles.labgeo, handles.pairs.shiftA, handles.pairs.shiftB);
% Discrepancy in theta
diff --git a/4_grains/gtCalculateGrain.m b/4_grains/gtCalculateGrain.m
index c6513b91..9c018f56 100644
--- a/4_grains/gtCalculateGrain.m
+++ b/4_grains/gtCalculateGrain.m
@@ -305,8 +305,10 @@ eta = gtGeoEtaFromDiffVec(dveclab, labgeo);
dvecsam = gtGeoLab2Sam(dveclab, rot_l2s, labgeo, samgeo, true);
% Let's convert grain center Sam -> Lab
+[~, shifts_sam] = gtMatchGetSampleShifts(parameters, om);
+
rot_s2l = permute(rot_l2s, [2 1 3]);
-gcsam_v = gcsam(ones(1, size(rot_l2s, 3)), :);
+gcsam_v = gcsam(ones(1, size(rot_l2s, 3)), :) + shifts_sam;
gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
% Lorentz factor for the given projections (note that eq. 4.1 from page 36
@@ -417,7 +419,7 @@ function sfDisplayFigure(parameters, app, grain, det_ind)
if (app.color)
ondet = grain.allblobs.detector(1).ondet;
om = grain.allblobs.omega(ondet);
- uvw = grain.allblobs.detector(1).uvw(ondet);
+ uvw = grain.allblobs.detector(1).uvw(ondet, :);
cm = inputwdefault(['Choose the type of colormap: omega coloring '...
'(o) or hkl coloring (h): [o]'], 'o');
diff --git a/5_reconstruction/gtAstraPrepareAbsorptionStack.m b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
index 89e44a2b..a60e0160 100644
--- a/5_reconstruction/gtAstraPrepareAbsorptionStack.m
+++ b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
@@ -68,7 +68,7 @@ end
% shifted such that the axis of rotation is in the center of the images
% The origin of the sample coordinate system is on the rotation axis - in
% the middle of the illuminated area
-fprintf('Building geometry..')
+fprintf('Building geometry..\n')
imgs_nums = 0 : rec_abs.interval : totproj-1;
Omega_rad = angle_rad(imgs_nums + 1);
Omega_deg = Omega_rad * 180 / pi;
@@ -76,9 +76,10 @@ Omega_deg = Omega_rad * 180 / pi;
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);
geom = gtGeoProjForReconstruction([], Omega_deg', [], bbpos, [], ...
- detgeo, labgeo, samgeo, recgeo, 'ASTRA_absorption');
+ detgeo, labgeo, samgeo, recgeo, 'ASTRA_absorption', gr_shift);
if (~isempty(conf.rot_angles))
rot_tensor_tot = eye(3);
for ii_a = 1:numel(conf.rot_angles)
@@ -95,7 +96,7 @@ 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: ');
+fprintf('\b\b: Done.\nBuilding sinogram: \n');
c = tic();
abs_files_dir = fullfile(acq.dir, '1_preprocessing', 'abs');
info = [];
diff --git a/5_reconstruction/gtAstra_FBP2D.m b/5_reconstruction/gtAstra_FBP2D.m
index e32c92bc..bb0aab30 100644
--- a/5_reconstruction/gtAstra_FBP2D.m
+++ b/5_reconstruction/gtAstra_FBP2D.m
@@ -1,6 +1,6 @@
-function [abs, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_offset)
+function [absvol, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_offset)
% GTASTRA_FBP2D
-% abs = gtAstra_FBP2D(parameters, proj)
+% absvol = gtAstra_FBP2D(parameters, proj)
% -------------------------------------
%
%
@@ -47,10 +47,12 @@ function [abs, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_offset
% projection geometry (type, #rows, #columns, vectors)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
proj_geom = astra_create_proj_geom('parallel', 1.0, hsize_zp, 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')
disp('shifting projections');
- proj_geom.option.ExtraDetectorOffset = ones(1, nproj) * extra_detector_offset;
+ proj_geom.option.ExtraDetectorOffset = proj_geom.option.ExtraDetectorOffset + ones(1, nproj) * extra_detector_offset;
end
vol_geom = astra_create_vol_geom(hsize, hsize);
@@ -69,7 +71,7 @@ function [abs, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_offset
cfg.FilterType = 'Ram-Lak';
cfg.option.GPUindex = 0;
- abs = zeros(xvol, yvol, zvol, 'single');
+ absvol = zeros(xvol, yvol, zvol, 'single');
for ii = 1:vsize
% Not sure if we really need to re-create the algorithm each time..
@@ -86,7 +88,7 @@ function [abs, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_offset
astra_mex_algorithm('run', alg_id);
% retrieve reconstruction from ASTRA lib
- abs(:, :, vsize-ii+1) = astra_mex_data2d('get', recon_id);
+ absvol(:, :, vsize-ii+1) = astra_mex_data2d('get', recon_id);
astra_mex_algorithm('delete', alg_id);
end
@@ -102,7 +104,7 @@ function [abs, res_norm] = gtAstra_FBP2D(parameters, proj, extra_detector_offset
res_norm = 1;
if (isfield(proj, 'padding'))
- abs = abs( ...
+ absvol = absvol( ...
proj.padding+1:end-proj.padding, ...
proj.padding+1:end-proj.padding, ...
proj.padding+1:end-proj.padding );
diff --git a/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m b/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
index d719c2b1..f6011046 100644
--- a/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
+++ b/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
@@ -11,7 +11,11 @@ function [proj, verts] = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, om
samgeo = parameters.samgeo;
recgeo = parameters.recgeo(det_ind);
detgeo = parameters.detgeo(det_ind);
-
+ if (isfield(parameters, 'clone') & parameters.clone.active)
+ cloning = true;
+ else
+ cloning = false;
+ end
num_projs = size(diff_beam_dirs, 1);
spot_shifts = gtFwdSimGetStackShifts(stackUSize, stackVSize, bb, true);
@@ -31,21 +35,22 @@ 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);
proj_geom = gtGeoProjForReconstruction(diff_beam_dirs, ...
omegas, gr_center, bbpos_det_grain, [], ...
- detgeo, labgeo, samgeo, recgeo, 'ASTRA_grain');
+ detgeo, labgeo, samgeo, recgeo, 'ASTRA_grain', gr_shift_sam);
% Calculate the projection geometry for the full images:
% - Vector describing full projection geometry (full images, grain at
% nominal position in sample volume)
proj_geom_full = gtGeoProjForReconstruction( ...
diff_beam_dirs, omegas, [], bbpos_det_full, ...
- [], detgeo, labgeo, samgeo, recgeo, 'ASTRA_full');
+ [], detgeo, labgeo, samgeo, recgeo, 'ASTRA_full', gr_shift_sam);
% Geometry for extinction spots
proj_geom_abs = gtGeoProjForReconstruction([], omegas, [], ...
bbpos_det_abs, [], detgeo, labgeo, samgeo, recgeo, ...
- 'ASTRA_absorption');
+ 'ASTRA_absorption', gr_shift_lab);
proj = gtFwdSimProjDefinition('fwd_sim');
@@ -59,7 +64,7 @@ function [proj, verts] = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, om
proj.centerpix = gtGeoSam2Sam(gr_center, samgeo, recgeo, true);
- use_polyhedron = numel(find(selected)) >= 3;
+ use_polyhedron = ~cloning & (numel(find(selected)) >= 3);
if (use_polyhedron)
% This should behave better with vertical detector
% (if no bug is introduced with it :D)
diff --git a/zUtil_Geo/gtGeoDiffVecInLabPair.m b/zUtil_Geo/gtGeoDiffVecInLabPair.m
index 9e2e41fb..7fe50c25 100644
--- a/zUtil_Geo/gtGeoDiffVecInLabPair.m
+++ b/zUtil_Geo/gtGeoDiffVecInLabPair.m
@@ -26,18 +26,21 @@
% labgeo.detdirv - detector V direction (unit row vector)
% labgeo.pixelsizeu - pixel size U (same unit as position values)
% labgeo.pixelsizev - pixel size V (same unit as position values)
-%
+% shiftA - shift of sample in LAB reference at rotation
+% positions corresponding to detUVA (nx3)
+% shiftB - shift of sample in LAB reference at roatation
+% postions corresponding to detUVB (nx3)
% OUTPUT diffvec - normalised diffraction vectors in LAB coordinates;
% directions point towards A; size (n,3)
% labA - Lab coordinates of points A (n,3)
% theta - Bragg angles theta in degrees (n,1)
-%
+% output values account for possible shifts of the sample
-function [diffvec,labA,theta] = gtGeoDiffVecInLabPair(detUVA,detUVB,labgeo)
+function [diffvec,labA,theta] = gtGeoDiffVecInLabPair(detUVA,detUVB,labgeo, shiftA, shiftB)
% Coordinate transformation from Detector to Lab
-labA = gtGeoDet2Lab(detUVA, labgeo, 0);
-labB = gtGeoDet2Lab(detUVB, labgeo, 0);
+labA = gtGeoDet2Lab(detUVA, labgeo, 0) - shiftA;
+labB = gtGeoDet2Lab(detUVB, labgeo, 0) - shiftB;
% Mirror B-s on the rotation axis while staying in the same LAB reference:
labBmirr = gtMathsMirrorPointsOnAxis(labB, labgeo.rotpos, labgeo.rotdir);
diff --git a/zUtil_Geo/gtGeoDiffVecInSample.m b/zUtil_Geo/gtGeoDiffVecInSample.m
index 88f55948..da1bdd4d 100644
--- a/zUtil_Geo/gtGeoDiffVecInSample.m
+++ b/zUtil_Geo/gtGeoDiffVecInSample.m
@@ -51,7 +51,7 @@ function [diffvec, samAXYZ, theta] = gtGeoDiffVecInSample(detAUV, omega, samBXYZ
samAXYZ = gtGeoLab2Sam(labA, omega, labgeo, samgeo, 0, 0);
% Diffraction vectors from points samXYZ to point samA:
- diffvec = samAXYZ - samBXYZ;
+ diffvec = samAXYZ - samBXYZ(ones(size(samAXYZ, 1), 1), :);
% Normalise diffvec-s (seems the fastest way for nx3 vectors):
dnorm = sqrt(diffvec(:, 1).^2 + diffvec(:, 2).^2 + diffvec(:, 3).^2) ;
diff --git a/zUtil_Geo/gtGeoProjForReconstruction.m b/zUtil_Geo/gtGeoProjForReconstruction.m
index 3b73283b..4644cdda 100644
--- a/zUtil_Geo/gtGeoProjForReconstruction.m
+++ b/zUtil_Geo/gtGeoProjForReconstruction.m
@@ -1,5 +1,5 @@
function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,...
- bbpos, bboff, detgeo, labgeo, samgeo, recgeo, rectype)
+ bbpos, bboff, detgeo, labgeo, samgeo, recgeo, rectype, vol_center_shift)
% GTGEOPROJFORRECONSTRUCTION Projection coordinates required for
% reconstructions.
%
@@ -34,6 +34,9 @@ function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,.
% recgeo = REConstruction reference as in parameters.recgeo
% rectype = <char> reconstruction geometry:
% 'ASTRA_grain', 'ASTRA_absorption' or 'ASTRA_full'
+% vol_center_shift = shift of grain center in SAMPLE reference (nx3)
+% accounting for sample movements relative to the
+% the rotation axis (see also gtMatchGetSampleShifts)
%
% OUTPUT:
% proj_geom = input for the specified reconstruction geometry:
@@ -52,6 +55,7 @@ function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,.
detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
end
+
% Bounding box center U,V coordinates
bbcent_det = [ ...
bbpos(:, 1) + (bbpos(:, 3) - 1) / 2 + bboff(:, 1) , ...
@@ -68,6 +72,12 @@ function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,.
ones_omega = ones(length(omega), 1);
end
+ if (~exist('vol_center_shift', 'var'))
+ vol_center_shift = zeros(numel(ones_omega), 3);
+ else
+ vol_center_shift_rec = gtGeoSam2Sam(vol_center_shift, samgeo, recgeo, false, false);
+ end
+
% Bounding box center X,Y,Z coordinates in LAB and REC reference
bbcent_lab = gtGeoDet2Lab(bbcent_det, detgeo, 0);
bbcent_rec = gtGeoLab2Sam(bbcent_lab, rot_omega, labgeo, recgeo, 0, 0);
@@ -88,6 +98,7 @@ 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'))
@@ -95,11 +106,11 @@ function proj_geom = gtGeoProjForReconstruction(projvec_sam, omega, vol_center,.
else
vol_center_rec = gtGeoSam2Sam(vol_center, samgeo, recgeo, false, false);
end
- vol_center_rec = vol_center_rec(ones_omega, :);
-
+ 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;
end
proj_geom = [projvec_rec, bbcent_rec, detdiru_rec, detdirv_rec];
diff --git a/zUtil_Indexter/gtINDEXCreateInputData.m b/zUtil_Indexter/gtINDEXCreateInputData.m
index 00a08c74..dbe902cf 100644
--- a/zUtil_Indexter/gtINDEXCreateInputData.m
+++ b/zUtil_Indexter/gtINDEXCreateInputData.m
@@ -62,7 +62,7 @@ cryst.ACM = gtIndexAngularConsMatrix(cryst); % in degrees
% Average pixel size (needed for calculating distance tolerances)
-cryst.pixelsize = 0.5*(parameters.labgeo.pixelsizeu + parameters.labgeo.pixelsizev);
+cryst.pixelsize = 0.5*(parameters.detgeo.pixelsizeu + parameters.detgeo.pixelsizev);
--
GitLab