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

for i = 1:size(hklsp,2)
    %pl0(i, :) = gtHKL2Cart(hklsp(:, i)', cryst); 
    % !!! should be changed to:
    pl0(i, :) = gtCrystHKL2Cartesian(hklsp(:,i), cryst.latticepar)';   
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 i=1:size(pl_g,1)

    %hkl(i,:) = hklt(thetatypesp(i),:);
    hkl(:,i) = hklt(:,thetatypesp(i))';
    pl       = pl_g(i,:)';
    sinth    = sind(thetasp(i));
    
    % 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(:,i)'; hkl(:,i)'; hkl(:,i)'; hkl(:,i)'];
        grain.allblobs.hklsp = [grain.allblobs.hklsp; hklsp(:,i)'; -hklsp(:,i)'; hklsp(:,i)'; -hklsp(:,i)'];
        grain.allblobs.theta = [grain.allblobs.theta; thetasp(i); thetasp(i); thetasp(i); thetasp(i)];
        grain.allblobs.thetatype = [grain.allblobs.thetatype; thetatypesp(i); thetatypesp(i); thetatypesp(i); thetatypesp(i)];
        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