diff --git a/zUtil_Deformation/GtOrientationSampling.m b/zUtil_Deformation/GtOrientationSampling.m
index 7eb655c03c6a18a8ad7bf20835f0425c71e5ede0..467c03734d506887c26072374e65fb464cc825cc 100644
--- a/zUtil_Deformation/GtOrientationSampling.m
+++ b/zUtil_Deformation/GtOrientationSampling.m
@@ -739,6 +739,405 @@ classdef GtOrientationSampling < handle
fprintf('Done in %g seconds.\n', toc(c))
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)
sel = self.ondet(self.included(self.selected));