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

ODF-uvw-solver: added use of w-shape-functions 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 106362d0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/GtGrainODFuvwSolver.m
+ 178
− 15
View file @ fbc0cfb8
......@@ -19,14 +19,31 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
blob_dephs = arrayfun(@(x)size(x.intm, 3), bls);
blob_dephs = reshape(blob_dephs, [], 1);
delta_omegas = self.sampler.get_omega_deviations();
blob_lims = cat(1, bls(:).bbwim);
with_shape_functions = ~isempty(self.shape_functions);
if (with_shape_functions)
for ii = numel(self.shape_functions):-1:1
proj_lims(:, :, ii) = cat(1, self.shape_functions{ii}(:).bbwim);
end
proj_lims = [min(proj_lims(:, 1, :), [], 3), max(proj_lims(:, 2, :), [], 3)];
delta_omegas = proj_lims(:, 2) - proj_lims(:, 1) + 1;
else
[delta_omegas, proj_lims] = self.sampler.get_omega_deviations(with_shape_functions);
end
chosen_depts = max(blob_dephs, delta_omegas(self.sampler.selected));
num_ws = max(chosen_depts) + 2;
self.size_sino = [size(bls(1).intm, 1), size(bls(1).intm, 2), num_ws, num_blobs];
self.sino = zeros(self.size_sino);
self.pre_paddings = floor((num_ws - blob_dephs) / 2) + 1;
% In the case that the sampled orientations determine the shift, we
% have to take it into account!
additional_shift = blob_lims(:, 1) - proj_lims(:, 1);
additional_shift(additional_shift < 0) = 0;
self.pre_paddings = floor((num_ws - chosen_depts) / 2) + 1 + additional_shift;
for ii_b = 1:num_blobs
ints_interval = self.pre_paddings(ii_b):(self.pre_paddings(ii_b) + blob_dephs(ii_b) -1);
......@@ -69,11 +86,16 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
end
end
function build_projection_matrices(self)
bls = self.sampler.bl(self.sampler.selected);
function build_projection_matrices_sf_none(self)
ref_gr = self.sampler.get_reference_grain();
ref_omind = ref_gr.allblobs.omind;
ref_inc = self.sampler.ondet(self.sampler.included);
ref_sel = self.sampler.selected;
om_step = gtGetOmegaStepDeg(self.parameters, self.sampler.detector_index);
bls = self.sampler.bl(ref_sel);
fprintf('Computing projection matrices: ')
c = tic();
num_ws = self.get_num_ws();
bls_bbws = cat(1, bls(:).bbwim);
......@@ -99,19 +121,15 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
self.St = cell(num_voxels, 1);
shift_grs = cell(num_voxels, 1);
ref_gr = self.sampler.get_reference_grain();
% computing shifts
for ii_v = 1:num_voxels
shift_grs{ii_v} = ref_gr;
shift_grs{ii_v}.center = self.voxel_centers(ii_v, :);
end
ref_omind = ref_gr.allblobs.omind;
ref_inc = self.sampler.ondet(self.sampler.included);
shift_grs = gtCalculateGrain_p(shift_grs, self.parameters, ...
'ref_omind', ref_omind, 'included', ref_inc);
ref_uv = ref_gr.allblobs.detector(det_ind).uvw(ref_inc(self.sampler.selected), 1:2);
ref_uv = ref_gr.allblobs.detector(det_ind).uvw(ref_inc(ref_sel), 1:2);
for ii_v = 1:num_voxels
b_us = [];
......@@ -121,14 +139,14 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
b_is = [];
b_os = [];
shift_uv = shift_grs{ii_v}.allblobs.detector(det_ind).uvw(bl_selected, 1:2) - ref_uv;
shift_uv = shift_grs{ii_v}.allblobs.detector(det_ind).uvw(ref_sel, 1:2) - ref_uv;
num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
for ii_o = 1:num_orients
ab = grid_gr{ii_o}.allblobs;
% W part
ws = ab.omega(self.sampler.selected) / om_step;
ws = ab.omega(ref_sel) / om_step;
min_ws = floor(ws);
max_ws = min_ws + 1;
......@@ -137,7 +155,7 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
min_w_cs = 1 - max_w_cs;
% UV part
uv = ab.detector(det_ind).uvw(self.sampler.selected, 1:2) + shift_uv;
uv = ab.detector(det_ind).uvw(ref_sel, 1:2) + shift_uv;
num_pos = size(uv, 1);
min_uvs = floor(uv);
......@@ -202,7 +220,152 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
fprintf(repmat('\b', [1 num_chars]));
end
fprintf('Done in %f seconds.\n', toc(c));
end
function build_projection_matrices_sf_w(self)
ref_gr = self.sampler.get_reference_grain();
ref_omind = ref_gr.allblobs.omind;
ref_inc = self.sampler.ondet(self.sampler.included);
ref_sel = self.sampler.selected;
bls = self.sampler.bl(self.sampler.selected);
num_ws = self.get_num_ws();
bls_bbws = cat(1, bls(:).bbwim);
min_w_conds = bls_bbws(:, 1) - self.pre_paddings + 1;
max_w_conds = min_w_conds + num_ws - 1;
bls_bbus = cat(1, bls(:).bbuim);
bls_bbvs = cat(1, bls(:).bbvim);
grid_gr = self.sampler.get_orientations();
num_orients = numel(grid_gr);
det_ind = self.sampler.detector_index;
% We make now the approximation that by moving in XYZ inside the
% sample, the projection moves like if the incidence was always
% perpendicular. In reality, given small tilts of the detector or
% small scattering angles, this is hardly violated.
self.chose_voxel_centers();
num_voxels = size(self.voxel_centers, 1);
self.S = cell(num_voxels, 1);
self.St = cell(num_voxels, 1);
shift_grs = cell(num_voxels, 1);
% computing shifts
for ii_v = 1:num_voxels
shift_grs{ii_v} = ref_gr;
shift_grs{ii_v}.center = self.voxel_centers(ii_v, :);
end
shift_grs = gtCalculateGrain_p(shift_grs, self.parameters, ...
'ref_omind', ref_omind, 'included', ref_inc);
ref_uv = ref_gr.allblobs.detector(det_ind).uvw(ref_inc(ref_sel), 1:2);
fprintf('\b\b: ')
for ii_v = 1:num_voxels
b_us = [];
b_vs = [];
b_ws = [];
b_cs = [];
b_is = [];
b_os = [];
shift_uv = shift_grs{ii_v}.allblobs.detector(det_ind).uvw(ref_sel, 1:2) - ref_uv;
num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
for ii_o = 1:num_orients
ab = grid_gr{ii_o}.allblobs;
sf = self.shape_functions{ii_o};
w_cs = cat(3, sf(:).intm);
w_cs = reshape(w_cs, [], 1);
num_blobs = numel(sf);
ws = cell(num_blobs, 1);
w_is = cell(num_blobs, 1);
for ii_b = 1:num_blobs
ws{ii_b} = sf(ii_b).bbwim(1):sf(ii_b).bbwim(2);
w_is{ii_b} = ii_b(ones(sf(ii_b).bbsize(3), 1));
end
ws = reshape([ws{:}], [], 1);
w_is = cat(1, w_is{:});
% UV part
uv = ab.detector(det_ind).uvw(ref_sel, 1:2) + shift_uv;
min_uvs = floor(uv);
max_uvs = min_uvs + 1;
max_uv_cs = uv - min_uvs;
min_uv_cs = 1 - max_uv_cs;
% Acceptance criteria
ok_u_mins = (min_uvs(:, 1) >= bls_bbus(:, 1)) & (min_uvs(:, 1) <= bls_bbus(:, 2));
ok_u_maxs = (max_uvs(:, 1) >= bls_bbus(:, 1)) & (max_uvs(:, 1) <= bls_bbus(:, 2)) & (max_uv_cs(:, 1) > eps('single'));
ok_v_mins = (min_uvs(:, 2) >= bls_bbvs(:, 1)) & (min_uvs(:, 2) <= bls_bbvs(:, 2));
ok_v_maxs = (max_uvs(:, 2) >= bls_bbvs(:, 1)) & (max_uvs(:, 2) <= bls_bbvs(:, 2)) & (max_uv_cs(:, 2) > eps('single'));
u_oks = [ok_u_mins, ok_u_maxs];
v_oks = [ok_v_mins, ok_v_maxs];
u_cs = [min_uv_cs(:, 1), max_uv_cs(:, 1)];
v_cs = [min_uv_cs(:, 2), max_uv_cs(:, 2)];
pos_us = [(min_uvs(:, 1) - bls_bbus(:, 1) + 1), (max_uvs(:, 1) - bls_bbus(:, 1) + 1)];
pos_vs = [(min_uvs(:, 2) - bls_bbvs(:, 1) + 1), (max_uvs(:, 2) - bls_bbvs(:, 1) + 1)];
pos_ws = ws - min_w_conds(w_is) + 1;
uvw_oks = u_oks(w_is, [1 1 2 2]) & v_oks(w_is, [1 2 1 2]);
uvw_cs = u_cs(w_is, [1 1 2 2]) .* v_cs(w_is, [1 2 1 2]) .* w_cs(:, [1 1 1 1]);
pos_us = pos_us(w_is, [1 1 2 2]);
pos_vs = pos_vs(w_is, [1 2 1 2]);
pos_ws = pos_ws(:, [1 1 1 1]);
indx = find(uvw_oks);
indx_mod = mod(indx - 1, numel(w_is)) + 1;
wrong_w = ws < min_w_conds(w_is) | ws > max_w_conds(w_is);
if (any(wrong_w))
ii_o
find(wrong_w)
[min_w_conds(w_is(wrong_w)), ws(wrong_w), max_w_conds(w_is(wrong_w)), ...
(max_w_conds(w_is(wrong_w)) - min_w_conds(w_is(wrong_w)) + 1), ...
(ws(wrong_w) - min_w_conds(w_is(wrong_w)) + 1) ]
end
b_us = [b_us; pos_us(indx)];
b_vs = [b_vs; pos_vs(indx)];
b_ws = [b_ws; pos_ws(indx)];
b_cs = [b_cs; uvw_cs(indx)];
b_is = [b_is; w_is(indx_mod)];
b_os = [b_os; ii_o(ones(numel(indx), 1))];
end
% sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
sino_indx = b_us ...
+ self.size_sino(1) * (b_vs - 1 ...
+ self.size_sino(2) * (b_ws - 1 ...
+ self.size_sino(3) * (b_is - 1)));
self.S{ii_v} = sparse( ...
sino_indx, b_os, b_cs, ...
numel(self.sino), num_orients);
self.St{ii_v} = self.S{ii_v}';
fprintf(repmat('\b', [1 num_chars]));
end
end
function vol = get_volume(self)
......
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