Shape-functions: introduced fwd-projected sf generation support

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

zUtil_Deformation/GtGrainODFuvwSolver.m

@@ -21,10 +21,10 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
             self.chose_voxel_centers();
 
             num_vox_centers = size(self.voxel_centers, 1);
+            recgeo = self.parameters.recgeo(self.sampler.detector_index);
 
             switch (self.shape_functions_type)
                 case 'nw2uvw'
-                    recgeo = self.parameters.recgeo(self.sampler.detector_index);
 
                     sf_nw = gtDefShapeFunctionsCreateNW(self.sampler, 2, ...
                         'data_type', self.data_type);
@@ -52,10 +52,10 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
                             fprintf('(%3d/%03d) ', ii_v, num_vox_centers);
                             disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
 
-                            self.shape_functions{ii_v} = gtDefShapeFunctionsFwdProjUVW( ...
-                                disp_sampler, 'factor', 15, ...
-                                'recenter_sf', false, ...
-                                'data_type', self.data_type);
+                            self.shape_functions{ii_v} = gtDefShapeFunctionsFwdProj( ...
+                                disp_sampler, 'recenter_sf', false, ...
+                                'data_type', self.data_type, ...
+                                'rspace_voxel_size', recgeo.voxsize * self.rspace_downscaling);
                         end
                         sfs = self.shape_functions;
                         save([shape_functions_name '.mat'], 'sfs', '-v7.3')

zUtil_Deformation/gtDefFwdProjGvdm.m

@@ -68,7 +68,8 @@ function bl = gtDefFwdProjGvdm(grain, ref_sel, gv, fedpars, parameters, det_ind,
         uvw_shifts = round(or_uvw)';
     end
 
-    linear_interp = ~(isfield(fedpars, 'projector') && strcmpi(projector, 'nearest'));
+    linear_interp = ~(isfield(fedpars, 'projector') ...
+        && strcmpi(fedpars.projector, 'nearest'));
 
     num_oversampling = size(gv.d, 3);
     gvpow = gv.pow(1, uinds);

zUtil_Deformation/gtDefFwdProjGvdm2NW.m

+function bl = gtDefFwdProjGvdm2NW(grain, ref_sel, gv, fedpars, parameters, eta_step, det_ind, verbose)
+
+    if (~exist('det_ind', 'var'))
+        det_ind = 1;
+    end
+    if (~exist('verbose', 'var'))
+        verbose = true;
+    end
+
+    o = GtConditionalOutput(verbose);
+
+    o.fprintf('Forward projection (%s):\n', mfilename)
+
+    nbl = numel(ref_sel);
+    uinds = gv.used_ind;
+    nv = numel(uinds);
+
+    labgeo = parameters.labgeo;
+    samgeo = parameters.samgeo;
+    om_step = gtGetOmegaStepDeg(parameters, det_ind);
+
+    g = gtMathsRod2OriMat(grain.R_vector);
+
+    if (strcmpi(fedpars.defmethod, 'rod_rightstretch'))
+        if (isfield(grain.allblobs, 'plcry'))
+            pls_orig = grain.allblobs.plcry(ref_sel, :);
+        else
+            % We bring the plane normals back to the status of pl_cry
+            pls_orig = grain.allblobs.plorig(ref_sel, :);
+            pls_orig = gtVectorLab2Cryst(pls_orig, g);
+        end
+    else
+        pls_orig = grain.allblobs.plorig(ref_sel, :);
+    end
+
+    sinths = grain.allblobs.sintheta(ref_sel);
+    ominds = grain.allblobs.omind(ref_sel);
+    ref_ws = grain.allblobs.detector(det_ind).uvw(ref_sel, 3);
+    w_shifts = round(ref_ws);
+    ref_ns = grain.allblobs.eta(ref_sel);
+    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, false);
+
+    linear_interp = ~(isfield(fedpars, 'projector') ...
+        && strcmpi(fedpars.projector, 'nearest'));
+
+    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);
+
+    %%% Computation of indices
+    o.fprintf('   * Computing indices and bbsizes: ')
+    t = tic();
+
+    nw = cell(nbl, 1);
+    nw_min = zeros(nbl, 2);
+    nw_max = zeros(nbl, 2);
+
+    valid_voxels = false(nv, nbl);
+
+    for ii_b = 1:nbl
+        num_chars = o.fprintf('%03d/%03d', ii_b, nbl);
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Calculate new detector coordinates
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+        % 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!
+
+        [pl_samd, drel] = gtStrainPlaneNormals(pl_orig, defT); % unit column vectors
+
+        % 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));
+
+        valid_voxels(:, ii_b) = ~isnan(om);
+
+        % Delete those where no reflection occurs
+        if (any(isnan(om)))
+            inds_bad = find(isnan(om));
+            gvd(:, inds_bad(1:min(10, numel(inds_bad))))
+            warning('gtFedFwdProjExact:bad_R_vectors', ...
+                'No diffraction from some elements after deformation (%d over %d) for blob %d.', ...
+                numel(inds_bad), numel(om), ii_b)
+        end
+
+        nw_bl = om ./ om_step;
+        nw_bl = gtGrainAnglesTabularFix360deg(nw_bl, ref_ws(ii_b), parameters);
+        nw_bl = nw_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);
+
+        nw{ii_b} = [n_bl; nw_bl];
+
+        nw_min(ii_b, :) = min(nw{ii_b}, [], 2);
+        nw_max(ii_b, :) = max(nw{ii_b}, [], 2);
+
+        o.fprintf(repmat('\b', [1, num_chars]));
+    end
+    o.fprintf('Done in %g s\n', toc(t));
+
+    o.fprintf('   * Computing max BBox size and feasibilities:\n')
+
+    if (linear_interp)
+        blob_nw_sizes = ceil(nw_max) - floor(nw_min) + 1;
+        blob_orig_nw_shifts = 1 - floor(nw_min);
+    else
+        blob_nw_sizes = round(nw_max) - round(nw_min) + 1;
+        blob_orig_nw_shifts = 1 - round(nw_min);
+    end
+
+    o.fprintf('     Computed Blob N sizes: [%s]\n', sprintf(' %d', blob_nw_sizes(:, 1)));
+    o.fprintf('     Computed Blob W sizes: [%s]\n', sprintf(' %d', blob_nw_sizes(:, 2)));
+
+    bl = gtFwdSimBlobDefinition('sf_nw', nbl);
+
+    %%% Blob projection
+    o.fprintf('   * Projecting volumes: ')
+    t = tic();
+    for ii_b = 1:nbl
+        num_chars = o.fprintf('%03d/%03d', ii_b, nbl);
+
+        % Detector coordinates U,V in blob
+        nw_bl = nw{ii_b} + blob_orig_nw_shifts(ii_b * ones(1, nv), :)';
+
+        % Let's now filter valid voxels
+        nw_bl = nw_bl(:, valid_voxels(:, ii_b))';
+        gvpow_v = reshape(gvpow(valid_voxels(:, ii_b)), [], 1);
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Sum intensities
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Distribute value over 8 pixels
+
+        if (linear_interp)
+            [inds_nw, ints_nw] = gtMathsGetInterpolationIndices(nw_bl, ...
+                gvpow_v, fedpars.bltype);
+        else
+            inds_nw = round(nw_bl);
+            ints_nw = gvpow_v;
+        end
+
+        % Accumulate all intensities in the blob voxel-wise
+        %   'uvw' needs to be nx1; 'ints' is now 1xnx8
+        try
+            bl(ii_b).intm = accumarray(inds_nw, ints_nw, blob_nw_sizes(ii_b, :));
+        catch mexc
+            fprintf(['\n\nERROR\nThe error is probably caused by the' ...
+                ' projected intensities falling outside the blob volume.' ...
+                '\nTry increasing the blob volume padding:' ...
+                ' fedpars.detector(det_ind).blobsizeadd, or FIXING' ...
+                ' YOUR CODE!!\n\n'])
+            disp('Blob ID:')
+            disp(ii_b)
+            disp('Blob size:')
+            disp(blob_nw_sizes(ii_b))
+            disp('Min projected NW coordinate:')
+            fprintf('\t%g\t%g \t(rounded: %d %d )\n', ...
+                min(nw_bl, [], 1), min(inds_nw, [], 1))
+            disp('Max projected W coordinate:')
+            fprintf('\t%g \t(rounded: %d )\n', ...
+                max(nw_bl, [], 1), max(inds_nw, [], 1))
+            new_exc = GtFedExceptionFwdProj(mexc, det_ind, [0 0 0]);
+            throw(new_exc)
+        end
+
+        bl(ii_b).bbsize = blob_nw_sizes(ii_b, :);
+
+        bl(ii_b).bbnim = ([nw_min(ii_b, 1), nw_max(ii_b, 1)] + n_shifts(ii_b)) .* eta_step;
+        im_w_low_lim = w_shifts(ii_b) - blob_orig_nw_shifts(ii_b, 2) + 1;
+        bl(ii_b).bbwim = [im_w_low_lim, im_w_low_lim + bl(ii_b).bbsize(2) - 1];
+
+        o.fprintf(repmat('\b', [1, num_chars]));
+    end
+    o.fprintf('Done in %g s\n', toc(t));
+end

zUtil_Deformation/gtDefFwdProjGvdm2W.m

+function bl = gtDefFwdProjGvdm2W(grain, ref_sel, gv, fedpars, parameters, det_ind, verbose)
+
+    if (~exist('det_ind', 'var'))
+        det_ind = 1;
+    end
+    if (~exist('verbose', 'var'))
+        verbose = true;
+    end
+
+    o = GtConditionalOutput(verbose);
+
+    o.fprintf('Forward projection (%s):\n', mfilename)
+
+    nbl = numel(ref_sel);
+    uinds = gv.used_ind;
+    nv = numel(uinds);
+
+    labgeo = parameters.labgeo;
+    samgeo = parameters.samgeo;
+    om_step = gtGetOmegaStepDeg(parameters, det_ind);
+
+    g = gtMathsRod2OriMat(grain.R_vector);
+
+    if (strcmpi(fedpars.defmethod, 'rod_rightstretch'))
+        if (isfield(grain.allblobs, 'plcry'))
+            pls_orig = grain.allblobs.plcry(ref_sel, :);
+        else
+            % We bring the plane normals back to the status of pl_cry
+            pls_orig = grain.allblobs.plorig(ref_sel, :);
+            pls_orig = gtVectorLab2Cryst(pls_orig, g);
+        end
+    else
+        pls_orig = grain.allblobs.plorig(ref_sel, :);
+    end
+
+    sinths = grain.allblobs.sintheta(ref_sel);
+    ominds = grain.allblobs.omind(ref_sel);
+    ref_ws = grain.allblobs.detector(det_ind).uvw(ref_sel, 3);
+    w_shifts = round(ref_ws);
+
+    % 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, false);
+
+    linear_interp = ~(isfield(fedpars, 'projector') ...
+        && strcmpi(fedpars.projector, 'nearest'));
+
+    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);
+
+    %%% Computation of indices
+    o.fprintf('   * Computing indices and bbsizes: ')
+    t = tic();
+
+    w = cell(nbl, 1);
+    w_min = zeros(nbl, 1);
+    w_max = zeros(nbl, 1);
+
+    valid_voxels = false(nv, nbl);
+
+    for ii_b = 1:nbl
+        num_chars = o.fprintf('%03d/%03d', ii_b, nbl);
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Calculate new detector coordinates
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+        % 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!
+
+        [pl_samd, drel] = gtStrainPlaneNormals(pl_orig, defT); % unit column vectors
+
+        % New sin(theta)
+        sinth = sinths(ii_b) ./ drel;
+
+        % Predict omega angles: 4 for each plane normal
+        om = gtFedPredictOmegaMultiple(pl_samd, sinth, labgeo.beamdir', ...
+            labgeo.rotdir', labgeo.rotcomp, ominds(ii_b));
+
+        valid_voxels(:, ii_b) = ~isnan(om);
+
+        % Delete those where no reflection occurs
+        if (any(isnan(om)))
+            inds_bad = find(isnan(om));
+            gvd(:, inds_bad(1:min(10, numel(inds_bad))))
+            warning('gtFedFwdProjExact:bad_R_vectors', ...
+                'No diffraction from some elements after deformation (%d over %d) for blob %d.', ...
+                numel(inds_bad), numel(om), ii_b)
+        end
+
+        w_bl = om ./ om_step;
+        w_bl = gtGrainAnglesTabularFix360deg(w_bl, ref_ws(ii_b), parameters);
+        w_bl = w_bl - w_shifts(ii_b);
+
+        w{ii_b} = w_bl;
+
+        w_min(ii_b, :) = min(w{ii_b});
+        w_max(ii_b, :) = max(w{ii_b});
+
+        o.fprintf(repmat('\b', [1, num_chars]));
+    end
+    o.fprintf('Done in %g s\n', toc(t));
+
+    o.fprintf('   * Computing max BBox size and feasibilities:\n')
+
+    if (linear_interp)
+        blob_w_sizes = ceil(w_max) - floor(w_min) + 1;
+        blob_orig_w_shifts = 1 - floor(w_min);
+    else
+        blob_w_sizes = round(w_max) - round(w_min) + 1;
+        blob_orig_w_shifts = 1 - round(w_min);
+    end
+
+    o.fprintf('     Computed Blob W sizes: [%s]\n', sprintf(' %d', blob_w_sizes));
+
+    bl = gtFwdSimBlobDefinition('sf_w', nbl);
+
+    %%% Blob projection
+    o.fprintf('   * Projecting volumes: ')
+    t = tic();
+    for ii_b = 1:nbl
+        num_chars = o.fprintf('%03d/%03d', ii_b, nbl);
+
+        % Detector coordinates U,V in blob
+        w_bl = w{ii_b} + blob_orig_w_shifts(ii_b * ones(1, nv))';
+
+        % Let's now filter valid voxels
+        w_bl = w_bl(:, valid_voxels(:, ii_b))';
+        gvpow_v = reshape(gvpow(valid_voxels(:, ii_b)), [], 1);
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Sum intensities
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Distribute value over 8 pixels
+
+        if (linear_interp)
+            [inds_w, ints_w] = gtMathsGetInterpolationIndices(w_bl, ...
+                gvpow_v, fedpars.bltype);
+        else
+            inds_w = round(w_bl);
+            ints_w = gvpow_v;
+        end
+
+        % Accumulate all intensities in the blob voxel-wise
+        %   'uvw' needs to be nx1; 'ints' is now 1xnx8
+        try
+            bl(ii_b).intm = accumarray(inds_w, ints_w, [blob_w_sizes(ii_b), 1]);
+        catch mexc
+            fprintf(['\n\nERROR\nThe error is probably caused by the' ...
+                ' projected intensities falling outside the blob volume.' ...
+                '\nTry increasing the blob volume padding:' ...
+                ' fedpars.detector(det_ind).blobsizeadd, or FIXING' ...
+                ' YOUR CODE!!\n\n'])
+            disp('Blob ID:')
+            disp(ii_b)
+            disp('Blob size:')
+            disp(blob_w_sizes(ii_b))
+            disp('Min projected W coordinate:')
+            fprintf('\t%g \t(rounded: %d )\n', ...
+                min(w_bl, [], 1), min(inds_w, [], 1))
+            disp('Max projected W coordinate:')
+            fprintf('\t%g \t(rounded: %d )\n', ...
+                max(w_bl, [], 1), max(inds_w, [], 1))
+            new_exc = GtFedExceptionFwdProj(mexc, det_ind, [0 0 0]);
+            throw(new_exc)
+        end
+
+        bl(ii_b).bbsize = blob_w_sizes(ii_b);
+
+        im_low_lim = w_shifts(ii_b) - blob_orig_w_shifts(ii_b) + 1;
+        bl(ii_b).bbwim = [im_low_lim, im_low_lim + bl(ii_b).bbsize - 1];
+
+        o.fprintf(repmat('\b', [1, num_chars]));
+    end
+    o.fprintf('Done in %g s\n', toc(t));
+end

zUtil_Deformation/gtDefShapeFunctionsFwdProjUVW.m → zUtil_Deformation/gtDefShapeFunctionsFwdProj.m

-function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, varargin)
+function shape_funcs = gtDefShapeFunctionsFwdProj(sampler, varargin)
 
     conf = struct( ...
+        'shape_function_type', 'uvw', ...
         'recenter_sf', true, ...
         'data_type', 'double', ...
-        'factor', 1, ...
+        'factor', 5, ...
         'rspace_oversampling', 1, ...
         'rspace_voxel_size', [1 1 1], ...
-        'interpolated_voxel', false);
+        'interpolated_voxel', false, ...
+        'projector', 'nearest');
     conf = parse_pv_pairs(conf, varargin);
 
     num_detectors = numel(sampler.parameters.detgeo);
@@ -22,15 +24,20 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, varargin)
         'bltype', conf.data_type, ...
         'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
         'defmethod', 'rod_rightstretch', ...
-        'detector', repmat( fed_pars_detector, [num_detectors, 1]));
+        'detector', repmat( fed_pars_detector, [num_detectors, 1]), ...
+        'projector', conf.projector );
 
-    voxel_factor = conf.factor([1 1 1]);
-    num_sub_voxels = prod(voxel_factor);
+    voxel_size = sampler.stats.sampling.gaps;
+
+    space_res = tand(sampler.estimate_maximum_resolution() ./ 2);
+    voxel_sampling = max(ceil(voxel_size ./ space_res * conf.factor), 5);
+
+    num_sub_voxels = prod(voxel_sampling);
     ones_nsv = ones(num_sub_voxels, 1);
 
     gv = struct('pcs', [], 'ind', [], 'd', [], 'pow', []);
 
-    [gv.ind(:, 1), gv.ind(:, 2), gv.ind(:, 3)] = ind2sub(voxel_factor, 1:num_sub_voxels);
+    [gv.ind(:, 1), gv.ind(:, 2), gv.ind(:, 3)] = ind2sub(voxel_sampling, 1:num_sub_voxels);
     gv.ind = gv.ind';
 
     ref_gr = sampler.get_reference_grain();
@@ -40,10 +47,10 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, varargin)
     voxel_size = sampler.stats.sampling.gaps;
     if (conf.interpolated_voxel)
         half_voxel_size = voxel_size; % Enlarging to neighbours
-        steps = voxel_size ./ voxel_factor * 2;
+        steps = voxel_size ./ voxel_sampling * 2;
     else
         half_voxel_size = voxel_size ./ 2;
-        steps = voxel_size ./ voxel_factor;
+        steps = voxel_size ./ voxel_sampling;
     end
     half_steps = steps ./ 2;
     beg_pos = - half_voxel_size + half_steps;
@@ -83,7 +90,7 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, varargin)
     num_ors = numel(sampler.lattice(1).gr);
     shape_funcs = cell(num_ors, 1);
 
-    fprintf('Computing UVW shape functions: ')
+    fprintf('Computing %s shape functions: ', upper(conf.shape_function_type))
     c = tic();
     for ii_g = 1:num_ors
         num_chars = fprintf('%03d/%03d', ii_g, num_ors);
@@ -91,10 +98,20 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, varargin)
         or = sampler.lattice(1).gr{ii_g};
 
         or_gv = gv;
-        or_gv.d = gtMathsRodSum(or.R_vector', or_gv.d);
-
-        shape_funcs{ii_g} = gtDefFwdProjGvdm(or, sel_bls, or_gv, ...
-            fedpars, sampler.parameters, det_ind, false);
+        or_gv.d = gtMathsRodSum(gv.d, or.R_vector');
+
+        switch (lower(conf.shape_function_type))
+            case 'w'
+                shape_funcs{ii_g} = gtDefFwdProjGvdm2W(or, sel_bls, or_gv, ...
+                    fedpars, sampler.parameters, det_ind, false);
+            case 'nw'
+                eta_step = 2 * atand(space_res(1));
+                shape_funcs{ii_g} = gtDefFwdProjGvdm2NW(or, sel_bls, or_gv, ...
+                    fedpars, sampler.parameters, eta_step, det_ind, false);
+            case 'uvw'
+                shape_funcs{ii_g} = gtDefFwdProjGvdm(or, sel_bls, or_gv, ...
+                    fedpars, sampler.parameters, det_ind, false);
+        end
 
         fprintf(repmat('\b', [1, num_chars]));
     end