Added 7_fed2 folder with vectorised fwd projection functions.

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

7_fed2/gtFedPredictDiffFacMultiple.m

+function [dfac, clab] = gtFedPredictDiffFacMultiple(Rottens, dlab, csam, ...
+                        detrefpos, detnorm)
+% GTFEDPREDICTDIFFFACMULTIPLE Diffraction vector scalar multiplicator
+% to calculate where the beam hits the detector. Vectorised for speed.
+%
+% [dfac, clab] = gtFedPredictDiffFac(Rottens, dlab, csam, detrefpos, detnorm)
+%
+% -------------------------------------------------------------------
+%
+% Diffraction vector dlab doesn't have to be normalised. 'csam' and 
+% 'detrefpos' must be in identical units.
+%
+% INPUT
+%   Rottens - rotation matrices that transform the Sample coordinates 'csam'
+%             into Lab coordinates (3x3xn or 3x(3xn); for column vectors);
+%             if empty, no rotation will be applied
+%   dlab    - diffraction vectors (diffracted beam direction) in Lab
+%             reference (3xn)
+%   csam    - positions in Sample reference (or in Lab if Rottens is empty)
+%             to be projected along 'dlab'; i.e. usually the grain centers;
+%             (3xn)
+%   detrefpos - any point in the detector plane (3x1)
+%   detnorm - detector plane normal pointing towards beam (3x1; unit vector)
+%
+% OUTPUT
+%   dfac    - scalar multiplicators; NaN where projection point away from detector
+%             (1xn)
+%   clab    - Lab coordinates of 'csam' (3xn)
+%
+
+% Position of centre after rotation (Sample -> Lab coord.)
+if isempty(Rottens)
+    clab = csam;
+else
+    csamt = reshape(csam,1,[]);
+    
+    rc = reshape(reshape(Rottens,3,[]).*csamt([1 1 1],:),9,[]);
+    
+    clab(3,:) = sum(rc([3 6 9],:),1);
+    clab(2,:) = sum(rc([2 5 8],:),1);
+    clab(1,:) = sum(rc([1 4 7],:),1);
+end
+
+% the scalar multiplicator of the diffraction vector
+dfac = detnorm'*(detrefpos(:,ones(1,size(clab,2))) - clab) ./ (detnorm'*dlab);
+
+% beam won't hit the detector if dfac is negative
+dfac(dfac <= 0) = NaN;
+
+end % of function
\ No newline at end of file

7_fed2/gtFedPredictDiffVecMultiple.m

+function dlab = gtFedPredictDiffVecMultiple(pllab, beamdir)
+% GTFEDPREDICTDIFFVECMULTIPLE Diffraction vectors from plane normals and 
+% beam directions in the Lab reference. Vectorised for speed.
+%
+% dlab = gtFedPredictDiffVecMultiple(pllab, beamdir)
+%
+% --------------------------------------------------------------
+% 
+% Returns the diffraction vectors for a set of plane normals and beam 
+% directions. All given in the Lab reference (i.e. turned into the omega 
+% rotational position). Output 'dlab' is NOT normalised!
+%
+% INPUT
+%   pllab   - normalised plane normals (3xn)
+%   beamdir - normalised beam direction(s) (3x1) or (3xn)
+%
+% OUTPUT
+%   dlab    - diffraction vectors (diffracted beam directions);
+%             NOT normalised! (3xn)
+
+ 
+if size(beamdir,1)==3 && size(beamdir,2)==1
+    dotpro = (2*beamdir')*pllab;
+    dlab   = beamdir(:,ones(1,size(pllab,2))) - dotpro([1 1 1],:).*pllab;
+else
+    dotpro = 2*sum(beamdir.*pllab, 1);
+    dlab   = beamdir - dotpro([1 1 1],:).*pllab;
+end
+
+% % Another way could be using sin(theta), but it's inefficient unless 
+% % the sign is known in advance. 
+% if isempty(sinth)
+% ...   
+% else
+%     if size(beamdir,1)==3 && size(beamdir,2)==1
+%         % Need to consider the plane normal pointing opposite to the beam:
+%         pllabsign = -sign(beamdir'*pllab);
+%         pllab     = pllabsign([1 1 1],:).*pllab;
+%         
+%         % Diffraction vector:
+%         dlab = 2*sinth([1 1 1],:).*pllab + beamdir(:,ones(1,size(pllab,2))) ;
+%         
+%     else
+%         % Need to consider the plane normal pointing opposite to the beam:
+%         pllabsign = -sign(sum(beamdir.*pllab,1));
+%         pllab     = pllabsign([1 1 1],:).*pllab;
+%         
+%         % Diffraction vector:
+%         dlab = 2*sinth([1 1 1],:).*pllab + beamdir ;
+%     end
+% end
+
+
+end % end of function
\ No newline at end of file

7_fed2/gtFedPredictOmegaMultiple.m

+function [om, plrot, plsigned, rot, omind] = gtFedPredictOmegaMultiple(...
+         pl, sinth, beamdir, rotdir, rotcomp, omind)
+% GTFEDPREDICTOMEGAMULTIPLE Omega angles where reflection occur from the 
+% the plane normals (Bragg's law is satisifed). Vectorised for speed.
+% 
+% [om, pllab, plsigned, rot, omind] = gtFedPredictOmegaMultiple(...
+%     pl, sinth, beamdir, rotdir, rotcomp, omind)
+%
+% -------------------------------------------------------------------
+%
+% Predicts the four omega angles where diffraction occurs for multiple 
+% plane normals 'pl '. Rotating 'pl' around 'rotdir' will bring them 
+% in a position where the angle between 'pl' and 'beamdir' will be theta 
+% (the Bragg angle). 
+% All inputs have to be given in the same reference. 
+%
+% The output 'plsigned' is as follows: 
+%   plsigned(:,:,1:2) =  pl
+%   plsigned(:,:,3:4) = -pl
+% These vectors result in rotated plane normals such that their dot products 
+% with the beam direction is always negative, i.e they satisfy: 
+%   plrot(:,i,j) = rot(:,:,i,j)*plsigned(:,i,j)  , 1<=j<=4          
+%   beamdir'*plrot(:,i,j) = -sinth
+%
+% The omega angles are arranged in the output as follows:
+%   - om(1,:)-om(3,:) are one Friedel pair, and om(2,:)-om(4,:) are the 
+%     other Friedel pair in an orthogonal setup ('beamdir' perpendicular
+%     to 'rotdir')
+%   - om1 is always smaller than om2
+% Therefore the omega indices are NOT necessarily ORDERED and are returned 
+% in 'omind'. These indices can be re-used later as input for speed-up, so 
+% that only the relevant omega angles are (re)computed.
+%
+%
+% INPUT
+%    pl      - normalised unit plane normals (3xn)
+%    sinth   - sin(theta) corresponding to 'pl' (1x1 or 1xn)
+%    beamdir - normalised beam direction (3x1 or 3x1)
+%    rotdir  - normalised rotation axis direction (right-handed) (3x1)
+%    rotcomp - rotational matrix components for COLUMN vectors !
+%              (pl_rotated = Srot*pl)
+%    omind   - omega index (1..4) to be computed; empty to get 
+%              all four; (1x1 or 1xn)
+%
+% OUTPUT
+%    om       - the four omegas in degrees where reflections appear;
+%               contains NaN where no solution exists (1xn or 4xn)
+%    plrot    - the plane normals rotated in the diffraction positions;
+%               no. 3 and 4 are inverted after rotation (pointing oppposite
+%               to beam) to ease calculation at a later stage;
+%               size (3xn or 3xnx4)
+%    plsigned - the original input plane normals with no. 3 and 4 inverted;
+%                 plrot(:,i,j) = rot(:,:,i,j)*plsigned(:,i,j) the  is where (3xn or 3xnx4)               
+%    rot      - rotation matrices (3x3xn or 3x3xnx4)
+%    omind    - omega indices (1xn or 4xn) 
+%
+
+
+% % Polychromatic case
+% if rotdir == [0 0 0]
+% 	
+% 	om = 0;
+% 	plrot = -pl*sign(pl'*beamdir);
+% 	plsigned = plrot;
+% 	
+% 	return
+% end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Quadratic equation for omega angles
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+nn = size(pl,2);
+
+if size(beamdir,1)==3 && size(beamdir,2)==1
+    A = beamdir'*(rotcomp.cos*pl);
+    B = beamdir'*(rotcomp.sin*pl);
+    C = beamdir'*(rotcomp.const*pl);
+else
+    A = sum(beamdir.*(rotcomp.cos*pl),1);
+    B = sum(beamdir.*(rotcomp.sin*pl),1);
+    C = sum(beamdir.*(rotcomp.const*pl),1);
+end
+
+D = A.^2 + B.^2 ;
+
+
+if ~isempty(omind)
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% Single omega index
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    if length(omind)==1
+        omind = omind(1,ones(1,nn));
+    end
+        
+    ssp = ((omind==1) | (omind==2));
+    ss  = ssp - ~ssp;
+    ss3 = ss([1 1 1],:);
+
+    CC = C + ss.*sinth;
+    DD = D - CC.^2;
+    E  = sqrt(DD);
+    ok = DD > 0;
+   
+    ssp = ((omind==2) | (omind==3));
+    ss  = ssp - ~ssp;
+    ss  = ss.*ok;
+    
+    om = 2*atand((-B + ss.*E) ./ (CC - A));
+    om(~ok) = NaN;
+            
+    % Limit range
+    om = mod(om,360);
+    
+    % ROTATION MATRICES AND PLANE NORMALS IN LAB
+    
+    % Get rotation matrix and multpily the input plane normals to get
+    % them in the diffracting position
+    rot = gtMathsRotationTensor(om, rotcomp);
+       
+    % expand for multiplication
+    plt = reshape(pl,1,[]);
+    plt = plt([1 1 1],:);
+    
+    plrot = reshape(rot,3,[]).*plt;
+    plrot = plrot(:,1:3:end) + plrot(:,2:3:end) + plrot(:,3:3:end);
+    plrot = reshape(plrot,3,[]);
+    
+    plrot    = ss3.*plrot;
+    plsigned = ss3.*pl;
+    
+else
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %% All four omega indices
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    om = NaN(4,nn);
+    
+    % Diffraction condition 1
+    % Pl_rot and beamdir dotproduct equals -sin(th)
+    CC = C + sinth;
+    DD = D - CC.^2;
+    ok = DD > 0;
+    E  = sqrt(DD(ok));
+    
+    om(1,ok) = 2*atand((-B(ok)-E)./(CC(ok)-A(ok)));
+    om(2,ok) = 2*atand((-B(ok)+E)./(CC(ok)-A(ok)));
+    
+    
+    % Diffraction condition 2
+    % Pl_rot and beamdir dotproduct equals +sin(th)
+    CC = C - sinth;
+    DD = D - CC.^2;
+    ok = DD > 0;
+    E  = sqrt(DD(ok));
+    
+    om(3,ok) = 2*atand((-B(ok)+E)./(CC(ok)-A(ok)));
+    om(4,ok) = 2*atand((-B(ok)-E)./(CC(ok)-A(ok)));
+    
+    
+    % Limit range
+    om = mod(om,360);
+    
+    % Omegea indices
+    omind = [1; 2; 3; 4];
+    omind = omind(:,ones(1,nn));
+
+    
+    % For output conventions:
+    % make sure om1 is smaller then om2
+    chom = om(1,:) > om(2,:);
+    om(1:2,chom)    = [om(2,chom); om(1,chom)];
+    omind(1:2,chom) = [omind(2,chom); omind(1,chom)];
+   
+    
+    % ROTATION MATRICES AND PLANE NORMALS IN LAB
+
+    plsigned        = zeros(3,nn,4);
+    plsigned(:,:,1) = pl;
+    plsigned(:,:,2) = pl;
+    plsigned(:,:,3) = -pl;
+    plsigned(:,:,4) = -pl;
+    
+    % Get rotation matrices and multpily the input plane normals to get
+    % them in the diffracting position (avoid looping)
+    rot(:,:,:,4) = gtMathsRotationTensor(om(4,:), rotcomp);
+    rot(:,:,:,3) = gtMathsRotationTensor(om(3,:), rotcomp);
+    rot(:,:,:,2) = gtMathsRotationTensor(om(2,:), rotcomp);
+    rot(:,:,:,1) = gtMathsRotationTensor(om(1,:), rotcomp);
+    
+    
+    plt = reshape(pl,1,[]);
+    plt = plt([1 1 1],:);
+    
+    plrot = reshape(rot,3,[]).*[plt, plt, plt, plt];
+    
+    plrot = plrot(:,1:3:end) + plrot(:,2:3:end) + plrot(:,3:3:end);
+    
+    plrot = reshape(plrot,3,nn,[]);
+      
+    
+    % Exchange om3 and om4 if needed, so that om1-om3 and om2-om4 are the
+    % opposite reflections (Friedel pairs in case the rotation axis is
+    % perpendicular to the beam).
+    % Midplane of setup is defined by beam direction and rotation axis.
+    % This is needed to consider consistently all possible setups where 
+    % the beam and rotation axis are not perpendicular.
+    
+    % Cross product beamdir and rotdir
+    br = [beamdir(2,:)*rotdir(3) - beamdir(3,:)*rotdir(2); ...
+          beamdir(3,:)*rotdir(1) - beamdir(1,:)*rotdir(3); ...
+          beamdir(1,:)*rotdir(2) - beamdir(2,:)*rotdir(1)];
+    
+    % Dot product of br and pllab
+    if size(beamdir,1)==3 && size(beamdir,2)==1
+        dot1 = br'*plrot(:,:,1);
+        dot3 = br'*plrot(:,:,3);
+    else
+        dot1 = sum(br.*plrot(:,:,1),1);
+        dot3 = sum(br.*plrot(:,:,3),1);
+    end
+    
+    % Indices to be swapped
+    chom = (((dot1.*dot3) > 0) & ~isnan(dot1) & ~isnan(dot3)) ;
+    
+    % Swap 3rd and 4th
+    om([3,4],chom)      = om([4,3],chom);
+    omind([3,4],chom)   = omind([4,3],chom);
+    plrot(:,chom,[3,4]) = plrot(:,chom,[4,3]);
+    rot(:,:,chom,[3,4]) = rot(:,:,chom,[4,3]);
+    
+    
+    % Change all pllab to opposite direction
+    plrot(:,:,[3,4]) = -plrot(:,:,[3,4]);
+    
+    
+end
+
+
+end % of function
\ No newline at end of file

7_fed2/gtFedPredictUVMultiple.m

+function uv = gtFedPredictUVMultiple(Rottens, dlab, csam, detrefpos, ...
+                                     detnorm, Qdet, detrefuv)
+% GTFEDPREDICTUVMULTIPLE Predicts (u,v) peak positions on the detector
+% plane. Vectorised for speed.
+% 
+% uv = gtFedPredictUVMultiple(Rottens, dlab, csam, detrefpos, detnorm, ...
+%      Qdet, detrefuv)
+%
+% -------------------------------------------------------------
+%
+% Predicts peak positions [u,v] on the detector plane.
+% Note:
+%   Diffraction vectors 'dlab' don't have to be normalised.
+%   'csam','detrefpos','detrefuv' must be in identical units.
+%   Qdet is in units that scales lab units into detector pixels.
+%
+% INPUT
+%   Rottens - rotation matrices that transform the Sample coordinates 'csam'
+%             into Lab coordinates (3x3xn or 3x(3xn); for column vectors);
+%             if empty, no rotation will be applied
+%   dlab    - diffraction vectors (diffracted beam direction) in Lab
+%             reference (3xn)
+%   csam    - positions in Sample reference (or in Lab if Rottens is empty)
+%             to be projected along 'dlab'; i.e. usually the grain centers;
+%             (3xn)
+%   detrefpos - any point in the detector plane (3x1)
+%   detnorm   - detector plane normal pointing towards beam (3x1; unit vector)
+%   Qdet      - detector projection matrix (2x3)
+%   detrefuv  - detector reference position (u,v) (2x1)
+%
+% OUTPUT
+%   uv      - (u,v) positions on the detector plane (2xn)
+
+% Multiplicator for dlab
+%  'dfac' is NaN where the beam is diffracted away from the detector plane,
+%  i.e. where the spot is not detected
+[dfac, clab] = gtFedPredictDiffFacMultiple(Rottens, dlab, csam, detrefpos, detnorm);
+
+% u,v position on detector plane
+uv = detrefuv(:,ones(1,size(dlab,2))) + Qdet*(clab + dfac([1 1 1],:).*dlab - ...
+     detrefpos(:,ones(1,size(dlab,2))));
+
+end
\ No newline at end of file

7_fed2/gtFedPredictUVWMultiple.m

+function uvw = gtFedPredictUVWMultiple(Rottens, dlab, csam, detrefpos, ...
+                                       detnorm, Qdet, detrefuv, om, omstep)
+% GTFEDPREDICTUVWMULTIPLE Predicts (u,v,w) peak positions on the detector
+% plane. Vectorised for speed.
+% 
+% uv = gtFedPredictUVWMultiple(Rottens, dlab, csam, detrefpos, detnorm, ...
+%      Qdet, detrefuv, om, omstep)
+%
+% -------------------------------------------------------------
+%
+% Predicts peak positions [u,v] on the detector plane.
+% Note:
+%   Diffraction vectors 'dlab' don't have to be normalised.
+%   'csam','detrefpos','detrefuv' must be in identical units.
+%   Qdet is in units that scales lab units into detector pixels.
+%   'om' is in degrees; 'omstep' is degrees/image
+%
+% INPUT
+%   Rottens - rotation matrices that transform the Sample coordinates 'csam'
+%             into Lab coordinates (3x3xn or 3x(3xn); for column vectors);
+%             if empty, no rotation will be applied
+%   dlab    - diffraction vectors (diffracted beam direction) in Lab
+%             reference (3xn)
+%   csam    - positions in Sample reference (or in Lab if Rottens is empty)
+%             to be projected along 'dlab'; i.e. usually the grain centers;
+%             (3xn)
+%   detrefpos - any point in the detector plane (3x1)
+%   detnorm   - detector plane normal pointing towards beam (3x1; unit vector)
+%   Qdet      - detector projection matrix (2x3)
+%   detrefuv  - detector reference position (u,v) (2x1)
+%   om        - omega positions in degrees (1xn)
+%   omstep    - omega step size in degrees/image
+%
+% OUTPUT
+%   uvw      - (u,v,w) positions on the detector plane;
+%              [pixel, pixel, image no.] (3xn);
+
+uv = gtFedPredictUVMultiple(Rottens, dlab, csam, detrefpos, detnorm, ...
+                            Qdet, detrefuv);
+
+uvw = [ uv; mod(om,360)/omstep ];
+
+end 
\ No newline at end of file

7_fed2/gtFedRod2RotTensor.m

+function rot = gtFedRod2RotTensor(r)
+% GTFEDROD2ROTTENSOR Rotation tensors from Rodrigues vectors. Vectorised
+% for speed.
+% 
+% rot = gtFedRod2RotTensor(r)
+%
+% -----------------------------------------------------
+%
+% Transforms a set of Rodrigues vectors into rotation tensors. It uses the 
+% general mathematical definition, that is right-handed rotations 
+% around 'r' with alpha according to norm(r) = tan(alpha/2) . I.e. it does
+% NOT imply any crystallographic notation or definitions. 
+% 
+% INPUT
+%   r   = Rodrigues vectors (1x3 or 3x1 or 3xn)
+%
+% OUTPUT
+%   rot = rotation matrices (3x3xn)
+%
+
+
+% If one single entry
+if (numel(r) == 3)
+    
+    rs1 = r(1)^2;
+    rs2 = r(2)^2;
+    rs3 = r(3)^2;
+    
+    r1r2 = r(1)*r(2);
+    r1r3 = r(1)*r(3);
+    r2r3 = r(2)*r(3);
+    
+    rot = [rs1 - rs2 - rs3 + 1, 2*(r1r2 - r(3)),     2*(r1r3 + r(2)); ...
+           2*(r1r2 + r(3)),     rs2 - rs1 - rs3 + 1, 2*(r2r3 - r(1)); ...
+           2*(r1r3 - r(2)),     2*(r2r3 + r(1)),     rs3 - rs1 - rs2 + 1];
+    
+    rot = rot/(rs1 + rs2 + rs3 + 1);
+    
+    return
+end
+
+
+% Otherwise vectorize computation
+n   = size(r,2);
+rot = zeros(9, n, class(r));
+
+
+% Diagonals (1,1) (2,2) (3,3)
+r2 = r.^2;
+rot(1,:) = r2(1,:);
+rot(5,:) = r2(2,:);
+rot(9,:) = r2(3,:);
+
+% Off diagonals (1,2) (2,1)
+rr = r(1,:).*r(2,:);
+rot(2,:) = rr + r(3,:);
+rot(4,:) = rr - r(3,:);
+
+% Off diagonals (1,3) (3,1)
+rr = r(1,:).*r(3,:);
+rot(3,:) = rr - r(2,:);
+rot(7,:) = rr + r(2,:);
+
+% Off diagonals (2,3) (3,2)
+rr = r(2,:).*r(3,:);
+rot(6,:) = rr + r(1,:);
+rot(8,:) = rr - r(1,:);
+
+
+% Multiply all by two
+rot = rot*2;
+
+
+% Add term to diagonals
+r2s = sum(r2,1);
+rot(1,:) = rot(1,:) - r2s + 1;
+rot(5,:) = rot(5,:) - r2s + 1;
+rot(9,:) = rot(9,:) - r2s + 1;
+
+
+% Divide all
+r2s = 1 + r2s;
+rot = rot./r2s(ones(9,1),:);
+
+
+% Reshape for output
+rot = reshape(rot, 3, 3, []);
+
+
+end % of function
\ No newline at end of file