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

ShapeFunctions/NW: fixed bug related to approximation 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 2109d0bc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/gtDefShapeFunctionsCreateNW.m
+ 70
− 61
View file @ 528b7a5d
......@@ -5,6 +5,7 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
conf = struct( ...
'data_type', 'double', ...
'fast_inverse', false, ...
'fast_lines', false, ...
'compute_w_sfs', false );
conf = parse_pv_pairs(conf, varargin);
......@@ -43,6 +44,7 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
g = gtMathsRod2OriMat(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: ')
......@@ -169,19 +171,19 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
gr = sampler.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);
ominds = gr.allblobs.omind(ref_inds);
ssp = ((ominds == 1) | (ominds == 2));
ss12 = ssp - ~ssp;
h = ss12(:, [1 1 1]) .* pl_cry';
% 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);
......@@ -195,17 +197,9 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
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, ref_r_dir] = compute_rod_line(t, n, rot_w, h, p.labgeo, p.samgeo);
[ref_r0_n, ref_r_dir_n] = compute_rod_line(t, n_1, rot_w, h, p.labgeo, p.samgeo);
[ref_r0_w, ref_r_dir_w] = compute_rod_line(t, n, rot_w_1, h, p.labgeo, p.samgeo);
ref_r0_diff_n = ref_r0_n - ref_r0;
ref_r0_diff_w = ref_r0_w - ref_r0;
......@@ -292,41 +286,63 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
tot_vecs = oversampling_w * num_ws(ii_b) * oversampling_n * 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], :);
if (conf.fast_lines)
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);
% 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_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, :);
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, [])';
% 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, :);
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_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, [])';
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, [])';
d_vers_nw = vars_r_dir_ns + vars_r_dir_ws;
d_vers_nw = vars_r_dir_ns + vars_r_dir_ws;
% r should be defined = r0 + t * (h + y) / (1 + h . y)
r_dir = rep_r_dirs + d_vers_nw;
% r should be defined = r0 + t * (h + y) / (1 + h . y)
r_dir = rep_r_dirs + d_vers_nw;
else
y_lab = gtGeoPlLabFromThetaEta(t(ii_b * ones_ns), ns', p.labgeo);
rot_ws = gtMathsRotationTensor(ws, rotcomp_w);
rot_ws = reshape(rot_ws, 3, []);
y_sam = y_lab * rot_ws;
y_sam = reshape(y_sam, oversampling_n * num_ns(ii_b), 3, oversampling_w * num_ws(ii_b));
y_sam = permute(y_sam, [1 3 2]);
y_sam = reshape(y_sam, [], 3);
renorms = 1 + y_sam * h(ii_b, :)';
renorms = renorms(:, [1 1 1]);
hh = h(ii_b * ones_tvecs, :);
% r0 is defined = (h x y) / (1 + h . y)
r0 = gtMathsCross(hh, y_sam) ./ renorms - gr.R_vector(ones_tvecs, :);
% r is defined = r0 + t * (h + y) / (1 + h . y)
r_dir = (hh + y_sam) ./ renorms;
% r_dir = gtMathsNormalizeVectorsList(r_dir);
end
% Translating r0s to the closest point to the average R
% vector, on the line r
......@@ -401,26 +417,19 @@ function sf = gtDefShapeFunctionsCreateNW(sampler, eta_resoltion, varargin)
fprintf(repmat(' ', [1, num_chars]));
fprintf(repmat('\b', [1, num_chars]));
fprintf('Done in %g seconds.\n', toc(c))
end
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';
function [r0, r_dir] = compute_rod_line(thetas, etas, rot_w, h, labgeo, samgeo)
y_lab = gtGeoPlLabFromThetaEta(thetas, etas, labgeo);
y_sam = gtGeoLab2Sam(y_lab, rot_w, labgeo, samgeo, true);
renorms = 1 + sum(y_w .* h_12, 2);
renorms = renorms(:, [1 1 1]);
renorms = 1 + sum(y_sam .* h, 2);
renorms = renorms(:, [1 1 1]);
% r0 is defined = (h x y) / (1 + h . y)
r0 = gtMathsCross(h_12, y_w) ./ renorms;
% r0 is defined = (h x y) / (1 + h . y)
r0 = gtMathsCross(h, y_sam) ./ renorms;
% r is defined = r0 + t * (h + y) / (1 + h . y)
r_dir = (h_12 + y_w) ./ renorms;
r_dir = gtMathsNormalizeVectorsList(r_dir);
end
% r is defined = r0 + t * (h + y) / (1 + h . y)
r_dir = (h + y_sam) ./ renorms;
r_dir = gtMathsNormalizeVectorsList(r_dir);
end
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