From 8d42d06e85593526772db81338ecc8152dc192dc Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Fri, 6 May 2016 18:20:09 +0200
Subject: [PATCH] 6D-reconstruction: enabling ospace-super-sampling with the
use of W-shape-functions
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Deformation/Gt6DBlobReconstructor.m | 28 ++++++--
.../Gt6DReconstructionAlgorithmFactory.m | 44 ++++++++++--
zUtil_Deformation/GtOrientationSampling.m | 28 +++++---
.../gtDefShapeFunctionsCreateW.m | 70 ++++++++++++++-----
4 files changed, 133 insertions(+), 37 deletions(-)
diff --git a/zUtil_Deformation/Gt6DBlobReconstructor.m b/zUtil_Deformation/Gt6DBlobReconstructor.m
index c1a352c0..1fdf06b9 100644
--- a/zUtil_Deformation/Gt6DBlobReconstructor.m
+++ b/zUtil_Deformation/Gt6DBlobReconstructor.m
@@ -289,10 +289,18 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
num_bytes = 0;
float_bytes = GtBenchmarks.getSizeVariable(zeros(1, 1, 'single'));
num_geoms = self.get_number_geometries();
- for n = 1:num_geoms
- sinogram_size = [ ...
- self.proj_size(1), size(self.geometries{n}, 1), self.proj_size(2)];
- num_bytes = num_bytes + float_bytes * prod(sinogram_size);
+ if (isempty(self.ss_geometries))
+ for n = 1:num_geoms
+ sinogram_size = [ ...
+ self.proj_size(1), size(self.geometries{n}, 1), self.proj_size(2)];
+ num_bytes = num_bytes + float_bytes * prod(sinogram_size);
+ end
+ else
+ for n = 1:num_geoms
+ sinogram_size = [ ...
+ self.proj_size(1), size(self.ss_geometries{n}, 1), self.proj_size(2)];
+ num_bytes = num_bytes + float_bytes * prod(sinogram_size);
+ end
end
end
end
@@ -514,8 +522,16 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
fprintf(' + Backprojecting ones sinos (fake bproj)..');
c = tic();
num_geoms = self.get_number_geometries();
- for ii = 1:num_geoms
- self.bwd_weights{ii} = size(self.geometries{ii}, 1);
+ if (~isempty(self.geometries))
+ for ii = 1:num_geoms
+ self.bwd_weights{ii} = size(self.geometries{ii}, 1);
+ end
+ else
+ % If this happens we are using oversampling with
+ % shape-functions
+ for ii = 1:num_geoms
+ self.bwd_weights{ii} = numel(self.blobs);
+ end
end
fprintf('\b\b (%2.1f s)\n', toc(c));
end
diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
index 694214bf..0142bee3 100644
--- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
+++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
@@ -447,11 +447,18 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
sel_allb_idxs = sampler.ondet(sampler.included(sel_bls));
omegas = ref_gr.allblobs.omega(sel_allb_idxs);
- orients = sampler.get_orientations();
- orients = [orients{:}];
+ [orients, orients_ss] = sampler.get_orientations();
+
+ num_ospace_oversampling = prod(sampler.ss_factor);
+ if (num_ospace_oversampling == 1)
+ ors_struct = [orients{:}];
+ else
+ ors_struct = cat(4, orients_ss{:});
+ ors_struct = [ors_struct{:}];
+ end
% Array of structures that contain all the info relative to
% each orientation, for each blob
- grains_props = self.get_selected_props(orients, sel_bls);
+ grains_props = self.get_selected_props(ors_struct, sel_bls);
[w_tab, ref_ws] = self.get_shape_funcs_w_tab(shape_funcs, omegas);
[lims_blobs_w, lims_projs_w, paddings, ext_blobs_w_orig] = self.get_blob_lims(w_tab, bl, ref_ws);
@@ -494,14 +501,39 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
psf = {};
end
+ tot_orient = sampler.get_total_num_orientations();
+ if (num_ospace_oversampling == 1)
+ geometries_ss = {};
+ offsets_ss = {};
+ else
+ geometries = reshape(geometries, num_ospace_oversampling, []);
+ offsets = reshape(offsets, num_ospace_oversampling, []);
+
+ geometries_ss = cell(tot_orient, 1);
+ offsets_ss = cell(tot_orient, 1);
+ for ii = 1:tot_orient
+ geometries_ss{ii} = geometries(:, ii);
+
+ % merging all the orientation-space super-sampling
+ % projection coefficients into one structure per
+ % orientation. This is needed in order to fwd/bwd-project
+ % all the different sub-orientations together.
+ offsets_ss{ii} = self.merge_offsets(offsets(:, ii));
+ end
+
+ geometries = {};
+ offsets = {};
+ end
+
algo_params = struct( ...
'blobs', {blobs}, ...
'geometries', {geometries}, ...
'offsets', {offsets}, ...
- 'geometries_ss', {{}}, ...
- 'offsets_ss', {{}}, ...
+ 'geometries_ss', {geometries_ss}, ...
+ 'offsets_ss', {offsets_ss}, ...
'psf', {psf}, ...
- 'good_orients', {true(1, numel(geometries))} );
+ 'blobs_w_lims', lims_blobs_w, ...
+ 'good_orients', {true(1, tot_orient)} );
end
function [sub_blob_slices, proj_coeffs] = get_sub_blobs(self, blobs, slices_interp, padding)
diff --git a/zUtil_Deformation/GtOrientationSampling.m b/zUtil_Deformation/GtOrientationSampling.m
index f544e63a..e8344cba 100644
--- a/zUtil_Deformation/GtOrientationSampling.m
+++ b/zUtil_Deformation/GtOrientationSampling.m
@@ -197,15 +197,24 @@ classdef GtOrientationSampling < handle
gr = self.lattice(ii_o).gr{ii_g};
beg_pos = gr.R_vector - half_voxel_size + half_steps;
end_pos = gr.R_vector + half_voxel_size - half_steps;
- gr_ss = {};
- for ii_x = beg_pos(1):steps(1):end_pos(1)
- for ii_y = beg_pos(2):steps(2):end_pos(2)
- for ii_z = beg_pos(3):steps(3):end_pos(3)
- gr_ss{end+1} = struct( ...
- 'R_vector', {[ii_x, ii_y, ii_z]}, ...
+
+ range_x = beg_pos(1):steps(1):end_pos(1);
+ range_y = beg_pos(2):steps(2):end_pos(2);
+ range_z = beg_pos(3):steps(3):end_pos(3);
+
+ gr_ss = cell(self.ss_factor);
+ for ii_z = 1:self.ss_factor(3)
+ for ii_y = 1:self.ss_factor(2)
+ r_vecs = [ ...
+ range_x', ...
+ range_y(ii_y) * ones(self.ss_factor(1), 1), ...
+ range_z(ii_z) * ones(self.ss_factor(1), 1)];
+ for ii_x = 1:self.ss_factor(1)
+ gr_ss{ii_x, ii_y, ii_z} = struct( ...
+ 'R_vector', {r_vecs(ii_x, :)}, ...
'center', {gr.center}, ...
'phaseid', {gr.phaseid}, ...
- 'allblobs', {[]} ); %#ok<AGROW>
+ 'allblobs', {[]} );
end
end
end
@@ -248,7 +257,10 @@ classdef GtOrientationSampling < handle
num_lattices = numel(self.lattice);
or_size = cell(num_lattices, 1);
for ii = 1:num_lattices
- or_size{ii} = size(self.lattice(ii).gr);
+ or_size{ii} = [ ...
+ size(self.lattice(ii).gr, 1), ...
+ size(self.lattice(ii).gr, 2), ...
+ size(self.lattice(ii).gr, 3) ];
end
end
diff --git a/zUtil_Deformation/gtDefShapeFunctionsCreateW.m b/zUtil_Deformation/gtDefShapeFunctionsCreateW.m
index b40eabfe..dad67c49 100644
--- a/zUtil_Deformation/gtDefShapeFunctionsCreateW.m
+++ b/zUtil_Deformation/gtDefShapeFunctionsCreateW.m
@@ -7,30 +7,49 @@ function sf = gtDefShapeFunctionsCreateW(sampler)
omstep = gtGetOmegaStepDeg(p, sampler.detector_index);
omstep_rod = tand(omstep / 2);
- sub_space_size = sampler.stats.sampling.gaps;
-
- % Let's first compute the W extent
- half_rspace_sizes = sub_space_size / 2;
-
ref_incs = sampler.ondet(sampler.included);
ref_inds = sampler.selected;
num_blobs = numel(find(ref_inds));
ref_gr = sampler.get_reference_grain();
- pl_samd = ref_gr.allblobs.plorig(ref_incs(ref_inds), :);
+ if (isfield(ref_gr.allblobs, 'plcry'))
+ pl_cry = ref_gr.allblobs.plcry(ref_incs(ref_inds), :);
+ else
+ pl_samd = ref_gr.allblobs.plorig(ref_incs(ref_inds), :);
+ g = gtMathsRod2OriMat(ref_gr.R_vector');
+ pl_cry = gtVectorLab2Cryst(pl_samd, g)';
+ end
+% pl_labd = gtGeoSam2Lab(pl_samd, eye(3), self.parameters.labgeo, self.parameters.samgeo, true, false)';
+
+ sub_space_size = sampler.stats.sampling.gaps;
+ ospace_samp_size = sampler.get_orientation_sampling_size();
+ ospace_samp_size = ospace_samp_size{1};
+
+ num_ospace_oversampling = prod(sampler.ss_factor);
+ if (num_ospace_oversampling == 1)
+ ospace_samp_lattice = sampler.lattice.gr;
+ else
+ ospace_samp_lattice = cat(4, sampler.lattice_ss.gr{:});
+ ospace_samp_lattice = reshape(ospace_samp_lattice, [sampler.ss_factor, ospace_samp_size]);
+ ospace_samp_lattice = permute(ospace_samp_lattice, [1 4 2 5 3 6]);
+
+ ospace_samp_size = ospace_samp_size .* sampler.ss_factor;
+ sub_space_size = sub_space_size ./ sampler.ss_factor;
+
+ ospace_samp_lattice = reshape(ospace_samp_lattice, ospace_samp_size);
+ end
- 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)';
+ % Let's first compute the W extent
+ half_rspace_sizes = sub_space_size / 2;
fprintf('Computing bounds of W shape functions: ')
c = tic();
- tot_orientations = numel(sampler.lattice.gr);
+ tot_orientations = numel(ospace_samp_lattice);
- oris_bounds_size = [size(sampler.lattice.gr, 1), size(sampler.lattice.gr, 2), size(sampler.lattice.gr, 3)] + 1;
- oris_lims_min = sampler.lattice.gr{1, 1, 1}.R_vector - half_rspace_sizes;
- oris_lims_max = sampler.lattice.gr{end, end, end}.R_vector + half_rspace_sizes;
+ oris_bounds_size = ospace_samp_size + 1;
+ oris_lims_min = ospace_samp_lattice{1, 1, 1}.R_vector - half_rspace_sizes;
+ oris_lims_max = ospace_samp_lattice{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));
@@ -131,7 +150,7 @@ function sf = gtDefShapeFunctionsCreateW(sampler)
num_chars = fprintf('%05d/%05d', ii_g, tot_orientations);
end
- gr = sampler.lattice.gr{ii_g};
+ or = ospace_samp_lattice{ii_g};
% 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);
@@ -145,7 +164,7 @@ function sf = gtDefShapeFunctionsCreateW(sampler)
% 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_tot_sp, :)' + r_vecs;
+ r_vecs = or.R_vector(ones_tot_sp, :)' + r_vecs;
gs = gtMathsRod2OriMat(r_vecs);
% Unrolling and transposing the matrices at the same time
@@ -169,10 +188,10 @@ function sf = gtDefShapeFunctionsCreateW(sampler)
D = Ac .^ 2 + As .^ 2;
% Precomputing useful values, like sinthetas
- ominds = gr.allblobs.omind(ref_inds);
+ ominds = or.allblobs.omind(ref_inds);
ssp = ((ominds == 1) | (ominds == 2));
ss = ssp - ~ssp;
- sinths = ss .* gr.allblobs.sintheta(ref_inds);
+ sinths = ss .* or.allblobs.sintheta(ref_inds);
sinths = sinths(:, ones(tot_sampling_points, 1));
CC = Acc + sinths;
@@ -235,4 +254,21 @@ function sf = gtDefShapeFunctionsCreateW(sampler)
end
end
fprintf('Done in %g seconds.\n', toc(c))
+
+ if (num_ospace_oversampling > 1)
+ fprintf('Re-grouping oversampled W shape functions: ')
+ c = tic();
+ ospace_samp_size = sampler.get_orientation_sampling_size();
+ ospace_samp_size = ospace_samp_size{1};
+ tot_orient = sampler.get_total_num_orientations();
+
+ regroup = [ ...
+ sampler.ss_factor(1), ospace_samp_size(1), ...
+ sampler.ss_factor(2), ospace_samp_size(2), ...
+ sampler.ss_factor(3), ospace_samp_size(3) ];
+ sf = reshape(sf, regroup);
+ sf = permute(sf, [1 3 5 2 4 6]);
+ sf = reshape(sf, num_ospace_oversampling, tot_orient);
+ fprintf('Done in %g seconds.\n', toc(c))
+ end
end
--
GitLab