Skip to content
Snippets Groups Projects
Commit e060a5b4 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

OrientationSampling: added function for computing the nw shape functions 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent f2d4ce0a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/GtOrientationSampling.m
+ 399
− 0
View file @ e060a5b4
...@@ -739,6 +739,405 @@ classdef GtOrientationSampling < handle ...@@ -739,6 +739,405 @@ classdef GtOrientationSampling < handle
fprintf('Done in %g seconds.\n', toc(c)) fprintf('Done in %g seconds.\n', toc(c))
end end
function sf = compute_shape_functions_nw(self, compute_w_sfs)
% FUNCTION sf = compute_shape_functions_nw(self)
%
if (~exist('compute_w_sfs', 'var') || isempty(compute_w_sfs))
compute_w_sfs = false;
end
om_step = gtGetOmegaStepDeg(self.parameters, self.detector_index);
om_step_rod = tand(om_step / 2);
oversampling = 11;
step = 1 / oversampling;
shift = 0.5 - step / 2;
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 NW 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))
fprintf('Precomputing shared values..')
% 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
space_res = tand(self.estimate_maximum_resolution() ./ 2);
if (compute_w_sfs)
num_poins = max(ceil(sub_space_size ./ space_res ./ 2), 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_w = (pos_x(2) - pos_x(1)) * (pos_y(2) - pos_y(1)) / prod(sub_space_size);
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_grid_w = reshape(pos_x, 1, [], 1);
pos_x_grid_w = pos_x_grid_w(1, :, ones(1, num_poins(2)));
pos_y_grid_w = reshape(pos_y, 1, 1, []);
pos_y_grid_w = pos_y_grid_w(1, ones(1, num_poins(1)), :);
pos_z_grid_w = zeros(1, num_poins(1), num_poins(2));
tot_sampling_points_w = num_poins(1) * num_poins(2);
ones_tot_sp_w = ones(tot_sampling_points_w, 1);
end
ref_ws = self.ref_gr.allblobs.omega(self.ondet(self.included(ref_inds))) ./ om_step;
% 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, :) ./ om_step;
w_tabs_int = gtGrainAnglesTabularFix360deg(w_tabs_int, ref_ws, self.parameters);
w_tabs_int = reshape(w_tabs_int, [], oris_bounds_size(1), oris_bounds_size(2), oris_bounds_size(3));
eta_step = 2 * atand(space_res(1));
ref_ns = self.ref_gr.allblobs.eta(self.ondet(self.included(ref_inds)));
n_tabs_int = [allblobs(:).eta];
n_tabs_int = n_tabs_int(ref_inds, :);
n_tabs_int = gtGrainAnglesTabularFix360deg(n_tabs_int, ref_ns);
n_tabs_int = n_tabs_int ./ eta_step;
n_tabs_int = reshape(n_tabs_int, [], oris_bounds_size(1), oris_bounds_size(2), oris_bounds_size(3));
rotcomp_w = gtMathsRotationMatrixComp(self.parameters.labgeo.rotdir', 'col');
if (compute_w_sfs)
rotcomp_z = gtMathsRotationMatrixComp([0, 0, 1]', 'col');
% Precomputing A matrix components for all the used omegas
ref_round_int_ws = round(ref_ws + 0.5) - 0.5;
ref_round_ws = ref_round_int_ws .* om_step;
rot_samp_w = gtMathsRotationTensor(ref_round_ws, rotcomp_w);
beamdirs = self.parameters.labgeo.beamdir * reshape(rot_samp_w, 3, []);
beamdirs = reshape(beamdirs, 3, [])';
Ac_part_w = beamdirs * rotcomp_z.cos;
As_part_w = beamdirs * rotcomp_z.sin;
Acc_part_w = beamdirs * rotcomp_z.const;
Ac_part_w = Ac_part_w(:, :, ones_tot_sp_w);
As_part_w = As_part_w(:, :, ones_tot_sp_w);
Acc_part_w = Acc_part_w(:, :, ones_tot_sp_w);
end
fprintf('\b\b: Done in %g seconds.\n', toc(c))
sf = cell(tot_orientations, 1);
[o_x, o_y, o_z] = ind2sub(oris_bounds_size-1, 1:tot_orientations);
fprintf('Computing NW shape functions: ')
c = tic();
for ii_g = 1:tot_orientations
if (mod(ii_g, chunk_size/10) == 1)
currtime = toc(c);
num_chars = fprintf('%05d/%05d (%g s, ETA %g s)', ...
ii_g, tot_orientations, currtime, ...
toc(c) / (ii_g - 1) * tot_orientations - currtime);
end
gr = self.lattice.gr{ii_g};
t = gr.allblobs.theta(ref_inds);
cos_thetas = cosd(t);
tan_thetas = tand(t);
n = gr.allblobs.eta(ref_inds);
cos_etas = cosd(-n);
sin_etas = sind(-n);
w = gr.allblobs.omega(ref_inds);
rot_w = gtMathsRotationTensor(w, rotcomp_w);
n_1 = gr.allblobs.eta(ref_inds) + eta_step;
cos_etas_1 = cosd(-n_1);
sin_etas_1 = sind(-n_1);
rot_w_1 = gtMathsRotationTensor(w + om_step, rotcomp_w);
% 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);
n_tab_int = n_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);
n_tab_int = reshape(n_tab_int, [], 8);
lims_w_proj_g_int = round([min(w_tab_int, [], 2), max(w_tab_int, [], 2)]);
lims_n_proj_g_int = round([min(n_tab_int, [], 2), max(n_tab_int, [], 2)]);
num_ws = lims_w_proj_g_int(:, 2) - lims_w_proj_g_int(:, 1) + 1;
num_ns = lims_n_proj_g_int(:, 2) - lims_n_proj_g_int(:, 1) + 1;
% Precomputing useful values, like sinthetas
ominds = gr.allblobs.omind(ref_inds);
ssp = ((ominds == 1) | (ominds == 2));
ss12 = ssp - ~ssp;
[ref_r0, ref_r_dir] = compute_rod_line(cos_thetas, ...
tan_thetas, sin_etas, cos_etas, rot_w, pl_cry, ss12);
[ref_r0_n, ref_r_dir_n] = compute_rod_line(cos_thetas, ...
tan_thetas, sin_etas_1, cos_etas_1, rot_w, pl_cry, ss12);
[ref_r0_w, ref_r_dir_w] = compute_rod_line(cos_thetas, ...
tan_thetas, sin_etas, cos_etas, rot_w_1, pl_cry, ss12);
ref_r0_diff_n = ref_r0_n - ref_r0;
ref_r0_diff_w = ref_r0_w - ref_r0;
ref_r_dir_diff_n = ref_r_dir_n - ref_r_dir;
ref_r_dir_diff_w = ref_r_dir_w - ref_r_dir;
ref_inv_r_dir = 1 ./ ref_r_dir;
% Translating everything to the origin
ref_r0 = ref_r0 - gr.R_vector(ones(numel(t), 1), :);
if (compute_w_sfs)
num_dws = num_ws + 1;
num_dns = num_ns + 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_w_samp = [pos_x_grid_w; pos_y_grid_w; pos_z_grid_w];
r_vecs_w_samp = reshape(r_vecs_w_samp, 3, []);
r_vecs_w_samp = gr.R_vector(ones_tot_sp_w, :)' + r_vecs_w_samp;
gs_w = gtMathsRod2OriMat(r_vecs_w_samp);
% Unrolling and transposing the matrices at the same time
% because we need g^-1
gs_w = reshape(gs_w, 3, [])';
ys = gs_w * pl_cry;
ys = reshape(ys, 3, [], num_blobs);
ys = permute(ys, [3 1 2]);
% Starting here to compute initial shifts aligned with the
% rounded omega, by frst doing three scalar products
Ac = sum(Ac_part_w .* ys, 2);
As = sum(As_part_w .* ys, 2);
Acc = sum(Acc_part_w .* ys, 2);
Ac = reshape(Ac, num_blobs, tot_sampling_points_w);
As = reshape(As, num_blobs, tot_sampling_points_w);
Acc = reshape(Acc, num_blobs, tot_sampling_points_w);
D = Ac .^ 2 + As .^ 2;
sinths = ss12 .* gr.allblobs.sintheta(ref_inds);
sinths = sinths(:, ones(tot_sampling_points_w, 1));
CC = Acc + sinths;
DD = D - CC .^ 2;
E = sqrt(DD);
ok = DD > 0;
% Apparently it is inverted, compared to the omega prediction
% function. Should be investigated
ssp = ~((ominds == 1) | (ominds == 3));
ss13 = ssp - ~ssp;
ss13 = ss13(:, ones(tot_sampling_points_w, 1));
% Shift along z in orientation space, to align with the images
% in proections space
dz = (-As + ss13 .* ok .* E) ./ (CC - Ac);
ws_disps = [ ...
(lims_w_proj_g_int(:, 1) - 0.5 - ref_round_int_ws), ...
(lims_w_proj_g_int(:, 2) + 0.5 - ref_round_int_ws) ];
ws_disps = tand(ws_disps .* om_step / 2);
end
sf_bls_nw(1:num_blobs) = gtFwdSimBlobDefinition('sf_nw');
% This is the most expensive operation because the different
% sizes of the omega spread, which depends on the blob itself,
% don't allow to have vectorized operations on all the blobs at
% the same time
for ii_b = 1:num_blobs
ws = ((lims_w_proj_g_int(ii_b, 1)-shift):step:(lims_w_proj_g_int(ii_b, 2)+shift)) .* om_step;
ones_ws = ones(oversampling * num_ws(ii_b), 1);
ns = ((lims_n_proj_g_int(ii_b, 1)-shift):step:(lims_n_proj_g_int(ii_b, 2)+shift)) .* eta_step;
ones_ns = ones(oversampling * num_ns(ii_b), 1);
% It is always in x, because this is still in lab
% coords!
tot_vecs = oversampling * num_ws(ii_b) * oversampling * num_ns(ii_b);
ones_tvecs = ones(tot_vecs, 1);
vars_ns = (ns - n(ii_b)) / eta_step;
vars_ns = vars_ns([1 1 1], :);
vars_ws = (ws - w(ii_b)) / om_step;
vars_ws = vars_ws([1 1 1], :);
% The system has already been translated to the origin
rep_r0s = ref_r0(ii_b, :)';
rep_r0s = rep_r0s(:, ones_ns, ones_ws);
vars_r0_ns = vars_ns .* ref_r0_diff_n(ii_b * ones_ns, :)';
vars_r0_ns = vars_r0_ns(:, :, ones_ws);
vars_r0_ws = vars_ws .* ref_r0_diff_w(ii_b * ones_ws, :)';
vars_r0_ws = reshape(vars_r0_ws, 3, 1, []);
vars_r0_ws = vars_r0_ws(:, ones_ns, :);
% r0 should be defined = (h x y) / (1 + h . y)
approx_r0 = rep_r0s + vars_r0_ns + vars_r0_ws;
r0 = reshape(approx_r0, 3, [])';
rep_r_dirs = ref_r_dir(ii_b * ones_tvecs, :);
vars_r_dir_ns = vars_ns .* ref_r_dir_diff_n(ii_b * ones_ns, :)';
vars_r_dir_ns = vars_r_dir_ns(:, :, ones_ws);
vars_r_dir_ns = reshape(vars_r_dir_ns, 3, [])';
vars_r_dir_ws = vars_ws .* ref_r_dir_diff_w(ii_b * ones_ws, :)';
vars_r_dir_ws = reshape(vars_r_dir_ws, 3, 1, []);
vars_r_dir_ws = vars_r_dir_ws(:, ones_ns, :);
vars_r_dir_ws = reshape(vars_r_dir_ws, 3, [])';
% r should be defined = r0 + t * (h + y) / (1 + h . y)
r_dir = rep_r_dirs + vars_r_dir_ns + vars_r_dir_ws;
% Translating r0s to the closest point to the average R
% vector, on the line r
dot_r0_avgR = sum((0 - r0) .* r_dir, 2);
r0 = r0 + r_dir .* dot_r0_avgR(:, [1 1 1]);
% Let's now find directions and limits
minus_dirs = r_dir < 0;
u_lims = half_rspace_sizes(ones_tvecs, :);
u_lims(minus_dirs) = -u_lims(minus_dirs);
l_lims = -u_lims;
% Approximate inverse: taylor expansion to the 1st
% order
inv_r_dir = ref_inv_r_dir(ii_b * ones_tvecs, :);
inv_r_dir = (1 - (vars_r_dir_ns + vars_r_dir_ws) .* inv_r_dir) .* inv_r_dir;
len_to_ulims = (u_lims - r0) .* inv_r_dir;
len_to_llims = (l_lims - r0) .* inv_r_dir;
len_to_ulims = min(len_to_ulims, [], 2);
len_to_llims = max(len_to_llims, [], 2);
n_ints = len_to_ulims - len_to_llims;
n_ints(n_ints < 0) = 0;
n_ints = reshape(n_ints, [oversampling, num_ns(ii_b), oversampling, num_ws(ii_b)]);
n_ints = sum(sum(n_ints, 1), 3);
n_ints = reshape(n_ints, [num_ns(ii_b), num_ws(ii_b)]);
n_ints = n_ints ./ sum(n_ints(:));
if (compute_w_sfs)
d = (ws_disps(ii_b, 1):om_step_rod:ws_disps(ii_b, 2))';
d = d(:, ones(tot_sampling_points_w, 1)) + dz(ii_b * ones(num_dws(ii_b), 1), :);
w_upper_limits = min(d(2:end, :), half_rspace_sizes(3));
w_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
w_ints = w_upper_limits - w_lower_limits;
w_ints(w_ints < 0) = 0;
w_ints = renorm_coeff_w * sum(w_ints, 2)';
[sum(n_ints, 1), w_ints]
pause
end
sf_bls_nw(ii_b).bbnim = lims_n_proj_g_int(ii_b, :);
sf_bls_nw(ii_b).bbwim = lims_w_proj_g_int(ii_b, :);
sf_bls_nw(ii_b).intm = n_ints;
sf_bls_nw(ii_b).bbsize = [num_ns(ii_b) num_ws(ii_b)];
end
sf{ii_g} = sf_bls_nw;
if (mod(ii_g, chunk_size/10) == 1)
fprintf(repmat('\b', [1, num_chars]));
end
end
fprintf(repmat(' ', [1, num_chars]));
fprintf(repmat('\b', [1, num_chars]));
fprintf('Done in %g seconds.\n', toc(c))
function [r0, r_dir] = compute_rod_line(cos_thetas, tan_thetas, sin_etas, cos_etas, rot_w, pl_cry, ss12)
y = repmat(cos_thetas, [1 3]) .* [-tan_thetas, -sin_etas, cos_etas];
y = permute(y, [2 3 1]);
y = y(:, [1 1 1], :);
y = reshape(y, 3, []);
y_w = y .* reshape(rot_w, 3, []);
y_w = reshape(sum(y_w, 1), 3, [])';
h_12 = ss12(:, [1 1 1]) .* pl_cry';
renorms = 1 + sum(y_w .* h_12, 2);
renorms = renorms(:, [1 1 1]);
% r0 is defined = (h x y) / (1 + h . y)
r0 = gtMathsCross(h_12, y_w) ./ renorms;
% r is defined = r0 + t * (h + y) / (1 + h . y)
r_dir = (h_12 + y_w) ./ renorms;
r_dir = gtMathsNormalizeVectorsList(r_dir);
end
end
function resolution = estimate_maximum_resolution(self) function resolution = estimate_maximum_resolution(self)
sel = self.ondet(self.included(self.selected)); sel = self.ondet(self.included(self.selected));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment