From fbc0cfb847d55ed0bd7a8db85b10258b0acf782e Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Tue, 23 Feb 2016 18:55:14 +0100
Subject: [PATCH] ODF-uvw-solver: added use of w-shape-functions
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Deformation/GtGrainODFuvwSolver.m | 193 ++++++++++++++++++++++--
1 file changed, 178 insertions(+), 15 deletions(-)
diff --git a/zUtil_Deformation/GtGrainODFuvwSolver.m b/zUtil_Deformation/GtGrainODFuvwSolver.m
index 54c7128d..5058697f 100644
--- a/zUtil_Deformation/GtGrainODFuvwSolver.m
+++ b/zUtil_Deformation/GtGrainODFuvwSolver.m
@@ -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)
--
GitLab