gtCalculateGrain.m

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
%

if ~exist('parameters','var') || isempty(parameters)
    parameters = load('parameters.mat');
    parameters = parameters.parameters;
end

% set default values for optional arguments
app.color = false;
app.showfigure = false;
app.axis = false;
app.Markersize = 15;
app.Grayscale = [-300 500];
app.overlay = zeros(parameters.acq.ydet, parameters.acq.xdet);
app.usestrain = false;
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 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

% 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);

pixelsize = mean([labgeo.pixelsizeu labgeo.pixelsizev]);

detpos  = (labgeo.detrefpos./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' ./ 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);

% Plane normals in Cartesian CRYSTAL coordinates
for ii = 1:size(hklsp,2)
    pl_cry(ii, :) = gtCrystHKL2Cartesian(hklsp(:, ii), Bmat)';
end

% Compute all orientation matrices g (Reminder: Vc = g . Vs)
g = gtMathsRod2RotMat(Rvec);

% Express plane normals in cartesian SAMPLE CS
pl_sam = gtVectorCryst2Lab(pl_cry, g);

% 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_sam, 1)

    %hkl(i,:) = hklt(thetatypesp(i),:);
    hkl(:,ii) = hklt(:,thetatypesp(ii))';

    % 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, drel] = gtStrainPlaneNormals(pl_sam(ii,:)', defT); % unit vector

        % New sin(theta) in deformed state is inversely proportional to elongation
        sinth = sind(thetasp(ii))/drel;
        theta = asind(sinth);
    else
        pl    = pl_sam(ii,:)';
        theta = thetasp(ii);
        sinth = sind(theta);
    end

    % 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 matching 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; theta; theta; theta; theta];
        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_sam

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 % end of function