Added gtIndexFwdSimGrains and gtIndexFwdSimPairs.

Forward simulation of (u,v,w) diffraction spot positions on the detector.

Signed-off-by: Peter Reischig <peter.reischig@esrf.fr>

zUtil_Indexter/gtIndexFwdSimGrains.m

+function grain = gtIndexFwdSimGrains(grain, parameters)
+% GTINDEXFWDSIMGRAINS Forward simulates the spot positions for all 
+% reflections in all grains.
+% 
+% grain = gtIndexFwdSimGrains(grain, parameters)
+%
+% ---------------------------------------------------
+%
+% INPUT
+%   grain      - cell array containing grains as in index.mat
+%   parameters - parameters as in the parameters file 
+%
+% OUTPUT
+%   grain.fsim. - added field containing fwd simulation data:
+%    .cu        - U coordinate of center
+%    .cv        - U coordinate of center
+%    .cw        - W (image no.) coordinate of center
+%    .eta       - eta diffraction angle in degrees
+%    .sinth     - sin(theta)
+%    .plref     - plane normal in reference crystal (unstrained theoretical)
+%    .refind    - linear index of reflection in grain
+%    .thetatype - thetatype
+
+%tic
+
+if ~exist('parameters','var')
+    parameters = load('parameters.mat');
+    parameters = parameters.parameters;
+end
+
+
+labgeo = parameters.labgeo;
+samgeo = parameters.samgeo;
+cryst  = parameters.cryst;
+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 (or Rod2g) 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.cu = uv(1,keep);
+    grain{ii}.fsim.cv = uv(2,keep);
+    grain{ii}.fsim.cw = om(keep)/omstep;
+    
+    grain{ii}.fsim.eta   = gtGeoEtaFromDiffVec(dveclab(:,keep)', labgeo)';
+    grain{ii}.fsim.sinth = sinth4(keep);
+    grain{ii}.fsim.plref = plref4(:,keep);
+    grain{ii}.fsim.refind = refind4(:,keep);
+    grain{ii}.fsim.thetatype = thtype4(keep);
+    %grain{ii}.fsim.dveclab  = reshape(dveclab,3,[],4);
+    %grain{ii}.fsim.dvecsam  = reshape(dvecsam,3,[],4);
+  
+end
+
+%toc
+
+end % of main function
+
+
+

zUtil_Indexter/gtIndexFwdSimPairs.m

+function grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, ...
+                 latticepar, omstep, energy)
+% GTINDEXFWDSIMPAIRS Forward simulates spot positions on the detector 
+% for all spots indexed as Friedel pairs.
+% 
+% grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, ...
+%                            latticepar, omstep, energy)
+%
+% ---------------------------------------------------
+%
+% INPUT
+%   grain       - cell array containing grains as in index.mat
+%   labgeo      - as in parameters.labgeo
+%   samgeo      - as in parameters.samgeo
+%   drift       - sample drifts; if empty, no drifts are considered
+%   latticepar  - lattice parameters (1x6)
+%   omstep      - omega rotational stepping (degrees/image)
+%   energy      - beam energy in keV
+%
+% OUTPUT
+%   grain.fsimA - added field containing fwd simulation data for spots A
+%   grain.fsimB - added field containing fwd simulation data for spots B
+
+
+%%%%%%%%%%%%%%%%%%%
+%% Prepare input
+%%%%%%%%%%%%%%%%%%%
+
+% B matrix for transformation from (hkl) to Cartesian Crystal reference 
+Bmat = gtCrystHKL2CartesianMatrix(latticepar);
+
+% 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);
+
+                                 
+% Rotation tensor describing crystal orientation; calculate all grains 
+% at once for speed-up.
+% (gtCrystRodriguesVector and gtFedRod2RotTensor (or Rod2g) used together 
+% gives the transformation tensor from CRYSTAL TO SAMPLE coordinates)
+Rvec    = gtIndexAllGrainValues(grain,'R_vector',[],1,1:3);
+cry2sam = gtFedRod2RotTensor(Rvec');
+
+
+
+for ii = 1:length(grain)
+     
+    if ~isfield(grain{ii},'strain')
+        grain{ii}.strain.strainT = NaN(3,3);
+    end
+        
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% Theoretical plane normals in deformed state
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    % Coordinates of the plane normals in the Cartesian Crystal system
+    pl_cry          = gtCrystHKL2Cartesian(grain{ii}.hklsp, Bmat);
+    grain{ii}.plref = pl_cry;
+     
+    % Plane normals in the undrifted SAMPLE reference
+    pl_sam = cry2sam(:,:,ii) * pl_cry;
+    
+    % Impose strain on grain
+    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-spacings in the reference crystal for the given reflections
+    dsp_ref = gtCrystDSpacing(grain{ii}.hkl, Bmat);
+
+    % d-spacing after applying the strain tensor
+    dsp_d = drel .* dsp_ref;
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% Apply rotational drift
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    if ~isempty(drift)
+        % In case of drifts, need to calculate the omega angles for A and 
+        % B reflections separately.
+        
+        grain{ii}.fsimA = sfApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
+            grain{ii}.omegaA, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+        
+        grain{ii}.fsimB = sfApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
+            grain{ii}.omegaB, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+        
+    else
+        % In case of no drifts, the omega angles for reflections A and B 
+        % can be calculated together.
+
+        % Bragg angles theta in the deformed state
+        sinth_d = 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]';        
+        pld     = sam2lab*pl_samd;
+        
+        % Predict omega angles: 4 for each plane normal (the Friedel pairs are
+        % om1-om3 and om2-om4)
+        % Do it only once, for A and B at the same time.
+        [om4, pl_lab4, ~, rot4, omind4] = gtFedPredictOmegaMultiple(...
+            pld, sinth_d, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
+        
+        % Check if all reflections occur
+        noref = any(isnan(om4), 1);
+        if any(noref)
+            msg = sprintf('No reflection was predicted for some pairs in grain #%d', ii);
+            error(msg)
+        end
+      
+        
+        % Get fwd simulated data for reflections A and B.
+        % Find the correct omega index for the pairs where the difference is
+        % the smallest (shouldn't matter if grain.omega is the measured one
+        % or the one corrected for the pair).       
+        grain{ii}.fsimA = sfFwdSim(grain{ii}.omegaA, grain{ii}.center', ...
+            sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+        
+        grain{ii}.fsimB = sfFwdSim(grain{ii}.omegaB, grain{ii}.center', ...
+            sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+               
+    end
+  
+    % Extra info
+    if 1
+        grain{ii}.fsimA.drel  = drel;
+        grain{ii}.fsimB.drel  = drel;
+        grain{ii}.fsimA.plsam = pl_samd;
+        grain{ii}.fsimB.plsam = pl_samd;
+    end
+end
+
+
+end % of main function
+
+
+%%%%%%%%%%%%%%%%%%%%%%
+%% SUB-FUNCTIONS
+%%%%%%%%%%%%%%%%%%%%%%
+
+function fsim = sfApplyDrifts(pl, dsp, gc, om, drift, labgeo, ...
+                samgeo, rotcomp, Qdet, omstep, energy)
+% pl    - plane normals before drifts
+% sinth - sin(theta) in deformed state
+% om    - omega-s to interpolate the rotational drift
+% gc    - grain center in Sample reference before drift
+%
+% Use the measured omega position as an approximation to calculate
+% the drifts, because the fwd simulated is not yet known. As the
+% drifts are assumed smooth and slow, this is a good approximation.
+
+    % Plane normals, grain centers and sin(theta) after drifts
+    [pld, gcd, sinthd] = gtFitApplyDrifts(pl, gc, dsp, om, drift, energy);
+    
+    % The plane normals need to be brought in the Lab reference where the 
+    % beam direction and rotation axis are defined.
+    % Use the Sample -> Lab rotation assuming omega=0;
+    % (samgeo.dir norm is always 1, so sam2lab will be preserve vector norms)
+    sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz]';
+    pld     = sam2lab*pld;
+    
+    
+    % Predict omega angles: 4 for each plane normal (the Friedel pairs are
+    % om1-om3 and om2-om4)
+    [om4, pllab4, ~, rot4, omind4] = gtFedPredictOmegaMultiple(...
+        pld, sinthd, labgeo.beamdir', labgeo.rotdir', rotcomp, []);
+    
+    % Check if all reflections occur
+    noref = any(isnan(om4), 1);
+    if any(noref)
+        error('No reflection was predicted for some pairs in grain.');
+    end
+    
+    % Find the correct omega index for the pairs where the difference is
+    % the smallest (shouldn't matter if grain.omega is the measured one
+    % or the one corrected for the pair).
+    
+    % Get fwd simulated data for the observed reflections
+    fsim = sfFwdSim(om, gcd, ...
+           sinthd, om4, pllab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+    
+end
+
+
+function fsim = sfFwdSim(om_mes, gcsam, sinth, om4, pllab4, rot4, ...
+                omind4, labgeo, samgeo, Qdet, omstep)
+    % Provides the fwd simulation data for a given omega index set.
+
+    % Absolute differences
+    dom = om4 - om_mes([1 1 1 1],:);
+    
+    % Find the indices with the smallest absolute difference; consider 
+    % errors and periodicity of 360deg in omega
+    [~, ind] = min( abs( [dom; dom - 360; dom + 360] ) , [], 1);
+    
+    % Predicted omega is 360deg more
+    omhigh = ( (5 <= ind) & (ind <= 8) );
+
+    % Predicted omega is 360deg less
+    omlow  = ( (9 <= ind) & (ind <= 12) );
+    
+    % Correct indices
+    ind(omhigh) = ind(omhigh) - 4;
+    ind(omlow)  = ind(omlow)  - 8;
+
+    % Get final omega angles (can be <0) using the corresponding linear
+    % indices of omind4
+    linind          = ind + (0:length(ind)-1)*4;
+    fsim.om         = om4(linind);
+    fsim.om(omhigh) = fsim.om(omhigh) - 360; 
+    fsim.om(omlow)  = fsim.om(omlow)  + 360; 
+
+    % Omega index 
+    fsim.omind = omind4(linind);
+
+    % Omegas in a column vector
+    ind1 = (ind==1);
+    ind2 = (ind==2);
+    ind3 = (ind==3);
+    ind4 = (ind==4);
+   
+    % Extract the correct pllab vectors according to omind
+    fsim.pllab         = zeros(3, length(ind));
+    fsim.pllab(:,ind1) = pllab4(:,ind1,1);
+    fsim.pllab(:,ind2) = pllab4(:,ind2,2);
+    fsim.pllab(:,ind3) = pllab4(:,ind3,3);
+    fsim.pllab(:,ind4) = pllab4(:,ind4,4);
+
+    % Extract the correct rotation matrices according to omind
+    fsim.rot           = zeros(3, 3, length(ind));
+    fsim.rot(:,:,ind1) = rot4(:,:,ind1,1);
+    fsim.rot(:,:,ind2) = rot4(:,:,ind2,2);
+    fsim.rot(:,:,ind3) = rot4(:,:,ind3,3);
+    fsim.rot(:,:,ind4) = rot4(:,:,ind4,4);
+   
+        
+    % Diffraction vector
+    fsim.dveclab = gtFedPredictDiffVecMultiple(fsim.pllab, labgeo.beamdir');
+    fsim.dvecsam = gtGeoLab2Sam(fsim.dveclab', fsim.om, labgeo, samgeo, 1, 1)';
+
+    % Make sure grain center is extended
+    gcsam = gcsam(:,ones(size(fsim.dveclab,2),1));
+    
+    % Absolute detector coordinates U,V in pixels 
+    uv = gtFedPredictUVMultiple(fsim.rot, fsim.dveclab, gcsam, labgeo.detrefpos', ...
+         labgeo.detnorm', Qdet, [labgeo.detrefu; labgeo.detrefv]);
+
+     
+    % 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 
+    fsim.uvw = [uv; fsim.om/omstep];
+   
+end
+
+