Shape-functions: updated NW sf to use diffractometer

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Deformation/gtDefFwdProjGvdm2NW.m

@@ -11,25 +11,47 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
 
     o.fprintf('Forward projection (%s):\n', mfilename)
 
+    om_step = gtAcqGetOmegaStep(parameters, det_ind);
     nbl = numel(ref_sel);
+
     uinds = gv.used_ind;
     nv = numel(uinds);
 
+    use_diffractometer = (isfield(parameters, 'diffractometer') ...
+        && numel(parameters.diffractometer) >= det_ind ...
+        && det_ind > 1);
+
     labgeo = parameters.labgeo;
     samgeo = parameters.samgeo;
-    om_step = gtAcqGetOmegaStep(parameters, det_ind);
+    if (use_diffractometer)
+        diff_acq = parameters.diffractometer(det_ind);
+        diff_ref = parameters.diffractometer(1);
 
-    g = gtMathsRod2OriMat(grain.R_vector);
+        rotcomp = gtMathsRotationMatrixComp(diff_acq.axes_rotation', 'col');
+        rotdir = diff_acq.axes_rotation';
+
+        instr_t = gtGeoDiffractometerTensor(diff_acq, 'sam2lab', ...
+            'reference_diffractometer', diff_ref);
+    else
+        rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
+        rotdir = labgeo.rotdir';
+    end
 
     if (strcmpi(fedpars.defmethod, 'rod_rightstretch'))
+        % In this case, the deformation is relative to [0, 0, 0], so we
+        % need crystal plane normals
         if (isfield(grain.allblobs, 'plcry'))
             pls_orig = grain.allblobs(det_ind).plcry(ref_sel, :);
         else
             % We bring the plane normals back to the status of pl_cry
             pls_orig = grain.allblobs(det_ind).plorig(ref_sel, :);
+
+            g = gtMathsRod2OriMat(grain.R_vector);
             pls_orig = gtVectorLab2Cryst(pls_orig, g);
         end
     else
+        % Here the deformation is relative to the average oriantion, so we
+        % need the undeformed sample plane normals
         pls_orig = grain.allblobs(det_ind).plorig(ref_sel, :);
     end
 
@@ -38,7 +60,7 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
         && ~isempty(fedpars.detector(det_ind).blobs_w_interp))
         blobs_w_interp = fedpars.detector(det_ind).blobs_w_interp;
         if (numel(blobs_w_interp) ~= nbl)
-            error('gtDefFwdProjGvdm2W:wrong_argument', ...
+            error([mfilename ':wrong_argument'], ...
                 'Number of "blobs_w_interp" (%d) doesn''t match with the number of blobs (%d)', ...
                 numel(blobs_w_interp), nbl)
         end
@@ -49,11 +71,20 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
     sinths = grain.allblobs(det_ind).sintheta(ref_sel);
     ominds = grain.allblobs(det_ind).omind(ref_sel);
 
+    if (~use_diffractometer)
+        % 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)
+        pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true);
+    end
+
     if (isfield(fedpars, 'ref_omega') && ~isempty(fedpars.ref_omega))
         ref_oms = fedpars.ref_omega(ref_sel);
     else
         ref_oms = grain.allblobs(det_ind).omega(ref_sel);
     end
+
     w_shifts = round(ref_oms ./ om_step ./ blobs_w_interp);
 
     if (isfield(fedpars, 'ref_eta') && ~isempty(fedpars.ref_eta))
@@ -61,13 +92,8 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
     else
         ref_ns = grain.allblobs(det_ind).eta(ref_sel);
     end
-    n_shifts = round(ref_ns ./ eta_step);
 
-    % 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)
-    pls_orig = gtGeoSam2Lab(pls_orig, eye(3), labgeo, samgeo, true);
+    n_shifts = round(ref_ns ./ eta_step);
 
     linear_interp = ~(isfield(fedpars, 'projector') ...
         && strcmpi(fedpars.projector, 'nearest'));
@@ -75,6 +101,7 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
     gvpow = gv.pow(1, uinds);
 
     gvd = gv.d(:, uinds);
+
     % Deformation tensor (relative to reference state) from its components
     defT = gtFedDefTensorFromComps(gvd, fedpars.dcomps, fedpars.defmethod, 0);
 
@@ -96,16 +123,25 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
 
         % New deformed plane normals and relative elongations (relative to
         % reference state)
-        pl_orig = pls_orig(ii_b, :)';  % ! use plcry and not plsam to keep omega order below!
+        % !!! use plcry and not plsam to keep omega order below !!!
+        pl_orig = pls_orig(ii_b, :)';
 
         [pl_samd, drel] = gtStrainPlaneNormals(pl_orig, defT); % unit column vectors
 
+        if (use_diffractometer)
+            % The plane norma ls 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)
+            pl_samd = gtGeoSam2Lab(pl_samd', instr_t, labgeo, samgeo, true)';
+        end
+
         % New sin(theta)
         sinth = sinths(ii_b) ./ drel;
 
         % Predict omega angles: 4 for each plane normal
-        [om, pl_lab] = gtFedPredictOmegaMultiple(pl_samd, sinth, labgeo.beamdir', ...
-            labgeo.rotdir', labgeo.rotcomp, ominds(ii_b));
+        [om, pl_lab] = gtFedPredictOmegaMultiple( ...
+            pl_samd, sinth, labgeo.beamdir', rotdir, rotcomp, ominds(ii_b));
 
         valid_voxels(:, ii_b) = ~isnan(om);
 
@@ -124,9 +160,11 @@ function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_s
         w_bl = w_bl / blobs_w_interp(ii_b);
         w_bl = w_bl - w_shifts(ii_b);
 
-        n_bl = gtGeoEtaFromDiffVec(pl_lab', parameters.labgeo)';
-        n_bl = gtGrainAnglesTabularFix360deg(n_bl, ref_ns(ii_b));
-        n_bl = n_bl ./ eta_step - n_shifts(ii_b);
+        eta = gtGeoEtaFromDiffVec(pl_lab', parameters.labgeo)';
+        eta = gtGrainAnglesTabularFix360deg(eta, ref_ns(ii_b));
+
+        n_bl = eta ./ eta_step;
+        n_bl = n_bl - n_shifts(ii_b);
 
         nw{ii_b} = [n_bl; w_bl];
 

zUtil_Deformation/gtDefFwdProjGvdm2W.m

@@ -60,7 +60,7 @@ function bl = gtDefFwdProjGvdm2W(grain, ref_sel, gv, fedpars, parameters, det_in
         && ~isempty(fedpars.detector(det_ind).blobs_w_interp))
         blobs_w_interp = fedpars.detector(det_ind).blobs_w_interp;
         if (numel(blobs_w_interp) ~= nbl)
-            error('gtDefFwdProjGvdm2W:wrong_argument', ...
+            error([mfilename ':wrong_argument'], ...
                 'Number of "blobs_w_interp" (%d) doesn''t match with the number of blobs (%d)', ...
                 numel(blobs_w_interp), nbl)
         end