diff --git a/zUtil_Deformation/GtOrientationSampling.m b/zUtil_Deformation/GtOrientationSampling.m
index 710671ebd26cd4cd4cc5ba78eb392a0093959797..6145aaad22df06e9a8647ad244879121dd070f38 100644
--- a/zUtil_Deformation/GtOrientationSampling.m
+++ b/zUtil_Deformation/GtOrientationSampling.m
@@ -510,6 +510,7 @@ classdef GtOrientationSampling < handle
%
omstep = gtGetOmegaStepDeg(self.parameters, self.detector_index);
+ omstep_rod = tand(omstep / 2);
sub_space_size = self.stats.sampling.gaps;
@@ -566,6 +567,28 @@ classdef GtOrientationSampling < handle
ref_ws = self.ref_gr.allblobs.omega(self.ondet(self.included(ref_inds))) ./ omstep;
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 = 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);
+
+ 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));
+
+ tot_sampling_points = num_poins(1) * num_poins(2);
+ ones_tot_sp = ones(tot_sampling_points, 1);
+
% 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();
@@ -575,57 +598,30 @@ classdef GtOrientationSampling < handle
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;
+ ref_round_int_ws = round(ref_ws + 0.5) - 0.5;
+ ref_round_ws = ref_round_int_ws .* omstep;
+ rot_w = gtMathsRotationTensor(ref_round_ws, rotcomp_w);
- rot_w = gtMathsRotationTensor(ws_lims, rotcomp_w);
+ beamdirs = self.parameters.labgeo.beamdir * reshape(rot_w, 3, []);
+ beamdirs = reshape(beamdirs, 3, [])';
- 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;
- Ac_part = beamdirs * rotcomp_z.cos;
- As_part = beamdirs * rotcomp_z.sin;
- Acc_part = beamdirs * rotcomp_z.const;
+ Ac_part = Ac_part(:, :, ones_tot_sp);
+ As_part = As_part(:, :, ones_tot_sp);
+ Acc_part = Acc_part(:, :, ones_tot_sp);
- 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);
@@ -637,8 +633,6 @@ classdef GtOrientationSampling < handle
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
@@ -647,61 +641,75 @@ classdef GtOrientationSampling < handle
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;
+ r_vecs = gr.R_vector(ones_tot_sp, :)' + 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);
+ ys = gs * 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 .* ys, 2);
+ As = sum(As_part .* ys, 2);
+ Acc = sum(Acc_part .* ys, 2);
+ Ac = reshape(Ac, num_blobs, tot_sampling_points);
+ As = reshape(As, num_blobs, tot_sampling_points);
+ Acc = reshape(Acc, num_blobs, tot_sampling_points);
+
+ D = Ac .^ 2 + As .^ 2;
+
+ % Precomputing useful values, like sinthetas
ominds = gr.allblobs.omind(ref_inds);
ssp = ((ominds == 1) | (ominds == 2));
ss = ssp - ~ssp;
sinths = ss .* gr.allblobs.sintheta(ref_inds);
+ sinths = sinths(:, ones(tot_sampling_points, 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));
ss = ssp - ~ssp;
+ ss = ss(:, ones(tot_sampling_points, 1));
- 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);
+ % Shift along z in orientation space, to align with the images
+ % in proections space
+ dz = (-As + ss .* ok .* E) ./ (CC - Ac);
- Ac = Ac_g_part * pls_blob;
- As = As_g_part * pls_blob;
- Acc = Acc_g_part * pls_blob;
+ 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 .* omstep / 2);
- D = Ac .^ 2 + As .^ 2;
+ num_ds = lims_w_proj_g_int(:, 2) - lims_w_proj_g_int(:, 1) + 2;
- CC = Acc + sinths(ii_b);
- DD = D - CC .^ 2;
- E = sqrt(DD);
- ok = DD > 0;
+ bbwims = cell(num_blobs, 1);
+ intms = cell(num_blobs, 1);
+ bbsizes = cell(num_blobs, 1);
- % 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);
+ % 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
+ d = (ws_disps(ii_b, 1):omstep_rod:ws_disps(ii_b, 2))';
+ d = d(:, ones(tot_sampling_points, 1)) + dz(ii_b * ones(num_ds(ii_b), 1), :);
upper_limits = min(d(2:end, :), half_rspace_sizes(3));
lower_limits = max(d(1:end-1, :), -half_rspace_sizes(3));
@@ -714,10 +722,17 @@ classdef GtOrientationSampling < handle
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)];
+ bbwims{ii_b} = lims_w_proj_g_int(ii_b, :);
+ intms{ii_b} = renorm_coeff * ints;
+ bbsizes{ii_b} = [1, 1, numel(ints)];
end
+
+ sf_bls_w(1:num_blobs) = gtFwdSimBlobDefinition();
+ [sf_bls_w(:).bbwim] = deal(bbwims{:});
+ [sf_bls_w(:).intm] = deal(intms{:});
+ [sf_bls_w(:).bbsize] = deal(bbsizes{:});
+ sf{ii_g} = sf_bls_w;
+
if (mod(ii_g, chunk_size) == 1)
fprintf(repmat('\b', [1, num_chars]));
end