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
%
%
%% 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';
detrefu = labgeo.detrefu;
detrefv = labgeo.detrefv;
Qdet = gtFedDetectorProjectionTensor(detdiru,detdirv,1,1);
tn = cross(detdiru,detdirv);
detpos = (labgeo.detrefpos./acq.pixelsize)';
detposcorner = detpos - detrefu*detdiru - detrefv*detdirv; % the (0,0) corner of the detector in Lab system (edge of detector pixel)
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(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 = [];
grain.allblobs.proj_geom = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute diffraction angles and detector intersection points
rotcomp = gtFedRotationMatrixComp(rotdir);
rotcomp_maths = gtMathsRotationMatrixComp(rotdir','row');
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);
S = gtMathsRotationTensor(om, rotcomp_maths);
% 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;
s1a = S(:,:,1)';
s1b = S(:,:,3)';
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;
s1a = S(:,:,3)';
s1b = S(:,:,1)';
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;
s2a = S(:,:,2)';
s2b = S(:,:,4)';
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;
s2a = S(:,:,4)';
s2b = S(:,:,2)';
pl2a = pls(:,4)';
pl2b = pls(:,2)';
pllab2a = pllab(:,4)';
pllab2b = pllab(:,2)';
d2a = d4';
d2b = d2';
end
% Eta angles
eta1a = gtGeoEtaFromDiffVec(d1a, labgeo); % takes row vector
eta1b = gtGeoEtaFromDiffVec(d1b, labgeo);
eta2a = gtGeoEtaFromDiffVec(d2a, labgeo);
eta2b = gtGeoEtaFromDiffVec(d2b, labgeo);
% 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
% add projection geometry for grain center of mass (using Rottens'
% instead Rottens in order to perfrom negative rotation)
r_th1a = s1a'*beamdir + 2*sinth*pl1a';
r_th1b = s1b'*beamdir + 2*sinth*pl1b';
r_th2a = s2a'*beamdir + 2*sinth*pl2a';
r_th2b = s2b'*beamdir + 2*sinth*pl2b';
detpos1a = detposcorner + uv1a(1)*detdiru + uv1a(2)*detdirv;
detpos1b = detposcorner + uv1b(1)*detdiru + uv1b(2)*detdirv;
detpos2a = detposcorner + uv2a(1)*detdiru + uv2a(2)*detdirv;
detpos2b = detposcorner + uv2b(1)*detdiru + uv2b(2)*detdirv;
proj_geom1a = [r_th1a; s1a'*detpos1a - csam; s1a'*detdiru; s1a'*detdirv]';
proj_geom1b = [r_th1b; s1b'*detpos1b - csam; s1b'*detdiru; s1b'*detdirv]';
proj_geom2a = [r_th2a; s2a'*detpos2a - csam; s2a'*detdiru; s2a'*detdirv]';
proj_geom2b = [r_th2b; s2b'*detpos2b - csam; s2b'*detdiru; s2b'*detdirv]';
% 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.proj_geom = [grain.allblobs.proj_geom; proj_geom1a; proj_geom1b; proj_geom2a; proj_geom2b];
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