From 21156f8869481a06f945b541de1332731a495787 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Wed, 17 Feb 2016 16:26:29 +0100
Subject: [PATCH] ODF Solvers: round of improvements, and added real-space
sampling to uvw solver
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Deformation/GtGrainODFuvwSolver.m | 235 ++++++++++++------
zUtil_Deformation/GtGrainODFwSolver.m | 133 ++++------
.../gtDefComputeGrainODFNearField.m | 2 +-
3 files changed, 205 insertions(+), 165 deletions(-)
diff --git a/zUtil_Deformation/GtGrainODFuvwSolver.m b/zUtil_Deformation/GtGrainODFuvwSolver.m
index f0fad513..b1351326 100644
--- a/zUtil_Deformation/GtGrainODFuvwSolver.m
+++ b/zUtil_Deformation/GtGrainODFuvwSolver.m
@@ -1,22 +1,26 @@
classdef GtGrainODFuvwSolver < GtGrainODFwSolver
properties
+ voxel_centers = [];
end
methods (Access = public)
- function self = GtGrainODFuvwSolver(parameters)
- self = self@GtGrainODFwSolver(parameters);
+ function self = GtGrainODFuvwSolver(parameters, varargin)
+ self = self@GtGrainODFwSolver(parameters, varargin{:});
end
end
methods (Access = public) % Low Level API
- function build_sinogram(self, bls, oversize)
+ function build_sinogram(self)
+ bls = self.sampler.bl(self.sampler.selected);
+
num_blobs = numel(bls);
blob_dephs = arrayfun(@(x)size(x.intm, 3), bls);
blob_dephs = reshape(blob_dephs, [], 1);
- num_ws = max(blob_dephs) + 2;
- if (exist('oversize', 'var'))
- num_ws = round(num_ws * oversize);
- end
+
+ delta_omegas = self.sampler.get_omega_deviations();
+ 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);
@@ -29,22 +33,41 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
self.sino(:, :, ints_interval, ii_b) = masked_blob;
end
self.sino = reshape(self.sino, [], 1);
+
+ if (self.verbose)
+ fprintf('Sino size: [%d, %d, %d, %d]\n', self.size_sino);
+ end
end
- function deg_res = guess_resolution(self, ref_gr, sampler)
- deg_res_omega = guess_resolution@GtGrainODFwSolver(self, ref_gr, sampler);
+ function chose_voxel_centers(self, over_sampling)
+ if (~exist('over_sampling', 'var') || isempty(over_sampling))
+ over_sampling = 1;
+ end
+
+ recgeo = self.parameters.recgeo(self.sampler.detector_index);
+ ref_gr = self.sampler.get_reference_grain();
+ proj = ref_gr.proj(self.sampler.detector_index);
+ r_space_vol_size = [proj.vol_size_y, proj.vol_size_x, proj.vol_size_z];
+
+ max_dists_from_center = (r_space_vol_size - 1) / 2;
+ num_voxels = prod(r_space_vol_size);
+
+ self.voxel_centers = zeros(num_voxels, 3);
+
+ counter = 1;
- bls = sampler.bl(sampler.selected);
+ for ss_x = linspace(-max_dists_from_center(1), max_dists_from_center(1), r_space_vol_size(1))
+ for ss_y = linspace(-max_dists_from_center(2), max_dists_from_center(2), r_space_vol_size(2))
+ for ss_z = linspace(-max_dists_from_center(3), max_dists_from_center(3), r_space_vol_size(3))
- bl_centers = [sum(cat(1, bls(:).bbuim), 2), sum(cat(1, bls(:).bbvim), 2)] / 2;
- gc_det_pos = gtGeoGrainCenterSam2Det(ref_gr.center, ref_gr.allblobs.omega, self.parameters);
- dists = sqrt(sum((bl_centers - gc_det_pos) .^ 2, 2));
+ offset_voxel = [ss_x, ss_y, ss_z] .* recgeo.voxsize;
- % Higest resolution in ETA
- deg_res_eta = min(2 .* atan2d(0.5, dists));
+ self.voxel_centers(counter, :) = ref_gr.center + offset_voxel;
- % We ake the lowest resolution between the two
- deg_res = max(deg_res_eta, deg_res_omega);
+ counter = counter + 1;
+ end
+ end
+ end
end
function build_projection_matrices(self, bls, bl_selected)
@@ -53,12 +76,10 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
else
om_step = self.parameters.labgeo.omstep;
end
- fprintf('Computing projection matrices..')
+ fprintf('Computing projection matrices: ')
c = tic();
num_ws = self.get_num_ws();
- num_orients = numel(self.grid_gr);
-
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;
@@ -66,87 +87,139 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
bls_bbus = cat(1, bls(:).bbuim);
bls_bbvs = cat(1, bls(:).bbvim);
- b_us = [];
- b_vs = [];
- b_ws = [];
- b_cs = [];
- b_is = [];
- b_os = [];
-
- for ii_o = 1:num_orients
- ab = self.grid_gr{ii_o}.allblobs;
- % W part
- ws = ab.omega(bl_selected) / om_step;
-
- min_ws = floor(ws);
- max_ws = min_ws + 1;
-
- max_w_cs = ws - min_ws;
- min_w_cs = 1 - max_w_cs;
-
- % UV part
- if (isfield(ab, 'detector'))
- uv = ab.detector(1).uvw(bl_selected, 1:2);
- else
- uv = ab.uvw(bl_selected, 1:2);
- end
- min_uvs = floor(uv);
- max_uvs = min_uvs + 1;
+ 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);
+ 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(bl_selected), 1:2);
+
+ 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(bl_selected, 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(bl_selected) / om_step;
- max_uv_cs = uv - min_uvs;
- min_uv_cs = 1 - max_uv_cs;
+ min_ws = floor(ws);
+ max_ws = min_ws + 1;
- % Acceptance criteria
- ok_w_mins = (min_ws >= min_w_conds) & (min_ws <= max_w_conds);
- ok_w_maxs = (max_ws >= min_w_conds) & (max_ws <= max_w_conds) & (max_w_cs > eps('single'));
+ max_w_cs = ws - min_ws;
+ min_w_cs = 1 - max_w_cs;
- 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'));
+ % UV part
+ uv = ab.detector(det_ind).uvw(bl_selected, 1:2) + shift_uv;
- 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'));
+ min_uvs = floor(uv);
+ max_uvs = min_uvs + 1;
- u_oks = [ok_u_mins, ok_u_maxs];
- v_oks = [ok_v_mins, ok_v_maxs];
- w_oks = [ok_w_mins, ok_w_maxs];
+ max_uv_cs = uv - min_uvs;
+ min_uv_cs = 1 - max_uv_cs;
- u_cs = [min_uv_cs(:, 1), max_uv_cs(:, 1)];
- v_cs = [min_uv_cs(:, 2), max_uv_cs(:, 2)];
- w_cs = [min_w_cs, max_w_cs];
+ % Acceptance criteria
+ ok_w_mins = (min_ws >= min_w_conds) & (min_ws <= max_w_conds);
+ ok_w_maxs = (max_ws >= min_w_conds) & (max_ws <= max_w_conds) & (max_w_cs > eps('single'));
- 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 = [(min_ws - bls_bbws(:, 1) + self.pre_paddings), (max_ws - bls_bbws(:, 1) + self.pre_paddings)];
+ 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'));
- for u_ii = 1:2
- for v_ii = 1:2
- for w_ii = 1:2
- indx = find(u_oks(:, u_ii) & v_oks(:, v_ii) & w_oks(:, w_ii));
- b_us = [b_us; pos_us(indx, u_ii)];
- b_vs = [b_vs; pos_vs(indx, v_ii)];
- b_ws = [b_ws; pos_ws(indx, w_ii)];
+ 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'));
- b_cs = [b_cs; u_cs(indx, u_ii) .* v_cs(indx, v_ii) .* w_cs(indx, w_ii)];
+ u_oks = [ok_u_mins, ok_u_maxs];
+ v_oks = [ok_v_mins, ok_v_maxs];
+ w_oks = [ok_w_mins, ok_w_maxs];
- b_is = [b_is; indx];
- b_os = [b_os; ii_o(ones(numel(indx), 1))];
+ u_cs = [min_uv_cs(:, 1), max_uv_cs(:, 1)];
+ v_cs = [min_uv_cs(:, 2), max_uv_cs(:, 2)];
+ w_cs = [min_w_cs, max_w_cs];
+
+ 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 = [(min_ws - bls_bbws(:, 1) + self.pre_paddings), (max_ws - bls_bbws(:, 1) + self.pre_paddings)];
+
+ for u_ii = 1:2
+ for v_ii = 1:2
+ for w_ii = 1:2
+ indx = find(u_oks(:, u_ii) & v_oks(:, v_ii) & w_oks(:, w_ii));
+ b_us = [b_us; pos_us(indx, u_ii)];
+ b_vs = [b_vs; pos_vs(indx, v_ii)];
+ b_ws = [b_ws; pos_ws(indx, w_ii)];
+
+ b_cs = [b_cs; u_cs(indx, u_ii) .* v_cs(indx, v_ii) .* w_cs(indx, w_ii)];
+
+ b_is = [b_is; indx];
+ b_os = [b_os; ii_o(ones(numel(indx), 1))];
+ end
end
end
end
- end
- sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
+ sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
- self.S = sparse( ...
- sino_indx, b_os, b_cs, ...
- numel(self.sino), num_orients);
+ self.S{ii_v} = sparse( ...
+ sino_indx, b_os, b_cs, ...
+ numel(self.sino), num_orients);
- self.St = self.S';
- fprintf('\b\b: Done in %f seconds.\n', toc(c));
+ self.St{ii_v} = self.S{ii_v}';
+
+ fprintf(repmat('\b', [1 num_chars]));
+ end
+ fprintf('Done in %f seconds.\n', toc(c));
+ end
+
+ function vol = get_volume(self)
+ vol = reshape(sum(self.volume, 2), self.size_volume);
end
end
methods (Access = protected)
+ function y = fp(self, x)
+ y = self.S{1} * x(:, 1);
+ for ii_s = 2:numel(self.S)
+ y = y + self.S{ii_s} * x(:, ii_s);
+ end
+ end
+
+ function y = bp(self, x)
+ for ii_s = numel(self.St):-1:1
+ y(:, ii_s) = self.St{ii_s} * x;
+ end
+ end
+
function num_ws = get_num_ws(self)
num_ws = self.size_sino(3);
end
diff --git a/zUtil_Deformation/GtGrainODFwSolver.m b/zUtil_Deformation/GtGrainODFwSolver.m
index c70710b7..54d632d0 100644
--- a/zUtil_Deformation/GtGrainODFwSolver.m
+++ b/zUtil_Deformation/GtGrainODFwSolver.m
@@ -2,7 +2,7 @@ classdef GtGrainODFwSolver < handle
properties
parameters;
- grid_gr;
+ sampler;
volume = [];
size_volume = [];
@@ -29,48 +29,37 @@ classdef GtGrainODFwSolver < handle
self = parse_pv_pairs(self, varargin);
end
- function vol = solve_synthetic(self, ref_gr, gvdm, algorithm, lambda)
- sampler = self.build_orientation_sampling_synthetic(ref_gr, 1.3);
- bls = sampler.bl(sampler.selected);
-
- self.build_sinogram(bls, 1.1);
- self.build_projection_matrices(bls, sampler.selected);
+ function vol = solve(self, ref_gr, varargin)
+ % FUNCTION vol = solve(self, ref_gr, varargin)
+ %
+ % INPUT (varargin):
+ % 'algorithm': 'sirt'
+ % 'lambda': 1e-2
+ % 'ospace_oversize': 1.1
+ % 'det_index': 1
+ %
+
+ conf = struct( ...
+ 'algorithm', 'sirt', ...
+ 'lambda', 1e-2, ...
+ 'ospace_oversize', 1.1, ...
+ 'det_index', 1 );
+ conf = parse_pv_pairs(conf, varargin);
+
+ self.build_orientation_sampling(ref_gr, conf.ospace_oversize, conf.det_index);
+ bls = self.sampler.bl(self.sampler.selected);
+
+ self.build_sinogram();
+ self.build_projection_matrices(bls, self.sampler.selected);
self.build_projection_weights();
- if (~exist('algorithm', 'var'))
- algorithm = 'cplsnn';
- end
- switch(lower(algorithm))
- case 'sirt'
- self.solve_sirt();
- case 'cplsnn'
- self.solve_cplsnn();
- case 'cplsl1nn'
- self.solve_cplsl1nn(lambda);
- case 'cpl1nn'
- self.solve_cpl1nn();
- end
- vol = self.get_volume();
- end
-
- function vol = solve_experimental(self, ref_gr, algorithm, lambda)
- sampler = self.build_orientation_sampling_experimental(ref_gr, 1.05);
- bls = sampler.bl(sampler.selected);
-
- self.build_sinogram(bls, 1.2);
- self.build_projection_matrices(bls, sampler.selected);
- self.build_projection_weights();
-
- if (~exist('algorithm', 'var'))
- algorithm = 'cplsnn';
- end
- switch(lower(algorithm))
+ switch(lower(conf.algorithm))
case 'sirt'
self.solve_sirt();
case 'cplsnn'
self.solve_cplsnn();
case 'cplsl1nn'
- self.solve_cplsl1nn(lambda);
+ self.solve_cplsl1nn(conf.lambda);
case 'cpl1nn'
self.solve_cpl1nn();
end
@@ -79,53 +68,24 @@ classdef GtGrainODFwSolver < handle
end
methods (Access = public) % Low Level API
- function sampler = build_orientation_sampling_synthetic(self, ref_gr, oversize)
- sampler = GtOrientationSampling(self.parameters, ref_gr, 'verbose', self.verbose);
- sampler.make_grid_estim_ODF_resoluion('cubic', true, oversize);
- self.grid_gr = sampler.get_orientations();
- self.size_volume = size(sampler.lattice.gr);
+ function build_orientation_sampling(self, ref_gr, oversize, det_index)
+ self.sampler = GtOrientationSampling(self.parameters, ref_gr, ...
+ 'verbose', self.verbose, 'detector_index', det_index);
+ self.sampler.make_grid_estim_ODF_resoluion('cubic', 0, oversize);
+ self.size_volume = size(self.sampler.lattice.gr);
end
- function sampler = build_orientation_sampling_experimental(self, ref_gr, oversize)
- sampler = GtOrientationSampling(self.parameters, ref_gr, 'verbose', self.verbose);
- sampler.make_grid_estim_ODF_resoluion('cubic', true, oversize);
- self.grid_gr = sampler.get_orientations();
- self.size_volume = size(sampler.lattice.gr);
- end
-
- function deg_res = guess_resolution(self, ref_gr, sampler)
- if (~isfield(self.parameters.labgeo, 'omstep'))
- omega_step = gtGetOmegaStepDeg(self.parameters);
- else
- omega_step = self.parameters.labgeo.omstep;
- end
-
- included = sampler.ondet(sampler.included);
- selected = included(sampler.selected);
- selected_etas = ref_gr.allblobs.eta(selected);
- inv_lorentz = abs(cosd(selected_etas));
-
- % Gives the highest resolution in OMEGA
- deg_res = min(omega_step .* inv_lorentz);
- end
+ function build_sinogram(self)
+ bls = self.sampler.bl(self.sampler.selected);
- function build_sinogram(self, bls, oversize)
num_blobs = numel(bls);
-% % Should be worked on!
-% ints_w = arrayfun(@(x){squeeze(sum(sum(x.intm, 1), 2))}, bls);
-% real_blob_limits = cellfun(@(x){...
-% [max(find(x, 1, 'first'), 1), ...
-% min(find(x, 1, 'last'), numel(x))]...
-% }, ints_w);
-% real_blob_limits = cat(1, real_blob_limits{:});
-%
-% blob_dephs = real_blob_limits(:, 2) - real_blob_limits(:, 1) + 1;
blob_dephs = arrayfun(@(x)size(x.intm, 3), bls);
blob_dephs = reshape(blob_dephs, [], 1);
- num_ws = max(blob_dephs) + 2;
- if (exist('oversize', 'var'))
- num_ws = round(num_ws * oversize);
- end
+
+ delta_omegas = self.sampler.get_omega_deviations();
+ chosen_depts = max(blob_dephs, delta_omegas(self.sampler.selected));
+
+ num_ws = max(chosen_depts) + 2;
self.sino = zeros(num_ws, num_blobs);
self.size_sino = size(self.sino);
@@ -138,6 +98,10 @@ classdef GtGrainODFwSolver < handle
self.sino(ints_interval, ii_b) = squeeze(sum(sum(masked_blob, 1), 2));
end
self.sino = reshape(self.sino, [], 1);
+
+ if (self.verbose)
+ fprintf('Sino size: [%d, %d]\n', self.size_sino);
+ end
end
function sino = get_sinogram(self)
@@ -146,7 +110,7 @@ classdef GtGrainODFwSolver < handle
function comp_sino = get_computed_sinogram(self)
if (~isempty(self.volume))
- comp_sino = self.S * self.volume;
+ comp_sino = self.fp(self.volume);
comp_sino = reshape(comp_sino, self.size_sino);
else
error('GtGrainODFSolver:no_reconstruction', ...
@@ -159,11 +123,13 @@ classdef GtGrainODFwSolver < handle
end
function or = get_orientations(self)
- or = reshape(self.grid_gr, self.size_volume);
+ grid_gr = self.sampler.get_orientations();
+ or = reshape(grid_gr, self.size_volume);
end
function r_vecs = get_R_vectors(self)
- r_vecs = [self.grid_gr{:}];
+ grid_gr = self.sampler.get_orientations();
+ r_vecs = [grid_gr{:}];
r_vecs = cat(1, r_vecs(:).R_vector);
end
@@ -177,8 +143,6 @@ classdef GtGrainODFwSolver < handle
c = tic();
num_ws = self.get_num_ws();
- num_orients = numel(self.grid_gr);
-
bls_bbws = cat(1, bls(:).bbwim);
min_conds = bls_bbws(:, 1) - self.pre_paddings + 1;
max_conds = min_conds + num_ws - 1;
@@ -188,8 +152,11 @@ classdef GtGrainODFwSolver < handle
b_is = [];
b_os = [];
+ grid_gr = self.sampler.get_orientations();
+ num_orients = numel(grid_gr);
+
for ii_o = 1:num_orients
- ab = self.grid_gr{ii_o}.allblobs;
+ ab = grid_gr{ii_o}.allblobs;
ws = ab.omega(bl_selected) / om_step;
min_ws = floor(ws);
diff --git a/zUtil_Deformation/gtDefComputeGrainODFNearField.m b/zUtil_Deformation/gtDefComputeGrainODFNearField.m
index ac0e4b22..032c78ec 100644
--- a/zUtil_Deformation/gtDefComputeGrainODFNearField.m
+++ b/zUtil_Deformation/gtDefComputeGrainODFNearField.m
@@ -16,7 +16,7 @@ function [odfw, odfw_R_vectors] = gtDefComputeGrainODFNearField(phase_id, grain_
sprintf('phase_%02d', phase_id));
gr = gtLoadGrain(phase_id, grain_id);
- odfw = sol.solve_experimental(gr);
+ odfw = sol.solve(gr);
odfw_R_vectors = sol.get_R_vectors();
--
GitLab