function grain = gtCalculateGrain(grain, parameters, varargin)
% GTCALCULATEGRAIN Calulate diffraction vectors and detector intersection points
%
% grain = gtCalculateGrain(grain, parameters, varargin)
% ----------------------------------------------------------
% Input and output arguments description, variable type (in <>) and default value (in {}) are written below.
%
%
% INPUT:
%
% grain = grain of interest (<struct>)
% parameters = parameters file (<struct> {parameters.mat})
% varargin = optional arguments list given by pairs
%
%
% OPTIONAL INPUT (as a list of pairs):
%
% showfigure = show or not the figure (<logical> {0})
% color = show or not the prediction spot positions coloured as
% omega (<logical> {0})
% axis = plot axis, image reference and beam direction in the
% image (<logical> {0})
% overlay = difspots image output from "gtPlaceDifspotinFull" (<double MxM>)
% M is the size of the overlay image (typically 2048)
% Markersize = marker size (<double> {15})
% Grayscale = scale of grays (<double 1x2> {[-300 500]})
%
%
% OUTPUT:
%
% grain = grain of interest (<struct>)
% added fields:
% - pllab <double Nx3>
% - hklsp <double Nx3>
% - dvec <double Nx3>
% - allblobs <struct>
%
% updated fields:
% - thetatype <double Nx1>
% - hkl <double Nx3>
% - pl <double Nx3>
% - theta <double Nx1>
% - eta <double Nx1>
% - omega <double Nx1>
%
% Simplified version of gtFedGenerateRandomGrain 15-03-2012 - WLudwig
%
%
% FUNCTIONS USED:
%[use]- gtConvEnergyToWavelength
%[use]- gtCrystSpacegroupReflections
%[use]- gtCrystSignedHKLs
%[use]- gtHKL2Cart
%[use]- gtFedDetectorProjectionTensor
%[use]- gtFedRotationMatrixComp
%[use]- gtFedPredictOmega
%[use]- gtFedRotationTensor
%[use]- gtFedPredictDiffVec
%[use]- gtMatchEtaOfPoint
%[use]- gtFedPredictUVW
%% set default values for optional arguments
app.color = 0;
app.showfigure = 0;
app.axis = 0;
app.Markersize = 15;
app.Grayscale = [-300 500];
app.overlay = zeros(parameters.acq.ydet, parameters.acq.xdet);
app = parse_pv_pairs(app,varargin);
disp('Initial settings')
disp(app)
if nargin < 1
disp('Usage: grain = gtCalculateGrain(grain, parameters, varargin)')
disp(' grain = grain{grainid}')
disp(' parameters = parameters file')
return
end
if ~exist('parameters','var') || isempty(parameters)
parameters=load('parameters.mat');
end
% Rodrigues vector (grain orientation)
Rvec = grain.R_vector;
disp(['grain.center: ' num2str(grain.center)]);
phase = grain.phaseID;
acq = parameters.acq;
cryst = parameters.cryst(phase);
labgeo = parameters.labgeo;
samgeo = parameters.samgeo;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detector geometry
beamdir = labgeo.beamdir';
rotdir = labgeo.rotdir';
detdiru = labgeo.detdiru';
detdirv = labgeo.detdirv';
Qdet = gtFedDetectorProjectionTensor(detdiru,detdirv,1,1);
tn = cross(detdiru,detdirv);
detpos = (labgeo.detrefpos./acq.pixelsize)';
uvorig = [labgeo.detrefu, labgeo.detrefv]';
csam = grain.center' / acq.pixelsize;
omstep = 180/acq.nproj;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Crystallography from parameters file
lambda = gtConvEnergyToWavelength(acq.energy);
hklsp = cryst.hklsp;
d0 = cryst.dspacingsp;
thetasp = cryst.thetasp;
hklt = cryst.hkl;
thetatypesp = cryst.thetatypesp;
disp('Translating Miller indices into normalized cartesian coordinates for plane normals')
Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
for ii = 1:size(hklsp,2)
%pl0(i, :) = gtHKL2Cart(hklsp(:, i)', cryst);
% !!! should be changed to:
pl0(ii, :) = gtCrystHKL2Cartesian(hklsp(:, ii), Bmat)';
end
gg0 = Rod2g(Rvec); % gg0 is the transformation tensor from SAMPLE TO CRYSTAL coordinates (Rotation)
pl_g=(gg0*pl0')'; % Plane normals in SAMPLE coord. rotated with gg0:
%% Initialse output variables
grain.allblobs.pl = [];
grain.allblobs.pllab = [];
grain.allblobs.hkl = [];
grain.allblobs.hklsp = [];
grain.allblobs.thetatype = [];
grain.allblobs.theta = [];
grain.allblobs.eta = [];
grain.allblobs.omega = [];
grain.allblobs.dvec = [];
grain.allblobs.uv = [];
grain.allblobs.uvw = [];
grain.allblobs.srot = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute diffraction angles and detector intersection points
rotcomp = gtFedRotationMatrixComp(rotdir);
for ii=1:size(pl_g,1)
%hkl(i,:) = hklt(thetatypesp(i),:);
hkl(:,ii) = hklt(:,thetatypesp(ii))';
pl = pl_g(ii,:)';
sinth = sind(thetasp(ii));
% Four omegas of the plane normal (1-3 and 2-4 are the two Friedel pairs):
% pls is plus or minus pl
[om,pllab,pls] = gtFedPredictOmega(pl, sinth, beamdir, rotdir, rotcomp);
% If reflection occurs
if ~isempty(om)
% Rotation tensors from omegas (one for each of the two Friedel pairs)
S1 = gtFedRotationTensor(om(1), rotcomp);
S2 = gtFedRotationTensor(om(2), rotcomp);
S3 = gtFedRotationTensor(om(3), rotcomp);
S4 = gtFedRotationTensor(om(4), rotcomp);
% Diffraction vectors
d1 = gtFedPredictDiffVec(pllab(:,1), sinth, beamdir); % takes coloumn vectors
d2 = gtFedPredictDiffVec(pllab(:,2), sinth, beamdir);
d3 = gtFedPredictDiffVec(pllab(:,3), sinth, beamdir);
d4 = gtFedPredictDiffVec(pllab(:,4), sinth, beamdir);
% Normalize vectors
d1 = d1/norm(d1);
d2 = d2/norm(d2);
d3 = d3/norm(d3);
d4 = d4/norm(d4);
% Following the convention of the mathcing output, the omega value
% smaller than 180deg (spot "a") is used in the pair data.
% The two Friedel pairs are 1a-1b and 2a-2b.
if om(1) < om(3)
om1a = om(1);
s1a = S1;
s1b = S3;
om1b = om(3);
pl1a = pls(:,1)';
pl1b = pls(:,3)';
pllab1a = pllab(:,1)';
pllab1b = pllab(:,3)';
d1a = d1';
d1b = d3';
else
om1a = om(3);
s1a = S3;
s1b = S1;
om1b = om(1);
pl1a = pls(:,3)';
pl1b = pls(:,1)';
pllab1a = pllab(:,3)';
pllab1b = pllab(:,1)';
d1a = d3';
d1b = d1';
end
if om(2) < om(4)
om2a = om(2);
om2b = om(4);
s2a = S2;
s2b = S4;
pl2a = pls(:,2)';
pl2b = pls(:,4)';
pllab2a = pllab(:,2)';
pllab2b = pllab(:,4)';
d2a = d2';
d2b = d4';
else
om2a = om(4);
om2b = om(2);
s2a = S4;
s2b = S2;
pl2a = pls(:,4)';
pl2b = pls(:,2)';
pllab2a = pllab(:,4)';
pllab2b = pllab(:,2)';
d2a = d4';
d2b = d2';
end
% Eta angles
eta1a = gtMatchEtaOfPoint(d1a); % takes row vector
eta1b = gtMatchEtaOfPoint(d1b);
eta2a = gtMatchEtaOfPoint(d2a);
eta2b = gtMatchEtaOfPoint(d2b);
% u,v coordinates on the detector
uv1a=gtFedPredictUVW(s1a,d1a',csam,detpos,tn,Qdet,uvorig,om1a,omstep);
uv1b=gtFedPredictUVW(s1b,d1b',csam,detpos,tn,Qdet,uvorig,om1b,omstep);
uv2a=gtFedPredictUVW(s2a,d2a',csam,detpos,tn,Qdet,uvorig,om2a,omstep);
uv2b=gtFedPredictUVW(s2b,d2b',csam,detpos,tn,Qdet,uvorig,om2b,omstep);
if isempty(uv1a)
uv1a = NaN(3,1);
end
if isempty(uv1b)
uv1b = NaN(3,1);
end
if isempty(uv2a)
uv2a = NaN(3,1);
end
if isempty(uv2b)
uv2b = NaN(3,1);
end
% fill the output structure with the two pairs
grain.allblobs.pl = [grain.allblobs.pl; pl1a; pl1b; pl2a; pl2b];
grain.allblobs.pllab = [grain.allblobs.pllab; pllab1a; pllab1b; pllab2a; pllab2b];
grain.allblobs.hkl = [grain.allblobs.hkl; hkl(:,ii)'; hkl(:,ii)'; hkl(:,ii)'; hkl(:,ii)'];
grain.allblobs.hklsp = [grain.allblobs.hklsp; hklsp(:,ii)'; -hklsp(:,ii)'; hklsp(:,ii)'; -hklsp(:,ii)'];
grain.allblobs.theta = [grain.allblobs.theta; thetasp(ii); thetasp(ii); thetasp(ii); thetasp(ii)];
grain.allblobs.thetatype = [grain.allblobs.thetatype; thetatypesp(ii); thetatypesp(ii); thetatypesp(ii); thetatypesp(ii)];
grain.allblobs.eta = [grain.allblobs.eta; eta1a; eta1b; eta2a; eta2b];
grain.allblobs.omega = [grain.allblobs.omega; om1a; om1b; om2a; om2b];
grain.allblobs.dvec = [grain.allblobs.dvec; d1a; d1b; d2a; d2b];
grain.allblobs.uv = [grain.allblobs.uv;uv1a(1),uv1a(2);uv1b(1),uv1b(2);uv2a(1),uv2a(2);uv2b(1),uv2b(2)];
grain.allblobs.uvw = [grain.allblobs.uvw;uv1a(1),uv1a(2),uv1a(3);uv1b(1),uv1b(2),uv1b(3);uv2a(1),uv2a(2),uv2a(3);uv2b(1),uv2b(2),uv2b(3)];
grain.allblobs.srot = [grain.allblobs.srot; s1a; s1b; s2a; s2b];
end % end if ~empty(om)
end % end for pl_g
if app.showfigure
figure
imshow(app.overlay, app.Grayscale);
hold on
if app.axis
detusize=acq.xdet;
detvsize=acq.ydet;
% Image frame
plot([0 detusize+1],[0 0],'k')
plot([0 detusize+1],[detvsize+1 detvsize+1],'k')
plot([0 0],[0 detvsize+1],'k')
plot([detusize+1 detusize+1],[0 detvsize+1],'k')
% Midlines
plot([0 detusize+1],[detvsize/2+0.5 detvsize/2+0.5],'-.k')
plot([detusize/2+0.5 detusize/2+0.5],[0 detvsize+1],'-.k')
% Set figure propeties
axis equal
set(gca,'YDir','reverse')
set(gca,'Position',[0.05 0.05 0.95 0.9])
xlim([-detusize*0.1 detusize*1.1])
ylim([-detvsize*0.1 detvsize*1.1])
xlabel('U direction')
ylabel('V direction')
% Beam direction in image
beamuv = Qdet*beamdir;
% Arrow in figure indicating beam direction
quiver(uvorig(1),uvorig(2),beamuv(1),beamuv(2),1000,'-k','Linewidth',3)
% Rotation axis direction in image
rotuv = Qdet*rotdir;
% Arrow in figure indicating rotation axis direction
quiver(uvorig(1),uvorig(2),rotuv(1),rotuv(2),1000,'-.k','Linewidth',3)
end
impixelinfo
om = grain.allblobs.omega;
uv = grain.allblobs.uv;
cmap = jet(1001);
omtoRGB = @(om) cmap(round((om)/(360)*1000+1),:);
if app.color
colormap(cmap)
colorbar
caxis([0 360])
else
colormap(gray);
end
for ii = 1 : length(uv)
if ~isempty(uv(ii,1)) % there might be no intersection when using offaxis detector configurations (vertical)
if app.color
plot(uv(ii,1),uv(ii,2),'or','MarkerSize',app.Markersize,'MarkerEdgeColor',omtoRGB(om(ii)));
else
plot(uv(ii,1),uv(ii,2),'ob','MarkerSize',app.Markersize);
end
end
end
end % showfigure
end