GtOrientationSampling: added orientation-space shape functions generation for the sampled grid

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

zUtil_Deformation/GtOrientationSampling.m

@@ -505,6 +505,226 @@ classdef GtOrientationSampling < handle
             end
         end
 
+        function sf = compute_shape_functions_w(self)
+        % FUNCTION sf = compute_shape_functions_w(self)
+        %
+
+            omstep = gtGetOmegaStepDeg(self.parameters, self.detector_index);
+
+            sub_space_size = self.stats.sampling.gaps;
+
+            % Let's first compute the W extent
+            half_rspace_sizes = sub_space_size / 2;
+
+            ref_inds = self.selected;
+            num_blobs = numel(ref_inds);
+
+            pl_samd = self.ref_gr.allblobs.plorig(self.ondet(self.included(ref_inds)), :);
+
+            g = gtMathsRod2OriMat(self.ref_gr.R_vector');
+            pl_cry = gtVectorLab2Cryst(pl_samd, g)';
+%             pl_labd = gtGeoSam2Lab(pl_samd, eye(3), self.parameters.labgeo, self.parameters.samgeo, true, false)';
+
+            fprintf('Computing bounds of W shape functions: ')
+            c = tic();
+            tot_orientations = numel(self.lattice.gr);
+
+            oris_bounds_size = [size(self.lattice.gr, 1), size(self.lattice.gr, 2), size(self.lattice.gr, 3)] + 1;
+            oris_lims_min = self.lattice.gr{1, 1, 1}.R_vector - half_rspace_sizes;
+            oris_lims_max = self.lattice.gr{end, end, end}.R_vector + half_rspace_sizes;
+
+            oris_steps_x = linspace(oris_lims_min(1), oris_lims_max(1), oris_bounds_size(1));
+            oris_steps_y = linspace(oris_lims_min(2), oris_lims_max(2), oris_bounds_size(2));
+            oris_steps_z = linspace(oris_lims_min(3), oris_lims_max(3), oris_bounds_size(3));
+
+            oris_bounds = cell(oris_bounds_size);
+            for dx = 1:oris_bounds_size(1)
+                for dy = 1:oris_bounds_size(2)
+                    for dz = 1:oris_bounds_size(3)
+                        r_vec = [oris_steps_x(dx), oris_steps_y(dy), oris_steps_z(dz)];
+                        oris_bounds{dx, dy, dz} = struct( ...
+                            'phaseid', self.ref_gr.phaseid, ...
+                            'center', self.ref_gr.center, 'R_vector', r_vec );
+                    end
+                end
+            end
+
+            chunk_size = 250;
+            for ii = 1:chunk_size:numel(oris_bounds)
+                end_chunk = min(ii+chunk_size, numel(oris_bounds));
+                num_chars_gr_ii = fprintf('%05d/%05d', ii, numel(oris_bounds));
+
+                oris_bounds(ii:end_chunk) = gtCalculateGrain_p( ...
+                    oris_bounds(ii:end_chunk), self.parameters, ...
+                    'ref_omind', self.ref_gr.allblobs.omind, ...
+                    'included', self.ondet(self.included) );
+
+                fprintf(repmat('\b', [1 num_chars_gr_ii]));
+            end
+            fprintf('Done in %g seconds.\n', toc(c))
+
+            ref_ws = self.ref_gr.allblobs.omega(self.ondet(self.included(ref_inds))) ./ omstep;
+
+            fprintf('Precomputing shared values..')
+            % There is no point in computing the same rotation matrices each
+            % time, so we pre compute them here and store them for later
+            c = tic();
+            % selecting Ws to find the interval
+            ori_bounds = [oris_bounds{:}];
+            allblobs = [ori_bounds(:).allblobs];
+            w_tabs_int = [allblobs(:).omega];
+            w_tabs_int = w_tabs_int(ref_inds, :) ./ omstep;
+            w_tabs_int = gtGrainAnglesTabularFix360deg(w_tabs_int, ref_ws, self.parameters);
+
+            lims_w_proj_int = [min(w_tabs_int, [], 2), max(w_tabs_int, [], 2)];
+            lims_w_proj_int = [floor(lims_w_proj_int(:, 1)), ceil(lims_w_proj_int(:, 2))];
+
+            w_tabs_int = reshape(w_tabs_int, [], oris_bounds_size(1), oris_bounds_size(2), oris_bounds_size(3));
+
+            precomp_matrix_prods = cell(num_blobs, 1);
+
+            rotcomp_z = gtMathsRotationMatrixComp([0, 0, 1]', 'col');
+
+            rotcomp_w = gtMathsRotationMatrixComp(self.parameters.labgeo.rotdir', 'col');
+
+            % Precomputing A matrix components for all the used omegas
+            for ii_b = 1:num_blobs
+                ws_lims = (lims_w_proj_int(ii_b, 1)-0.5):(lims_w_proj_int(ii_b, 2)+0.5);
+                ws_lims = ws_lims .* omstep;
+
+                rot_w = gtMathsRotationTensor(ws_lims, rotcomp_w);
+
+                beamdirs = self.parameters.labgeo.beamdir * reshape(rot_w, 3, []);
+                beamdirs = reshape(beamdirs, 3, [])';
+
+                Ac_part = beamdirs * rotcomp_z.cos;
+                As_part = beamdirs * rotcomp_z.sin;
+                Acc_part = beamdirs * rotcomp_z.const;
+
+                precomp_matrix_prods{ii_b} = struct( ...
+                    'w', [lims_w_proj_int(ii_b, 1) lims_w_proj_int(ii_b, 2)], ...
+                    'Ac', Ac_part, 'As', As_part, 'Acc', Acc_part );
+            end
+            fprintf('\b\b: Done in %g seconds.\n', toc(c))
+
+            space_res = self.estimate_maximum_resolution() ./ 5;
+            num_poins = max(ceil(sub_space_size ./ space_res), 5);
+            pos_x = linspace(-half_rspace_sizes(1), half_rspace_sizes(1), num_poins(1)+1);
+            pos_y = linspace(-half_rspace_sizes(2), half_rspace_sizes(2), num_poins(2)+1);
+
+            renorm_coeff = (pos_x(2) - pos_x(1)) * (pos_y(2) - pos_y(1)) / prod(sub_space_size);
+
+            % Forming the grid in orientation space, on the plane that is
+            % perpendicular to the axis of rotation (z-axis in the w case),
+            % where to sample the sub-regions of the orientation-space voxels
+            pos_x = (pos_x(1:end-1) + pos_x(2:end)) / 2;
+            pos_y = (pos_y(1:end-1) + pos_y(2:end)) / 2;
+
+            pos_x_exp = reshape(pos_x, 1, [], 1);
+            pos_x_exp = pos_x_exp(1, :, ones(1, num_poins(2)));
+            pos_y_exp = reshape(pos_y, 1, 1, []);
+            pos_y_exp = pos_y_exp(1, ones(1, num_poins(1)), :);
+            pos_z_exp = zeros(1, num_poins(1), num_poins(2));
+
+            sf = cell(tot_orientations, 1);
+
+            [o_x, o_y, o_z] = ind2sub(oris_bounds_size-1, 1:tot_orientations);
+
+            fprintf('Computing W shape functions: ')
+            c = tic();
+            for ii_g = 1:tot_orientations
+                if (mod(ii_g, chunk_size) == 1)
+                    num_chars = fprintf('%05d/%05d', ii_g, tot_orientations);
+                end
+
+                sf{ii_g}(1:num_blobs) = gtFwdSimBlobDefinition();
+
+                gr = self.lattice.gr{ii_g};
+
+                % Extracting ospace boundaries for the given voxel
+                w_tab_int = w_tabs_int(:, o_x(ii_g):o_x(ii_g)+1, o_y(ii_g):o_y(ii_g)+1, o_z(ii_g):o_z(ii_g)+1);
+                w_tab_int = reshape(w_tab_int, [], 8);
+
+                lims_w_proj_g_int = round([min(w_tab_int, [], 2), max(w_tab_int, [], 2)]);
+
+                % the extra index is needed for the last boundary point
+                As_inds = [ ...
+                    lims_w_proj_g_int(:, 1) - lims_w_proj_int(:, 1) ...
+                    lims_w_proj_g_int(:, 2) - lims_w_proj_int(:, 1) + 1] + 1;
+
+                % Forming the grid in orientation space, on the plane that is
+                % perpendicular to the axis of rotation (z-axis in the w case),
+                % where to sample the sub-regions of the orientation-space voxel
+                % that we want to integrate.
+                r_vecs = [pos_x_exp; pos_y_exp; pos_z_exp];
+                r_vecs = reshape(r_vecs, 3, []);
+                r_vecs = gr.R_vector(ones(1, size(r_vecs, 2)), :)' + r_vecs;
+
+                gs = gtMathsRod2OriMat(r_vecs);
+                % Unrolling and transposing the matrices at the same time
+                % because we need g^-1
+                gs = reshape(gs, 3, [])';
+
+                pls = gs * pl_cry;
+                pls = reshape(pls, 3, [], num_blobs);
+
+                ominds = gr.allblobs.omind(ref_inds);
+                ssp = ((ominds == 1) | (ominds == 2));
+                ss  = ssp - ~ssp;
+                sinths = ss .* gr.allblobs.sintheta(ref_inds);
+
+                % Apparently it is inverted, compared to the omega prediction
+                % function. Should be investigated
+                ssp = ~((ominds == 1) | (ominds == 3));
+                ss  = ssp - ~ssp;
+
+                for ii_b = 1:num_blobs
+                    A_inds = As_inds(ii_b, 1):As_inds(ii_b, 2);
+
+                    Ac_g_part = precomp_matrix_prods{ii_b}.Ac(A_inds, :);
+                    As_g_part = precomp_matrix_prods{ii_b}.As(A_inds, :);
+                    Acc_g_part = precomp_matrix_prods{ii_b}.Acc(A_inds, :);
+
+                    pls_blob = pls(:, :, ii_b);
+
+                    Ac = Ac_g_part * pls_blob;
+                    As = As_g_part * pls_blob;
+                    Acc = Acc_g_part * pls_blob;
+
+                    D = Ac .^ 2 + As .^ 2;
+
+                    CC = Acc + sinths(ii_b);
+                    DD = D - CC .^ 2;
+                    E  = sqrt(DD);
+                    ok = DD > 0;
+
+                    % finding what the positions in w-space, correspond to in
+                    % orientation-space, on the z-axis, at the sampled x and y
+                    % positions around in the considered volume.
+                    d = (-As + ss(ii_b) .* ok .* E) ./ (CC - Ac);
+
+                    upper_limits = min(d(2:end, :), half_rspace_sizes(3));
+                    lower_limits = max(d(1:end-1, :), -half_rspace_sizes(3));
+
+                    % integration is simple and linear: for every w and
+                    % xy-position, we consider the corresponding segment on the
+                    % z-direction
+                    ints = upper_limits - lower_limits;
+                    ints(ints < 0) = 0;
+                    ints = sum(ints, 2);
+                    ints = reshape(ints, 1, 1, []);
+
+                    sf{ii_g}(ii_b).bbwim = lims_w_proj_g_int(ii_b, :);
+                    sf{ii_g}(ii_b).intm = renorm_coeff * ints;
+                    sf{ii_g}(ii_b).bbsize = [1, 1, numel(ints)];
+                end
+                if (mod(ii_g, chunk_size) == 1)
+                    fprintf(repmat('\b', [1, num_chars]));
+                end
+            end
+            fprintf('Done in %g seconds.\n', toc(c))
+        end
+
         function resolution = estimate_maximum_resolution(self)
             sel = self.ondet(self.included(self.selected));