From b96820c0d69ba28eb0c1424926564c5c28d137ad Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Mon, 22 Feb 2016 17:03:57 +0100
Subject: [PATCH] GtOrientationSampling: added orientation-space shape
functions generation for the sampled grid
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Deformation/GtOrientationSampling.m | 220 ++++++++++++++++++++++
1 file changed, 220 insertions(+)
diff --git a/zUtil_Deformation/GtOrientationSampling.m b/zUtil_Deformation/GtOrientationSampling.m
index b7dbb170..fa0de10c 100644
--- a/zUtil_Deformation/GtOrientationSampling.m
+++ b/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));
--
GitLab