gtIndexFwdSimGrains.m

function grain = gtIndexFwdSimGrains(grain, parameters, fsimID)
% GTINDEXFWDSIMGRAINS  Forward simulates the spot positions for all 
%                      reflections in all grains.
%
%     grain = gtIndexFwdSimGrains(grain, [parameters], [fsimID])
%     ----------------------------------------------------------
%     It should be run after gtIndexAddInfo, like below:
%
%       index = load('4_grains/phase_##/index.mat')
%       grain = gtIndexAddInfo(index.grain, phaseid, parameters);
%
%     INPUT
%       grain      = <cell>       structure containing grain information
%       parameters = <struct>     parameters as in the parameters file 
%       fsimID     = <double>     forward simulation id - in case we want to forward simulate
%                                 on more than one detector {1}
%
%     OUTPUT
%       grain.fsim = <struct>     added field containing fwd simulation data:
%                 .cu        = <double>     U coordinate of center
%                 .cv        = <double>     U coordinate of center
%                 .cw        = <double>     W (image no.) coordinate of center
%                 .eta       = <double>     eta diffraction angle in degrees
%                 .sinth     = <double>     sin(theta)
%                 .plref     = <double>     plane normal in reference crystal (unstrained theoretical)
%                 .refind    = <double>     linear index of reflection in grain
%                 .thetatype = <double>     thetatype


if ~exist('parameters','var') || isempty(parameters)
    parameters = load('parameters.mat');
    parameters = parameters.parameters;
end
if ~exist('fsimID','var') || isempty(fsimID)
    fsimID = 1;
end
  
labgeo = parameters.labgeo;
if ~isfield(labgeo, 'detorig')
    labgeo.detorig = gtGeoDetOrig(labgeo);
end
if ~isfield(labgeo, 'detnorm')
    labgeo.detnorm = cross(labgeo.detdiru, labgeo.detidirv);
end

samgeo = parameters.samgeo;
cryst  = parameters.cryst(grain{1}.phaseid);
energy = parameters.acq.energy;
omstep = 180/parameters.acq.nproj;

% Scaling for Lab -> Detector transformation (mm -> pixels)
detscaleu = 1/labgeo.pixelsizeu;
detscalev = 1/labgeo.pixelsizev;

% Rotation tensor components
rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir','col');

% Detector projection matrix
Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru', labgeo.detdirv', ...
                                     detscaleu, detscalev);

% B matrix for transformation from (hkl) to Cartesian Crystal reference 
Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);

% Coordinates of the plane normals in the Cartesian Crystal system
pl_ref  = gtCrystHKL2Cartesian(double(cryst.hklsp), Bmat);
refind0 = 1:size(pl_ref,2);

% D-spacings in the reference crystal for the given reflections
dsp_ref = gtCrystDSpacing(double(cryst.hkl(:,cryst.thetatypesp)), Bmat);

                                 
% Rotation tensor describing crystal orientation; calculate all grains 
% at once for speed-up.
% (gtCrystRodriguesVector and gtFedRod2RotTensor used together 
% gives the transformation tensor from CRYSTAL TO SAMPLE coordinates)
Rvec    = gtIndexAllGrainValues(grain,'R_vector',[],1,1:3);
cry2sam = gtFedRod2RotTensor(Rvec');


% Parallel loop through grains
parfor ii = 1:length(grain)
                    
    % Plane normals in the Sample reference
    pl_sam = cry2sam(:,:,ii) * pl_ref;
    
    % Impose strain on grain
    if ~isfield(grain{ii},'strain')
        grain{ii}.strain.strainT = NaN(3,3);
    end
    
    if any(isnan(grain{ii}.strain.strainT(:)))
        %fprintf('No strain info in grain #%d. Using as measured plane normals.\n',ii)
        pl_samd = pl_sam;
        drel    = ones(1,size(pl_sam,2));
    else
        % New deformed plane normals and relative elongations
        [pl_samd, drel] = gtStrainPlaneNormals(pl_sam, grain{ii}.strain.strainT);
    end
    
    % d-spacing after applying the strain tensor
    dsp_d = drel .* dsp_ref;

    % Bragg angles theta in the deformed state
    sinth = 0.5*gtConvEnergyToWavelength(energy)./dsp_d;
    
    % 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)
    sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz]';
    pl_labd = sam2lab*pl_samd;
    
    
    % Predict omega angles: 4 for each (hkl) plane normal (the Friedel pairs are
    % om1-om3 and om2-om4); -(hkl) plane normals are assumed not to be in the list 
    [om4, pllab4, ~, rot4, omind4] = gtFedPredictOmegaMultiple(...
        pl_labd, sinth, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
    
    % warning('Check function here! Possible bug!')
    
    % Delete those where no reflection occurs
    todel              = (isnan(om4(1,:)) | isnan(om4(3,:)));
    om4(:,todel)       = [];
    omind4(:,todel)    = [];
    sinth(todel)       = [];
    pllab4(:,todel,:)  = [];
    rot4(:,:,todel,:)  = [];
    
    plref  = pl_ref(:,~todel);
    refind = refind0(~todel);
    thtype = cryst.thetatypesp(~todel);
    hklsp  = cryst.hklsp(:,~todel);
    
    om = reshape(om4',[],1)';
    plref4  = [plref, plref, plref, plref];
    sinth4  = [sinth, sinth, sinth, sinth];
    thtype4 = [thtype, thtype, thtype, thtype];
    refind4 = [refind, refind, refind, refind];
    
    % Diffraction vectors in Lab and Sample
    dveclab = gtFedPredictDiffVecMultiple(reshape(pllab4,3,[]), labgeo.beamdir');    
    %dvecsam = gtGeoLab2Sam(dveclab', om, labgeo,...
    %          parameters.samgeo, 1, 1)';
        
    % Etend grain center
    gc_sam = grain{ii}.center';
    gc_sam = gc_sam(:,ones(size(dveclab,2),1));
    
    % Absolute detector coordinates U,V in pixels 
    uv = gtFedPredictUVMultiple(rot4, dveclab, gc_sam, labgeo.detrefpos', ...
         labgeo.detnorm', Qdet, [labgeo.detrefu; labgeo.detrefv])

    % Keep only those spots that fall on the detector 
    keep = (uv(1,:) >= 1) & (labgeo.detsizeu >= uv(1,:)) & ...
           (uv(2,:) >= 1) & (labgeo.detsizev >= uv(2,:));
        
    % Absolute coordinates U,V,W (pixel,pixel,image); to account for 
    % measurement errors, do not correct for the image number
    % to be mod(om,360)/omstep
    
    grain{ii}.fsim(fsimID).cu        = uv(1,keep);
    grain{ii}.fsim(fsimID).cv        = uv(2,keep);
    grain{ii}.fsim(fsimID).cw        = om(keep)/omstep;
    grain{ii}.fsim(fsimID).com       = om(keep);

    grain{ii}.fsim(fsimID).eta       = gtGeoEtaFromDiffVec(dveclab(:,keep)', labgeo)';
    grain{ii}.fsim(fsimID).sinth     = sinth4(keep);
    grain{ii}.fsim(fsimID).plref     = plref4(:,keep);
    grain{ii}.fsim(fsimID).refind    = refind4(:,keep);
    grain{ii}.fsim(fsimID).thetatype = thtype4(keep);
    grain{ii}.fsim(fsimID).dveclab   = reshape(dveclab,3,[],4);
    %grain{ii}.fsim(fsimID).dvecsam   = reshape(dvecsam,3,[],4);

    
end % end parfor loop

%toc

end % of main function