Skip to content
Snippets Groups Projects
Commit 6d96227a authored by Nicola Vigano's avatar Nicola Vigano
Browse files

detgeo: another round of fixes


(thanks to Lorenzo for spotting some of them once again)

Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 57641df3
No related branches found
No related tags found
No related merge requests found
......@@ -278,7 +278,7 @@ function parameters = sfGetReflectionInfoFromFable(parameters, phase_id)
disp('Using fable calculation of theoretical reflections')
disp('The only missing information is the intensity')
% reflection list
list = gtCrystCalculateReflections(parameters.cryst(phase_id), parameters.labgeo, parameters.acq.energy, 'reflections_');
list = gtCrystCalculateReflections(parameters.cryst(phase_id), parameters.detgeo, parameters.acq.energy, 'reflections_');
parameters.cryst(phase_id) = gtAddMatFile(parameters.cryst(phase_id), list, true, false, false, false, false);
clear list
......
......@@ -200,7 +200,7 @@ else
if strcmpi(check,'y')
parameters.acq.rotu = gtFindRotationAxis(parameters.acq,bbdir);
parameters.acq.rotx = parameters.acq.rotu;
parameters.labgeo.detrefpos(2) = -(parameters.acq.rotu - parameters.acq.xdet / 2) * parameters.acq.pixelsize;
parameters.detgeo.detrefpos(2) = -(parameters.acq.rotu - parameters.acq.xdet / 2) * parameters.detgeo.pixelsizeu;
end
save(parameters_name,'parameters');
......
......@@ -42,9 +42,9 @@ totproj = totproj-1;
% One quarter of detector, 10 images
set(handles.u_start, 'string', '1')
set(handles.u_end, 'string', num2str(round(parameters.labgeo.detsizeu/2)))
set(handles.u_end, 'string', num2str(round(parameters.detgeo.detsizeu/2)))
set(handles.v_start, 'string', '1')
set(handles.v_end, 'string', num2str(round(parameters.labgeo.detsizev/2)))
set(handles.v_end, 'string', num2str(round(parameters.detgeo.detsizev/2)))
set(handles.image_start, 'string', num2str(round(totproj/2)))
set(handles.image_end, 'string', num2str(round(totproj/2)+9))
% start with the first image
......
function list = gtCrystCalculateReflections(cryst, labgeo, energy, filename)
function list = gtCrystCalculateReflections(cryst, detgeo, energy, filename)
% GTCRYSTCALCULATEREFLECTIONS Calculate the reflections list using python
% from xfab library in fable
% list = gtCrystCalculateReflections(cryst, labgeo, energy, [filename])
% list = gtCrystCalculateReflections(cryst, detgeo, energy, [filename])
% -------------------------------------------------------------------
% List of theoretical reflections for a given crystal system and spacegroup
%
......@@ -30,9 +30,9 @@ function list = gtCrystCalculateReflections(cryst, labgeo, energy, filename)
latticepar = cryst.latticepar;
sg = cryst.spacegroup;
tmp = labgeo.detanglemin/2; % theta
tmp = detgeo.detanglemin/2; % theta
minsintl = min(0.0001, sind(tmp)/gtConvEnergyToWavelength(energy));
tmp = labgeo.detanglemax/2;
tmp = detgeo.detanglemax/2;
maxsintl = sind(tmp)/gtConvEnergyToWavelength(energy);
if ~exist('filename','var') || isempty(filename)
......
function [diffvec, labAXYZ, theta] = gtGeoDiffVecInLab(detAUV, labBXYZ, labgeo)
% FUNCTION [diffvec, labAXYZ, theta] = gtGeoDiffVecInLab(detAUV, labBXYZ, labgeo)
%
% FUNCTION [diffvec,labAXYZ,theta] = gtGeoDiffVecInLab(detAUV,labBXYZ,labgeo)
%
% Given a set of points A on the detector with their Detector coordinates
% U,V, and another set of points B with their Lab coordinates X,Y,Z,
% Given a set of points A on the detector with their Detector coordinates
% U, V, and another set of points B with their Lab coordinates X, Y, Z,
% the function returns the corresponding diffraction vectors in the Lab
% reference. The diffraction vectors indicate the direction of the
% reference. The diffraction vectors indicate the direction of the
% diffracted beam paths.
% It also returns the Bragg angles and the Lab coordinates of points A.
% It also returns the Bragg angles and the Lab coordinates of points A.
%
% INPUT detAUV - Detector coordinates U,V of the points A (n,2)
% labBXYZ - Lab coordinates X,Y,Z of points B (n,3)
% INPUT detAUV - Detector coordinates U, V of the points A (n, 2)
% labBXYZ - Lab coordinates X, Y, Z of points B (n, 3)
%
% Geometry parameters in Lab coordinates:
%
% labgeo.rotpos - rotation axis position (row vector)
% labgeo.rotdir - rotation axis direction (unit row vector)
% labgeo.detorig - position of detector origin (row vector)
% labgeo.detdiru - detector U direction (unit row vector)
% 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)
% detgeo.detorig - position of detector origin (row vector)
% detgeo.detdiru - detector U direction (unit row vector)
% detgeo.detdirv - detector V direction (unit row vector)
% detgeo.pixelsizeu - pixel size U (same unit as position values)
% detgeo.pixelsizev - pixel size V (same unit as position values)
%
% OUTPUT diffvec - Lab coordinates of the normalised diffraction vectors;
% direction points from points B to A (n,3)
% labAXYZ - Lab coordinates of points A (n,3)
% theta - Bragg angles theta in degrees (n,1)
% direction points from points B to A (n, 3)
% labAXYZ - Lab coordinates of points A (n, 3)
% theta - Bragg angles theta in degrees (n, 1)
%
function [diffvec,labAXYZ,theta] = gtGeoDiffVecInLab(detAUV,labBXYZ,labgeo)
% Coordinate transformation from Detector to Lab
labAXYZ = gtGeoDet2Lab(detAUV,labgeo,0);
% Diffraction vectors from labXYZ to point A:
diffvec = labAXYZ - labBXYZ;
% Coordinate transformation from Detector to Lab
labAXYZ = gtGeoDet2Lab(detAUV, labgeo, 0);
% Normalise diffvec-s (seems the fastest way for nx3 vectors):
dnorm = sqrt(diffvec(:,1).^2 + diffvec(:,2).^2 + diffvec(:,3).^2) ;
% Diffraction vectors from labXYZ to point A:
diffvec = labAXYZ - labBXYZ;
diffvec(:,1) = diffvec(:,1)./dnorm ;
diffvec(:,2) = diffvec(:,2)./dnorm ;
diffvec(:,3) = diffvec(:,3)./dnorm ;
% Normalise diffvec-s (seems the fastest way for nx3 vectors):
dnorm = sqrt(diffvec(:, 1) .^ 2 + diffvec(:, 2) .^ 2 + diffvec(:, 3) .^2 );
diffvec(:, 1) = diffvec(:, 1) ./ dnorm;
diffvec(:, 2) = diffvec(:, 2) ./ dnorm;
diffvec(:, 3) = diffvec(:, 3) ./ dnorm;
if nargout == 3
% Use the dot products of beamdir and diffvec to get theta
theta = 0.5*acosd(diffvec*labgeo.beamdir');
if (nargout == 3)
% Use the dot products of beamdir and diffvec to get theta
theta = 0.5 * acosd(diffvec * labgeo.beamdir');
end
end
function [detUV,lab2det] = gtGeoLab2Det(labXYZ,labgeo,freevec_flag)
%
% FUNCTION [detUV,lab2det] = gtGeoLab2Det(labXYZ,labgeo,freevec_flag)
function [detUV, lab2det] = gtGeoLab2Det(labXYZ, detgeo, freevec_flag)
% FUNCTION [detUV, lab2det] = gtGeoLab2Det(labXYZ, detgeo, freevec_flag)
%
% LAB -> DETECTOR coordinate transformation using arbitrary geometry.
%
......@@ -8,51 +7,46 @@ function [detUV,lab2det] = gtGeoLab2Det(labXYZ,labgeo,freevec_flag)
% Detector coordinates. If labXYZ is not in the detector plane,
% its perpendicular projection onto the detector plane will be used.
%
% INPUT labXYZ - array of LAB coordinates (n,3)
% INPUT labXYZ - array of LAB coordinates (n, 3)
%
% Geometry parameters given in LAB coordinates:
% labgeo.detorig - position of detector origin: pixel {0,0};
% size (1,3)
% labgeo.detdiru - detector U direction (unit row vector)
% labgeo.detdirv - detector V direction (unit row vector)
% labgeo.pixelsizeu - pixel size U (same unit as position)
% labgeo.pixelsizev - pixel size V (same unit as position)
% detgeo.detorig - position of detector origin: pixel {0, 0};
% size (1, 3)
% detgeo.detdiru - detector U direction (unit row vector)
% detgeo.detdirv - detector V direction (unit row vector)
% detgeo.pixelsizeu - pixel size U (same unit as position)
% detgeo.pixelsizev - pixel size V (same unit as position)
%
% freevec_flag - 1 for free vectors, 0 for bound vectors
%
% OUTPUT detUV - detector coordinates in pixels (n,2)
% lab2det - transformation matrix from Lab to Detector (3,2)
% OUTPUT detUV - detector coordinates in pixels (n, 2)
% lab2det - transformation matrix from Lab to Detector (3, 2)
%
if isfield(detgeo, 'lab2det')
lab2det = detgeo.lab2det;
else
% Determine the det2lab transformation matrix from geometry.
% The u, v basis may not be orthonormal, not like the Lab reference which
% is an orthonormal basis by definition. Therefore here the decomposition
% of labXYZ into u and v vectors needs to be done adequately.
% The transformation matrix lab2det is the same as the pseudo-inverse
% of det2lab.
u = detgeo.detdiru;
v = detgeo.detdirv;
uv = u*v';
if isfield(labgeo,'lab2det')
lab2det = labgeo.lab2det;
else
% Determine the det2lab transformation matrix from geometry.
% The u,v basis may not be orthonormal, not like the Lab reference which
% is an orthonormal basis by definition. Therefore here the decomposition
% of labXYZ into u and v vectors needs to be done adequately.
% The transformation matrix lab2det is the same as the pseudo-inverse
% of det2lab.
u = labgeo.detdiru;
v = labgeo.detdirv;
uv = u*v';
a = (u-v*uv)/(1-uv*uv);
b = v-a*uv;
lab2det = [a/labgeo.pixelsizeu; b/labgeo.pixelsizev]';
end
% New coordinates
a = (u - v * uv) / (1 - uv * uv);
b = v - a * uv;
if freevec_flag % free vectors
detUV = labXYZ*lab2det;
lab2det = [a / detgeo.pixelsizeu; b / detgeo.pixelsizev]';
end
else % offset vectors
% New coordinates
detUV = (labXYZ - labgeo.detorig(ones(size(labXYZ,1),1),:)) * lab2det;
if (freevec_flag) % free vectors
detUV = labXYZ * lab2det;
else % offset vectors
detUV = (labXYZ - detgeo.detorig(ones(size(labXYZ, 1), 1), :)) * lab2det;
end
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment