-
Nicola Vigano authored
Signed-off-by:Nicola Vigano <nicola.vigano@esrf.fr>
Nicola Vigano authoredSigned-off-by:
gtCalculateGrain_p.m 15.31 KiB
function grain = gtCalculateGrain_p(grain, parameters, varargin)
% GTCALCULATEGRAIN_P Calulates diffraction vectors and detector intersection points
%
% grain = gtCalculateGrain_p(grain, parameters, varargin)
% -----------------------------------------------------
% Parallel version of gtCalculateGrain
%
% INPUT:
% grain = <struct>/<cell <struct> > grain(s) of interest
% parameters = <struct> parameters file
%
% OPTIONAL INPUT (as a list of pairs):
% showfigure = <logical> show or not the figure {false}
% color = <logical> show or not the prediction spot positions coloured as
% omega {false}
% axis = <logical> plot axis, image reference and beam direction in the
% image {false}
% overlay = <double> difspots image output from "gtPlaceDifspotinFull" (MxM)
% M is the size of the overlay image (typically 2048)
% markersize = <double> marker size {10}
% clims = <double> scale of grays {[-300 500]}
% usestrain = <logical> use strain info
%
% OUTPUT:
% grain = <struct>/<cell <struct> > grain(s) of interest
% Added fields:
% .pllab = <double>
% .hklsp = <double>
% .dvec = <double>
% .dvecsam = <double>
% .allblobs = <struct>
% Updated fields:
% .thetatype = <double>
% .hkl = <double>
% .pl = <double>
% .theta = <double>
% .eta = <double>
% .omega = <double>
%
% Simplified version of gtFedGenerateRandomGrain 15-03-2012 - WLudwig
if (nargin < 1)
disp('Usage: grain = gtCalculateGrain(grain, parameters, varargin)');
return
end
if (~exist('parameters','var') || isempty(parameters))
parameters = gtLoadParameters();
end
% set default values for optional arguments
app.color = false;
app.showfigure = false;
app.axis = false;
app.markersize = 10;
app.clims = [-300 500];
app.overlay = []; % zeros(parameters.acq.ydet, parameters.acq.xdet);
app.usestrain = false;
app.included = [];
app.ref_omind = [];
app.detectors = [];
app = parse_pv_pairs(app,varargin);
% If strain data is not available, disable strain
if (~isfield(grain, 'strain') || ~isfield(grain.strain, 'strainT') ...
|| isnan(grain.strain.strainT(1, 1)))
app.usestrain = false;
end
if (iscell(grain))
phase = unique(cellfun(@(x)x.phaseid, grain));
if (numel(phase) > 1)
error('gtCalculateGrain:wrong_argument', ...
'Cannot handle multiple phases!');
end
R_vec = cellfun(@(x){x.R_vector}, grain);
R_vec = cat(1, R_vec{:});
gcsam = cellfun(@(x){x.center}, grain);
gcsam = cat(3, gcsam{:});
num_grains = numel(grain);
else
% Rodrigues vector (grain orientation)
R_vec = grain.R_vector;
phase = grain.phaseid;
gcsam = grain.center;
num_grains = 1;
end
cryst = parameters.cryst(phase);
labgeo = parameters.labgeo;
samgeo = parameters.samgeo;
if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
labgeo.omstep = gtGetOmegaStepDeg(parameters);
end
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Detector geometry
if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
end
detgeo = parameters.detgeo;
if (~isempty(app.detectors))
detgeo = detgeo(app.detectors);
end
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Crystallography from parameters file
lambda = gtConvEnergyToWavelength(parameters.acq.energy);
%disp('Translating Miller indices into normalized cartesian coordinates for plane normals')
Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
% D-spacings in the reference crystal for the given reflections
dsp_ref = gtCrystDSpacing(double(cryst.hkl(:, cryst.thetatypesp)), Bmat);
% Plane normals in Cartesian CRYSTAL coordinates
plcry = gtCrystHKL2Cartesian(double(cryst.hklsp), Bmat);
% Compute all orientation matrices g (Reminder: Vc = g . Vs)
g = gtMathsRod2OriMat(R_vec');
% Express plane normals in cartesian SAMPLE CS
pl_sam = gtVectorCryst2Lab(plcry', g);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute diffraction angles and detector intersection points
rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
% Let's merge multple grains information for vectorization purposes
pl_samd = permute(pl_sam, [1 3 2]);
pl_samd = reshape(pl_samd, [], 3);
dsp_ref = dsp_ref(1, :, ones(1, num_grains));
dsp_ref = reshape(dsp_ref, 1, []);
num_plane_normals = size(pl_sam, 1);
% Impose strain on grain
if (app.usestrain)
% Deformation tensor (use strain tensor for now)
defT = grain.strain.strainT;
% New deformed plane normal and relative elongation (relative d-spacing)
[pl_samd, drel] = gtStrainPlaneNormals(pl_samd', defT); % unit vector
else
pl_samd = pl_samd';
drel = ones(1, num_plane_normals * num_grains);
end
% d-spacing after applying the strain tensor
dsp_d = drel .* dsp_ref;
% Bragg angles theta in the deformed state
sinth = 0.5 * lambda ./ dsp_d;
% New sin(theta) in deformed state is inversely proportional to elongation
theta = asind(sinth);
% The plane normals need to be brought in the Lab reference where the
% beam direction and rotation axis are defined.
% Use the Sample -> Lab orientation transformation assuming omega=0;
% (vector length preserved for free vectors)
pl_labd = gtGeoSam2Lab(pl_samd', eye(3), labgeo, samgeo, true, false)';
% Four omegas of the plane normal (1-3 and 2-4 are the two Friedel pairs):
% pls is plus or minus pl, and St is the rotation tensor per each omega
[om4, pllab4, pls4, rot_l2s_4, omind4] = gtFedPredictOmegaMultiple(pl_labd, ...
sinth, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
% Delete those where no reflection occurs
no_fwd_diffr = (isnan(om4(1, :)) | isnan(om4(3, :)));
% todel4 is 4 x (num_pn * num_gr) already!
no_fwd_diffr4 = no_fwd_diffr([1 1 1 1], :);
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Propagating the same common values for all the 4 omegas coming from a
% % single plane normal.
plorig4 = pl_samd(:, :, [1 1 1 1]);
if (size(cryst.hkl, 1) == 3 || size(cryst.hkl, 1) == 4)
hkl = cryst.hkl;
hklsp = cryst.hklsp;
else
hkl = cryst.hkl';
hklsp = cryst.hklsp';
end
num_hkl_indxs = size(hkl, 1);
thetatypesp = reshape(cryst.thetatypesp, [], 1);
thetatypesp = thetatypesp(:, ones(1, num_grains));
thetatypesp = reshape(thetatypesp, [], 1);
hkl = hkl(:, thetatypesp)';
hklsp = hklsp(:, :, ones(1, num_grains));
hklsp = reshape(hklsp, num_hkl_indxs, []);
hklsp = hklsp';
thetatypesp4 = thetatypesp(:, [1 1 1 1]);
hkl4 = hkl(:, :, [1 1 1 1]);
hklsp4 = hklsp(:, :, [1 1 1 1]);
hklsp4(:, :, [3 4]) = -hklsp4(:, :, [3 4]);
sinth4 = sinth([1 1 1 1], :);theta4 = theta([1 1 1 1], :);
if (isempty(app.ref_omind))
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Following the convention of the matching output, the omega value
% % smaller than 180deg (spot "a") is used in the pair data.
chind4 = om4(1, :) > om4(3, :);
pllab4(:, chind4, [1, 3]) = pllab4(:, chind4, [3, 1]);
pls4(:, chind4, [1, 3]) = pls4(:, chind4, [3, 1]);
om4([1, 3], chind4) = om4([3, 1], chind4);
omind4([1, 3], chind4) = omind4([3, 1], chind4);
rot_l2s_4(:, :, chind4, [1, 3]) = rot_l2s_4(:, :, chind4, [3, 1]);
hklsp4(chind4, :, [1, 3]) = hklsp4(chind4, :, [3, 1]);
chind4 = om4(2, :) > om4(4, :);
pllab4(:, chind4, [2, 4]) = pllab4(:, chind4, [4, 2]);
pls4(:, chind4, [2, 4]) = pls4(:, chind4, [4, 2]);
om4([2, 4], chind4) = om4([4, 2], chind4);
omind4([2, 4], chind4) = omind4([4, 2], chind4);
rot_l2s_4(:, :, chind4, [2, 4]) = rot_l2s_4(:, :, chind4, [4, 2]);
hklsp4(chind4, :, [2, 4]) = hklsp4(chind4, :, [4, 2]);
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Four omegas of the plane normal (1-3 and 2-4 are the two Friedel pairs):
% % But the two Friedel pairs are going to be 1a-1b and 2a-2b.
pllab4(:, :, [2, 3]) = pllab4(:, :, [3, 2]);
pls4(:, :, [2, 3]) = pls4(:, :, [3, 2]);
om4([2, 3], :) = om4([3, 2], :);
omind4([2, 3], :) = omind4([3, 2], :);
rot_l2s_4(:, :, :, [2, 3]) = rot_l2s_4(:, :, :, [3, 2]);
hklsp4(:, :, [2, 3]) = hklsp4(:, :, [3, 2]);
end
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Unfolding the 4-fold data structures
no_fwd_diffr = reshape(no_fwd_diffr4, [], num_grains);
om = reshape(om4, [], num_grains);
omind = reshape(omind4, [], num_grains);
sinth = reshape(sinth4, [], num_grains);
theta = reshape(theta4, [], num_grains);
% pllab4, pls4, and plorig4 are: 3 x (num_pn * num_gr) x 4
pllab = reshape(permute(pllab4, [1 3 2]), 3, [], num_grains);
pls = reshape(permute(pls4, [1 3 2]), 3, [], num_grains);
plorig = reshape(permute(plorig4, [1 3 2]), 3, [], num_grains);
rot_l2s = reshape(permute(rot_l2s_4, [1 2 4 3]), 3, 3, [], num_grains);
thetatypesp = reshape(thetatypesp4', [], num_grains);
hkl = permute(hkl4, [2 3 1]);
hkl = reshape(hkl, num_hkl_indxs, [], num_grains);
hklsp = permute(hklsp4, [2 3 1]);
hklsp = reshape(hklsp, num_hkl_indxs, [], num_grains);
if (~isempty(app.ref_omind))
grs_omind_ind = [ ...
find(omind == 1)'; find(omind == 2)'; ...
find(omind == 3)'; find(omind == 4)' ];
grs_omind_ind = reshape(grs_omind_ind, [], 1);
ref_omind_ind = [ ...
find(app.ref_omind == 1)'; find(app.ref_omind == 2)'; ...
find(app.ref_omind == 3)'; find(app.ref_omind == 4)' ];
ref_omind_ind = reshape(ref_omind_ind, [], 1);
ref_omind_ind = repmat(ref_omind_ind, [1 num_grains]) ...
+ repmat((0:num_grains-1) * num_plane_normals*4, [num_plane_normals*4, 1]);
ref_omind_ind = reshape(ref_omind_ind, [], 1);
no_fwd_diffr(ref_omind_ind) = no_fwd_diffr(grs_omind_ind);
om(ref_omind_ind) = om(grs_omind_ind);
omind(ref_omind_ind) = omind(grs_omind_ind);
sinth(ref_omind_ind) = sinth(grs_omind_ind);
theta(ref_omind_ind) = theta(grs_omind_ind);
pllab(:, ref_omind_ind) = pllab(:, grs_omind_ind);
pls(:, ref_omind_ind) = pls(:, grs_omind_ind);
plorig(:, ref_omind_ind) = plorig(:, grs_omind_ind);
rot_l2s(:, :, ref_omind_ind) = rot_l2s(:, :, grs_omind_ind);
thetatypesp(ref_omind_ind) = thetatypesp(grs_omind_ind);
hkl(:, ref_omind_ind) = hkl(:, grs_omind_ind);
hklsp(:, ref_omind_ind) = hklsp(:, grs_omind_ind);
end
if (~isempty(app.included))
no_fwd_diffr = no_fwd_diffr(app.included, :);
om = om(app.included, :);
omind = omind(app.included, :);
sinth = sinth(app.included, :);
theta = theta(app.included, :);
pllab = pllab(:, app.included, :);
pls = pls(:, app.included, :);
plorig = plorig(:, app.included, :);
rot_l2s = rot_l2s(:, :, app.included, :);
thetatypesp = thetatypesp(app.included, :);
hkl = hkl(:, app.included, :);
hklsp = hklsp(:, app.included, :);
end
no_fwd_diffr = reshape(no_fwd_diffr, [], 1);
om = reshape(om, [], 1);
pllab = reshape(pllab, 3, []);
plorig = permute(plorig, [2 1 3]);
pls = permute(pls, [2 1 3]);
rot_l2s = reshape(rot_l2s, 3, 3, []);
hkl = permute(hkl, [2 1 3]);
hklsp = permute(hklsp, [2 1 3]);
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Computing the remaining physical quantities
% Diffraction vectors % takes coloumn vectors
dveclab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
% Normalize vectors
dvec_norms = sqrt(sum(dveclab .* dveclab, 1));
dveclab = (dveclab ./ dvec_norms([1 1 1], :))';
% Eta angles % takes row vector
eta = gtGeoEtaFromDiffVec(dveclab, labgeo);
% Diffraction vector in Sample referencedvecsam = gtGeoLab2Sam(dveclab, om, labgeo, samgeo, true);
% Let's convert grain center Sam -> Lab
rot_s2l = permute(rot_l2s, [2 1 3]);
gcsam_v = permute(gcsam, [1 3 2]);
gcsam_v = gcsam_v(ones(1, size(rot_l2s, 3)/num_grains), :, :);
gcsam_v = reshape(gcsam_v, [], 3);
gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
% Lorentz factor for the given projections (note that eq. 4.1 from page 36
% in Henning's book is not correct)
Lfactor = 1 ./ abs(sind(eta));
num_dets = numel(detgeo);
ondet = cell(num_dets, 1);
uvw = cell(num_dets, 1);
for n = 1:num_dets
% u,v,w coordinates on the detector
uvw{n} = gtFedPredictUVWMultiple([], dveclab', gclab_v', ...
detgeo(n).detrefpos', detgeo(n).detnorm', detgeo(n).Qdet, ...
[detgeo(n).detrefu, detgeo(n).detrefv]', om', labgeo.omstep);
% Precompute which reflections will actually fall into the detector
ondet{n} = ~no_fwd_diffr' & ...
(uvw{n}(1, :) >= 1) & (detgeo(n).detsizeu >= uvw{n}(1, :)) & ...
(uvw{n}(2, :) >= 1) & (detgeo(n).detsizev >= uvw{n}(2, :));
end
if (num_grains > 1)
for n = 1:num_dets
uvw{n} = reshape(uvw{n}, 3, [], num_grains);
uvw{n} = permute(uvw{n}, [2 1 3]);
ondet{n} = reshape(ondet{n}, [], num_grains);
end
pllab = reshape(pllab, 3, [], num_grains);
pllab = permute(pllab, [2 1 3]);
eta = reshape(eta, [], num_grains);
Lfactor = reshape(Lfactor, [], num_grains);
om = reshape(om, [], num_grains);
dveclab = reshape(dveclab', 3, [], num_grains);
dvecsam = reshape(dvecsam', 3, [], num_grains);
dveclab = permute(dveclab, [2 1 3]);
dvecsam = permute(dvecsam, [2 1 3]);
rot_l2s = reshape(rot_l2s, 3, 3, [], num_grains);
for ii_gr = 1:num_grains
% Initialse output variables
for n = 1:num_dets
grain{ii_gr}.allblobs.detector(n) = struct( ...
'uvw', {uvw{n}(:, :, ii_gr)}, ...
'ondet', {ondet{n}(:, ii_gr)});
end
grain{ii_gr}.allblobs.plorig = plorig(:, :, ii_gr);
grain{ii_gr}.allblobs.pl = pls(:, :, ii_gr);
grain{ii_gr}.allblobs.pllab = pllab(:, :, ii_gr);
grain{ii_gr}.allblobs.hkl = hkl(:, :, ii_gr);
grain{ii_gr}.allblobs.hklsp = hklsp(:, :, ii_gr);
grain{ii_gr}.allblobs.sintheta = sinth(:, ii_gr);
grain{ii_gr}.allblobs.theta = theta(:, ii_gr);
grain{ii_gr}.allblobs.thetatype = thetatypesp(:, ii_gr); grain{ii_gr}.allblobs.eta = eta(:, ii_gr);
grain{ii_gr}.allblobs.lorentzfac = Lfactor(:, ii_gr);
grain{ii_gr}.allblobs.omega = om(:, ii_gr);
grain{ii_gr}.allblobs.omind = omind(:, ii_gr);
grain{ii_gr}.allblobs.dvec = dveclab(:, :, ii_gr);
grain{ii_gr}.allblobs.dvecsam = dvecsam(:, :, ii_gr);
grain{ii_gr}.allblobs.srot = rot_l2s(:, :, :, ii_gr);
grain{ii_gr}.full = [];
end
else
% Initialse output variables
for n = 1:num_dets
grain.allblobs.detector(n) = struct('uvw', uvw(n), 'ondet', ondet(n));
end
grain.allblobs.plorig = plorig;
grain.allblobs.pl = pls;
grain.allblobs.pllab = pllab';
grain.allblobs.hkl = hkl;
grain.allblobs.hklsp = hklsp;
grain.allblobs.sintheta = sinth;
grain.allblobs.theta = theta;
grain.allblobs.eta = eta;
grain.allblobs.thetatype = thetatypesp;
grain.allblobs.lorentzfac = Lfactor;
grain.allblobs.omega = om;
grain.allblobs.omind = omind;
grain.allblobs.dvec = dveclab;
grain.allblobs.dvecsam = dvecsam;
grain.allblobs.srot = rot_l2s;
grain.full = [];
end
end