-
Nicola Vigano authored
Orientation-Sampling: reduced sampling strategies to two types: numpoints and resolution, each with and without ODF estimation Signed-off-by:Nicola Vigano <nicola.vigano@esrf.fr>
Nicola Vigano authoredOrientation-Sampling: reduced sampling strategies to two types: numpoints and resolution, each with and without ODF estimation Signed-off-by:
GtOrientationSampling.m 58.35 KiB
classdef GtOrientationSampling < handle
properties
lattice = struct('gr', {});
lattice_ss = {};
ss_factor = [1 1 1];
stats;
R_vectors;
fedpars = [];
parameters = [];
ref_gr = [];
% redundant info that helps in forward compatibility
bl;
ondet;
included;
selected;
detector_index = 1;
verbose = false;
end
methods
function self = GtOrientationSampling(parameters, ref_gr, varargin)
self.parameters = parameters;
self.ref_gr = ref_gr;
self = parse_pv_pairs(self, varargin);
if (isfield(ref_gr.proj(self.detector_index), 'bl'))
self.bl = ref_gr.proj(self.detector_index).bl;
else
self.bl = ref_gr.bl;
end
if (isfield(ref_gr.proj(self.detector_index), 'selected'))
self.ondet = ref_gr.proj(self.detector_index).ondet;
self.included = ref_gr.proj(self.detector_index).included;
self.selected = ref_gr.proj(self.detector_index).selected;
else
self.ondet = ref_gr.ondet;
self.included = ref_gr.included;
self.selected = ref_gr.selected;
end
end
function make_grid_numpoints_estim_ODF(self, type, edge_points, even, oversize)
gvdm_tmp = self.guess_ODF_BB();
if (~exist('oversize', 'var'))
oversize = [];
end
if (~exist('even', 'var'))
even = [];
end
self.make_grid_numpoints(type, edge_points, gvdm_tmp, even, oversize);
end
function make_grid_numpoints(self, type, edge_points, gvdm, even, oversize)
if (~exist('oversize', 'var'))
oversize = 1.1;
end
if (exist('even', 'var') && ~isempty(even) && even)
[~, ~, dranges] = self.get_deformation_stats(gvdm);
[~, min_ii] = min(dranges); other_iis = [1:min_ii-1 min_ii+1:3];
min_spacing = dranges(min_ii) * oversize / min_edge_points;
other_edge_points = ceil(dranges(other_iis) * oversize / min_spacing);
edge_points(min_ii) = min_edge_points;
edge_points(other_iis(1)) = other_edge_points(1);
edge_points(other_iis(2)) = other_edge_points(2);
elseif (numel(edge_points) == 1)
edge_points = edge_points([1 1 1]);
end
self.make_grid(type, edge_points, gvdm, oversize)
end
% function make_res_even_simple_grid(self, type, min_edge_points, gvdm, oversize)
% [~, ~, dranges] = self.get_deformation_stats(gvdm);
%
% if (~exist('oversize', 'var'))
% oversize = 1.1;
% end
%
% max_resolution = self.estimate_maximum_resolution();
% % max_res_ratios = max_resolution / min(max_resolution);
% % dranges = dranges ./ max_res_ratios;
%
% [~, min_ii] = min(dranges);
% other_iis = [1:min_ii-1 min_ii+1:3];
% if (min_edge_points > 0)
% min_spacing = dranges(min_ii) * oversize / min_edge_points;
% other_edge_points = ceil(dranges(other_iis) * oversize / min_spacing);
% edge_points(min_ii) = min_edge_points;
% edge_points(other_iis(1)) = other_edge_points(1);
% edge_points(other_iis(2)) = other_edge_points(2);
% else
% size_bb = 2 * atand(dranges * oversize);
% chosen_res = max_resolution / (abs(min_edge_points) + (min_edge_points == 0));
% edge_points = ceil(size_bb ./ chosen_res) + 1;
% end
%
% self.make_grid(type, edge_points, gvdm, oversize)
% end
function make_grid_resolution_estim_ODF(self, type, resolution, oversize)
if (~exist('oversize', 'var'))
oversize = 1.1;
end
gvdm_tmp = self.guess_ODF_BB();
self.make_grid_resolution(type, resolution, gvdm_tmp, oversize)
end
function make_grid_resolution(self, type, resolution, gvdm, oversize)
[~, ~, dranges] = self.get_deformation_stats(gvdm);
if (~exist('oversize', 'var'))
oversize = 1.1;
end
resolution_rod = tand(resolution / 2);
bb_size_rod = dranges * oversize;
edge_points = max(ceil(bb_size_rod ./ resolution_rod) + 1, 2);
self.make_grid(type, edge_points, gvdm, oversize)
end
function make_grid(self, type, edge_points, gvdm, oversize)
[dmean, dcenters, dranges] = self.get_deformation_stats(gvdm);
self.print_sampling_details(dranges, oversize, edge_points)
self.lattice = struct('gr', {}); self.R_vectors = {};
type = lower(type);
switch(type)
case 'cubic'
sampling_ranges{3} = GtOrientationSampling.get_cubic_range(edge_points(3), oversize);
sampling_ranges{2} = GtOrientationSampling.get_cubic_range(edge_points(2), oversize);
sampling_ranges{1} = GtOrientationSampling.get_cubic_range(edge_points(1), oversize);
self.produce_R_vectors_even_simple_grid(sampling_ranges, dcenters, dranges)
case 'fcc'
case 'bcc'
sampling_ranges{3} = GtOrientationSampling.get_cubic_range(edge_points(3), oversize);
sampling_ranges{2} = GtOrientationSampling.get_cubic_range(edge_points(2), oversize);
sampling_ranges{1} = GtOrientationSampling.get_cubic_range(edge_points(1), oversize);
self.produce_R_vectors_simple_grid(sampling_ranges, dcenters, dranges, 1)
sampling_ranges{3} = GtOrientationSampling.get_bcc_center_range(edge_points(3), oversize);
sampling_ranges{2} = GtOrientationSampling.get_bcc_center_range(edge_points(2), oversize);
sampling_ranges{1} = GtOrientationSampling.get_bcc_center_range(edge_points(1), oversize);
self.produce_R_vectors_simple_grid(sampling_ranges, dcenters, dranges, 2)
otherwise
error('GtOrientationSampling:make_simple_grid', ...
'Unknown lattice structure: %s', type)
end
if (self.verbose && usejava('jvm'))
self.plot_sampling(gvdm, dmean);
end
% Slow part which is based on gtCalculateGrain: could be improved!
self.produce_geometry();
self.produce_stats(type);
if (self.verbose && usejava('jvm'))
self.print_stats();
end
if (self.verbose && usejava('jvm'))
self.plot_blob_coverage();
end
end
function make_supersampling_simple_grid(self, directions, factor)
self.ss_factor = [1 1 1];
self.ss_factor(directions) = factor;
orient_voxel_size = self.stats.sampling.gaps;
half_voxel_size = orient_voxel_size ./ 2;
steps = orient_voxel_size ./ self.ss_factor;
half_steps = steps ./ 2;
for ii_o = 1:numel(self.lattice)
size_or = size(self.lattice(ii_o).gr);
self.lattice_ss(ii_o).gr = cell(size_or);
for ii_g = 1:numel(self.lattice(ii_o).gr)
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;
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', {[]} );
end
end
end
self.lattice_ss(ii_o).gr{ii_g} = gr_ss;
end
end
if (self.verbose && usejava('jvm'))
self.plot_supersampling();
end
self.produce_geometry_supersampling();
end
function valid = is_valid(self, type, edge_points)
valid = strcmpi(self.stats.sampling.type, type) ...
&& all(size(self.lattice(1).gr) == edge_points([1 1 1]));
end
function [or, or_ss] = get_orientations(self)
or = {};
for ii = 1:numel(self.lattice)
or = [or(:); self.lattice(ii).gr(:)];
end
or_ss = {};
for ii = 1:numel(self.lattice_ss)
or_ss = [or_ss(:); self.lattice_ss(ii).gr(:)];
end
end
function ref_gr = get_reference_grain(self)
ref_gr = self.ref_gr;
end
function r_vectors = get_R_vectors(self)
r_vectors = self.R_vectors;
end
function or_size = get_orientation_sampling_size(self)
num_lattices = numel(self.lattice);
or_size = cell(num_lattices, 1);
for ii = 1:num_lattices
or_size{ii} = [ ...
size(self.lattice(ii).gr, 1), ...
size(self.lattice(ii).gr, 2), ...
size(self.lattice(ii).gr, 3) ];
end
end
function tot_orient = get_total_num_orientations(self)
tot_orient = sum(arrayfun(@(x)numel(x.gr), self.lattice));
end
function [avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = getAverageOrientations(self, vols)
[avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = ...
GtOrientationSampling.computeAverageOrientations(self.get_orientations(), vols);
end
function avg_R_vec = getAverageOrientation(self, vols)
avg_R_vec = GtOrientationSampling.computeAverageOrientation(self.get_orientations(), vols);
end
function odf = getODF(self, vols)
odf = GtOrientationSampling.computeODF(self.get_orientations(), vols);
end
function [gvdm_tmp, verts, all_plane_normals] = guess_ODF_BB(self, use_corners, use_etas)
if (~exist('use_corners', 'var') || isempty(use_corners))
use_corners = true;
end
if (~exist('use_etas', 'var') || isempty(use_etas))
use_etas = false;
end
omega_step = gtGetOmegaStepDeg(self.parameters, self.detector_index);
samgeo = self.parameters.samgeo;
labgeo = self.parameters.labgeo;
detgeo = self.parameters.detgeo(self.detector_index);
bls = self.bl(self.selected);
g = gtMathsRod2OriMat(self.ref_gr.R_vector);
sel = self.ondet(self.included(self.selected));
num_refl = numel(sel);
if (~num_refl)
error('GtOrientationSampling:guess_ODF_BB:wrong_argument', ...
'No reflections selected! So, no estimation can be performed!')
end
% Taking theoretical omegas instead of experimental ones!
w = self.ref_gr.allblobs.omega(sel);
n = self.ref_gr.allblobs.eta(sel);
t = self.ref_gr.allblobs.theta(sel);
gc_sam = self.ref_gr.center(ones(numel(w), 1), :);
% can we determine the maximum eta spread?
uv = self.ref_gr.allblobs.detector(self.detector_index).uvw(sel, 1:2);
gc_lab = gtGeoSam2Lab(gc_sam, w, labgeo, samgeo, false);
% Computing the etas according to the corners of the blobs
u_min_max = cat(1, bls(:).mbbu);
v_min_max = cat(1, bls(:).mbbv);
pos_corner_bl = [u_min_max(:, 1) - 0.5, v_min_max(:, 1) - 0.5];
pos_corner_br = [u_min_max(:, 1) - 0.5, v_min_max(:, 2) + 0.5];
pos_corner_tl = [u_min_max(:, 2) + 0.5, v_min_max(:, 1) - 0.5];
pos_corner_tr = [u_min_max(:, 2) + 0.5, v_min_max(:, 2) + 0.5];
diffvec_corner_bl = gtGeoDet2Lab(pos_corner_bl, detgeo, false) - gc_lab;
diffvec_corner_br = gtGeoDet2Lab(pos_corner_br, detgeo, false) - gc_lab;
diffvec_corner_tl = gtGeoDet2Lab(pos_corner_tl, detgeo, false) - gc_lab;
diffvec_corner_tr = gtGeoDet2Lab(pos_corner_tr, detgeo, false) - gc_lab;
n_bl = gtGeoEtaFromDiffVec(diffvec_corner_bl, labgeo);
n_br = gtGeoEtaFromDiffVec(diffvec_corner_br, labgeo);
n_tl = gtGeoEtaFromDiffVec(diffvec_corner_tl, labgeo);
n_tr = gtGeoEtaFromDiffVec(diffvec_corner_tr, labgeo);
n_m_corners = min([n_bl, n_br, n_tl, n_tr], [], 2);
n_p_corners = max([n_bl, n_br, n_tl, n_tr], [], 2);
% Computing the etas according to the size of the blob
gc_det_pos = gtGeoGrainCenterSam2Det(self.ref_gr.center, w, self.parameters);
cosnc = abs(cosd(n)) >= 0.5;
dn = zeros(size(n));
bls_bbm = cat(1, bls(:).mbbsize);
% 1st component is the horizontal, while the 2nd is the
% vertical !!
dn(cosnc) = bls_bbm(cosnc, 1) ./ abs(cosd(n(cosnc)));
dn(~cosnc) = bls_bbm(~cosnc, 2) ./ abs(sind(n(~cosnc)));
dn = dn ./ (2 * sqrt(sum((uv - gc_det_pos) .^ 2, 2)));
dn = atand(dn);
% Let's take a small deviation in eta, to find the plane that
% determines the limits in orientation-space
n_m_symm = n - dn;
n_p_symm = n + dn;
if (use_corners)
n_m = n_m_corners;
n_p = n_p_corners;
else
n_m = n_m_symm;
n_p = n_p_symm;
end
% We retrieve the signed y (instead of the unsigned one):
% y = self.ref_gr.allblobs.plorig(sel, :)';
y0_exp_lab = gtGeoPlLabFromThetaEta(t, n, labgeo);
y0_exp_sam = gtGeoLab2Sam(y0_exp_lab, w, labgeo, samgeo, true);
y = y0_exp_sam';
% % Useful to validate the previously introduced flip of eta
% ominds = self.ref_gr.allblobs.omind(sel);
% ssp = ((ominds == 1) | (ominds == 2));
% ss12 = ssp - ~ssp;
%
% y_l = self.ref_gr.allblobs.pllab(sel, :);
% y_s = self.ref_gr.allblobs.pl(sel, :);
% y_o = self.ref_gr.allblobs.plorig(sel, :);
% [y0_exp_lab, y_l, y0_exp_lab-y_l]
% [y0_exp_sam, y_s, y0_exp_sam-y_s, ss12(:, [1 1 1]) .* y_o]
% Should be identical to pl_cry
hhs = gtVectorLab2Cryst(y', g);
ym_exp_lab = gtGeoPlLabFromThetaEta(t, n_m, labgeo);
yp_exp_lab = gtGeoPlLabFromThetaEta(t, n_p, labgeo);
% r0s related to the extreeme omegas
r0ws_ex = zeros(2 * num_refl, 3);
% r01s related to the extreeme omegas
r0ns_ex = zeros(2 * num_refl, 3);
% versors of r (classic) related to the extreeme omegas
rcs_ex_w = zeros(2 * num_refl, 3);
% versors of the perpendicular to the rs in the extreeme omegas
% (~ parallel to eta)
rns_ex_w = zeros(2 * num_refl, 3);
% versors of the perpendicular to the rs in the extreeme etas
% (~ parallel to onegas)
rws_ex_w = zeros(2 * num_refl, 3);
sub_ws_abs = cat(1, bls(:).mbbw);
% sub_w_abs will contain the two bounding omegas for the given
% reflection.
sub_ws_abs = (sub_ws_abs + repmat([-0.5, 0.5], [num_refl, 1])) * omega_step;
for ii_r = 1:num_refl
sub_w_abs = sub_ws_abs(ii_r, :);
% the y vectors for the central eta, at the bounding omegas
y0_exp_sam = gtGeoLab2Sam(y0_exp_lab([ii_r ii_r], :), ...
sub_w_abs, labgeo, samgeo, true);
% the y vectors for the lower eta, at the bounding omegas
ym_exp_sam = gtGeoLab2Sam(ym_exp_lab([ii_r ii_r], :), ...
sub_w_abs, labgeo, samgeo, true);
% the y vectors for the upper eta, at the bounding omegas
yp_exp_sam = gtGeoLab2Sam(yp_exp_lab([ii_r ii_r], :), ...
sub_w_abs, labgeo, samgeo, true);
hh = hhs(ii_r, :);
norms_0c = (1 + hh * y0_exp_sam')'; norms_0c = norms_0c(:, [1 1 1]);
norms_0m = (1 + hh * ym_exp_sam')';
norms_0m = norms_0m(:, [1 1 1]);
norms_0p = (1 + hh * yp_exp_sam')';
norms_0p = norms_0p(:, [1 1 1]);
% r0 is defined = (h x y) / (1 + h . y)
r0c = gtMathsCross(hh, y0_exp_sam) ./ norms_0c;
r0m = gtMathsCross(hh, ym_exp_sam) ./ norms_0m;
r0p = gtMathsCross(hh, yp_exp_sam) ./ norms_0p;
% r is defined = r0 + t * (h + y) / (1 + h . y)
r_dir_c = (hh([1 1], :) + y0_exp_sam) ./ norms_0c;
r_dir_m = (hh([1 1], :) + ym_exp_sam) ./ norms_0m;
r_dir_p = (hh([1 1], :) + yp_exp_sam) ./ norms_0p;
r_dir_c = gtMathsNormalizeVectorsList(r_dir_c);
r_dir_m = gtMathsNormalizeVectorsList(r_dir_m);
r_dir_p = gtMathsNormalizeVectorsList(r_dir_p);
% Translating r0s to the closest point to the average R
% vector, on the line r
dot_r0c_avgR = sum((self.ref_gr.R_vector([1 1], :) - r0c) .* r_dir_c, 2);
dot_r0m_avgR = sum((self.ref_gr.R_vector([1 1], :) - r0m) .* r_dir_m, 2);
dot_r0p_avgR = sum((self.ref_gr.R_vector([1 1], :) - r0p) .* r_dir_p, 2);
r0c = r0c + r_dir_c .* dot_r0c_avgR(:, [1 1 1]);
r0m = r0m + r_dir_m .* dot_r0m_avgR(:, [1 1 1]);
r0p = r0p + r_dir_p .* dot_r0p_avgR(:, [1 1 1]);
% Delimiters in w
r0ws_ex(ii_r*2-1:ii_r*2, :) = r0c;
% Delimiters in n
r0ns_ex(ii_r*2-1:ii_r*2, :) = [sum(r0m, 1); sum(r0p, 1)] / 2;
rcs_ex_w(ii_r*2-1:ii_r*2, :) = r_dir_c;
% It will be used for recentering the r0w, on top of the
% average point and get a colinear plane norm to the
% distance of the plane to the average point
r_dir_n = r0p - r0m;
r_dir_n_norm = sqrt(sum(r_dir_n .^ 2, 2));
rns_ex_w(ii_r*2-1:ii_r*2, :) = r_dir_n ./ r_dir_n_norm(:, [1 1 1]);
% Same as the previous, but for the points r0n
r_dir_w = [r0m(end, :) - r0m(1, :); r0p(end, :) - r0p(1, :)];
r_dir_w_norm = sqrt(sum(r_dir_w .^ 2, 2));
rws_ex_w(ii_r*2-1:ii_r*2, :) = r_dir_w ./ r_dir_w_norm(:, [1 1 1]);
end
% Let's find the closest point to the average R vector on the plane
avg_R_vec_exp = self.ref_gr.R_vector(ones(size(r0ws_ex, 1), 1), :);
trasl_ex_points_w = avg_R_vec_exp - r0ws_ex;
trasl_ex_points_w = sum(trasl_ex_points_w .* rns_ex_w, 2) ./ sqrt(sum(rns_ex_w .^ 2, 2));
trasl_ex_points_w = r0ws_ex + trasl_ex_points_w(:, [1 1 1]) .* rns_ex_w;
if (use_etas)
trasl_ex_points_n = avg_R_vec_exp - r0ns_ex;
trasl_ex_points_n = sum(trasl_ex_points_n .* rws_ex_w, 2) ./ sqrt(sum(rws_ex_w .^ 2, 2));
trasl_ex_points_n = r0ns_ex + trasl_ex_points_n(:, [1 1 1]) .* rws_ex_w;
all_plane_normals = [ ...
trasl_ex_points_w - avg_R_vec_exp; ...
trasl_ex_points_n - avg_R_vec_exp ];
else
all_plane_normals = trasl_ex_points_w - avg_R_vec_exp;
end
% Computing the circumscribing polyhedron in Orientation-space
verts = gtMathsGetPolyhedronVerticesFromPlaneNormals(all_plane_normals);
verts = verts + self.ref_gr.R_vector(ones(size(verts, 1), 1), :);
% f = figure();
% ax = axes('parent', f);
% hold(ax, 'on')
% scatter3(ax, trasl_ex_points_w(:, 1), trasl_ex_points_w(:, 2), trasl_ex_points_w(:, 3), 'b')
% if (use_etas)
% scatter3(ax, trasl_ex_points_n(:, 1), trasl_ex_points_n(:, 2), trasl_ex_points_n(:, 3), 'y')
% end
% scatter3(ax, verts(:, 1), verts(:, 2), verts(:, 3), 'r')
% hold(ax, 'off')
% bbox_r0s = [min(trasl_ex_points, [], 1), max(trasl_ex_points, [], 1)];
% bbox_r0s = [bbox_r0s(1:3), bbox_r0s(4:6) - bbox_r0s(1:3)];
% or_bbox = bbox_r0s;
bbox_verts = [min(verts, [], 1), max(verts, [], 1)];
bbox_verts = [bbox_verts(1:3), bbox_verts(4:6) - bbox_verts(1:3)];
or_bbox = bbox_verts;
gvdm_tmp = [];
for cz = 0:1
for cy = 0:1
for cx = 0:1
gvdm_tmp = [gvdm_tmp, (or_bbox(1:3)' + [cx; cy; cz] .* or_bbox(4:6)')]; %#ok<AGROW>
end
end
end
end
function resolution = estimate_maximum_resolution(self)
sel = self.ondet(self.included(self.selected));
n = self.ref_gr.allblobs.eta(sel);
t = self.ref_gr.allblobs.theta(sel);
detgeo = self.parameters.detgeo(self.detector_index);
delta_pixel = abs(cosd(n)) .* detgeo.pixelsizeu + abs(sind(n)) .* detgeo.pixelsizev;
xy_resolution = atand(delta_pixel ./ (norm(detgeo.detrefpos) .* tand(2 .* t)));
xy_resolution = min(xy_resolution);
z_resolution = gtGetOmegaStepDeg(self.parameters, self.detector_index);
resolution = [xy_resolution, xy_resolution, z_resolution];
end
function num_interp = suggest_num_interp(self, z_res)
if (~exist('z_res', 'var') || isempty(z_res))
z_res = 2 * atand(self.stats.sampling.gaps(3));
end
om_step = gtGetOmegaStepDeg(self.parameters);
% The factor 1.5 allows at least a minimum amount of crosstalk among
% the orientations
num_interp = max(z_res / om_step * 1.5, 1); % / self.ss_factor(3)
end
function num_interp = minimum_num_interp(self, z_res)
if (~exist('z_res', 'var') || isempty(z_res))
z_res = 2 * atand(self.stats.sampling.gaps(3));
end
om_step = gtGetOmegaStepDeg(self.parameters);
num_interp = max(z_res / om_step, 1); % / self.ss_factor(3)
end
function sampler = regenerate_displaced(self, new_center)
detgeo = self.parameters.detgeo;
labgeo = self.parameters.labgeo;
samgeo = self.parameters.samgeo;
sampler = self;
% Let's update the reference grain sampler.ref_gr.center = new_center;
rot_s2l = permute(sampler.ref_gr.allblobs.srot, [2 1 3]);
gcsam_v = new_center(ones(1, size(rot_s2l, 3)), :);
gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
num_dets = numel(detgeo);
for n = 1:num_dets
sampler.ref_gr.allblobs.detector(n).uvw(:, 1:2) = gtFedPredictUVMultiple( ...
[], sampler.ref_gr.allblobs.dvec', gclab_v', ...
detgeo(n).detrefpos', detgeo(n).detnorm', detgeo(n).Qdet, ...
[detgeo(n).detrefu, detgeo(n).detrefv]')';
end
% And now let's take care of the sampled orientations
grid_gr = self.get_orientations();
grid_gr = [grid_gr{:}];
num_orients = numel(grid_gr);
ab = [grid_gr(:).allblobs];
rot_l2s = cat(3, ab(:).srot);
rot_s2l = permute(rot_l2s, [2 1 3]);
dveclab = cat(1, ab(:).dvec);
gcsam_v = new_center(ones(1, size(rot_s2l, 3)), :);
gclab_v = gtGeoSam2Lab(gcsam_v, rot_s2l, labgeo, samgeo, false);
uv_tot = cell(num_dets, 1);
for n = 1:num_dets
uv_tot{n} = gtFedPredictUVMultiple([], dveclab', gclab_v', ...
detgeo(n).detrefpos', detgeo(n).detnorm', detgeo(n).Qdet, ...
[detgeo(n).detrefu, detgeo(n).detrefv]');
uv_tot{n} = reshape(uv_tot{n}, 2, [], num_orients);
uv_tot{n} = permute(uv_tot{n}, [2 1 3]);
end
for ii_g = 1:num_orients
sampler.lattice(1).gr{ii_g}.center = new_center;
for n = 1:num_dets
sampler.lattice(1).gr{ii_g}.allblobs.detector(n).uvw(:, 1:2) = uv_tot{n}(:, :, ii_g);
end
end
end
function [delta_omegas, proj_lims] = get_omega_deviations(self, with_shape_functions)
if (~exist('with_shape_functions', 'var') || isempty(with_shape_functions))
with_shape_functions = false;
end
if (with_shape_functions)
sub_space_size = self.stats.sampling.gaps;
half_rspace_sizes = sub_space_size / 2;
min_bounds = self.lattice(1).gr{1, 1, 1}.R_vector - half_rspace_sizes;
max_bounds = self.lattice(1).gr{end, end, end}.R_vector + half_rspace_sizes;
bounds_x = [min_bounds(1), max_bounds(1)];
bounds_y = [min_bounds(2), max_bounds(2)];
bounds_z = [min_bounds(3), max_bounds(3)];
grs = cell(2, 2, 2);
for dx = 1:2
for dy = 1:2
for dz = 1:2
r_vec = [bounds_x(dx), bounds_y(dy), bounds_z(dz)];
grs{dx, dy, dz} = struct( ...
'phaseid', self.ref_gr.phaseid, ...
'center', self.ref_gr.center, 'R_vector', r_vec ); end
end
end
grs = gtCalculateGrain_p(grs, self.parameters, ...
'ref_omind', self.ref_gr.allblobs.omind, ...
'included', self.ondet(self.included(self.selected)) );
grs = [grs{:}];
else
grs = [self.lattice(1).gr{:}];
end
ab = [grs(:).allblobs];
w_tab = [ab(:).omega] / gtGetOmegaStepDeg(self.parameters, self.detector_index);
sel = self.ondet(self.included(self.selected));
ref_int_ws = self.ref_gr.allblobs.detector(self.detector_index).uvw(sel, 3);
w_tab = gtGrainAnglesTabularFix360deg(w_tab, ref_int_ws, self.parameters);
max_ws = ceil(max(w_tab, [], 2));
min_ws = floor(min(w_tab, [], 2));
proj_lims = [min_ws, max_ws];
delta_omegas = max_ws - min_ws + 1;
end
function [delta_uvs, proj_uv_lims] = get_uv_deviations(self)
grs = [self.lattice(1).gr{:}];
num_blobs = numel(find(self.selected));
num_oris = numel(grs);
uv_tab = zeros(num_blobs, 2, num_oris);
for ii_o = 1:num_oris
uv_tab(:, :, ii_o) = grs(ii_o).allblobs.detector(self.detector_index).uvw(:, 1:2);
end
min_proj_uv = reshape(min(floor(uv_tab), [], 3), [], 1, 2);
max_proj_uv = reshape(max(ceil(uv_tab), [], 3), [], 1, 2);
proj_uv_lims = [min_proj_uv, max_proj_uv];
delta_uvs = reshape(proj_uv_lims(:, 2, :) - proj_uv_lims(:, 1, :), [], 2);
end
function bytes = get_memory_consumption_geometry(self)
bytes = GtBenchmarks.getSizeVariable(self.lattice) ...
+ GtBenchmarks.getSizeVariable(self.lattice_ss) ...
+ GtBenchmarks.getSizeVariable(self.R_vectors);
end
function bytes = get_memory_consumption_graindata(self)
bytes = GtBenchmarks.getSizeVariable(self.ref_gr);
end
end
methods (Access = protected)
function produce_R_vectors_simple_grid(self, sampling_range, dcenters, dranges, o_ii)
if (~exist('o_ii', 'var'))
o_ii = 1;
end
num_comps = numel(sampling_range);
try
self.R_vectors{o_ii} = zeros(num_comps ^ 3, 3);
catch mexc
fprintf('Trying to allocate a [%d, 3] array.. but without enough memory!\n', ...
num_comps ^ 3)
rethrow(mexc)
end
points = zeros(num_comps, 3);
dcenters = dcenters(ones(1, num_comps), :);
dranges = dranges(ones(1, num_comps), :);
self.lattice(o_ii).gr = cell(num_comps, num_comps, num_comps);
l_ii = 0;
for z_ii = 1:num_comps
for y_ii = 1:num_comps
points(:, 1) = sampling_range(:);
points(:, 2) = sampling_range(y_ii);
points(:, 3) = sampling_range(z_ii);
R_vecs = dcenters + points .* dranges;
self.R_vectors{o_ii}(l_ii+1:l_ii+num_comps, :) = R_vecs;
l_ii = l_ii + num_comps;
end
end
end
function produce_R_vectors_even_simple_grid(self, sampling_ranges, dcenters, dranges, o_ii)
if (~exist('o_ii', 'var'))
o_ii = 1;
end
num_z_comps = numel(sampling_ranges{3});
num_y_comps = numel(sampling_ranges{2});
num_x_comps = numel(sampling_ranges{1});
try
self.R_vectors{o_ii} = zeros(num_z_comps * num_y_comps * num_x_comps, 3);
catch mexc
fprintf('Trying to allocate a [%d, 3] array.. but without enough memory!\n', ...
num_z_comps * num_y_comps * num_x_comps)
rethrow(mexc)
end
points = zeros(num_x_comps, 3);
dcenters = dcenters(ones(1, num_x_comps), :);
dranges = dranges(ones(1, num_x_comps), :);
self.lattice(o_ii).gr = cell(num_x_comps, num_y_comps, num_z_comps);
l_ii = 0;
for z_ii = 1:num_z_comps
for y_ii = 1:num_y_comps
points(:, 1) = sampling_ranges{1}(:);
points(:, 2) = sampling_ranges{2}(y_ii);
points(:, 3) = sampling_ranges{3}(z_ii);
R_vecs = dcenters + points .* dranges;
self.R_vectors{o_ii}(l_ii+1:l_ii+num_x_comps, :) = R_vecs;
l_ii = l_ii + num_x_comps;
end
end
end
function produce_geometry(self)
grain_center = self.get_grain_center();
ref_omind = self.ref_gr.allblobs.omind;
ref_included = self.ondet(self.included(self.selected));
fprintf('Computing geometries: ')
c = tic();
for o_ii = 1:numel(self.lattice)
num_orients = numel(self.lattice(o_ii).gr);
num_chars_o = fprintf('order %d/%d: (%d) ', o_ii-1, numel(self.lattice), num_orients);
order_gvdm = self.R_vectors{o_ii};
num_chars_gr = fprintf('creating initial structures..');
for ii = 1:num_orients
self.lattice(o_ii).gr{ii} = struct( ...
'R_vector', order_gvdm(ii, :), ...
'center', grain_center, ...
'phaseid', self.ref_gr.phaseid ); end
fprintf(repmat('\b', [1 num_chars_gr]));
fprintf(repmat(' ', [1 num_chars_gr]));
fprintf(repmat('\b', [1 num_chars_gr]));
switch(2)
case 3 % All in once
num_chars_gr = fprintf('gtCalculateGrain in parallel..');
self.lattice(o_ii).gr = gtCalculateGrain_p( ...
self.lattice(o_ii).gr, self.parameters, ...
'ref_omind', ref_omind, 'included', ref_included );
case 2 % In chunks
num_chars_gr = fprintf('gtCalculateGrain in chunks: ');
chunk_size = 1000;
for ii = 1:chunk_size:num_orients
end_chunk = min(ii+chunk_size, num_orients);
num_chars_gr_ii = fprintf('%05d/%05d', ii, num_orients);
self.lattice(o_ii).gr(ii:end_chunk) = gtCalculateGrain_p( ...
self.lattice(o_ii).gr(ii:end_chunk), self.parameters, ...
'ref_omind', ref_omind, 'included', ref_included );
fprintf(repmat('\b', [1 num_chars_gr_ii]));
end
fprintf(repmat(' ', [1 num_chars_gr_ii]));
fprintf(repmat('\b', [1 num_chars_gr_ii]));
case 1 % One by One
num_chars_gr = fprintf('gtCalculateGrain: ');
for ii = 1:num_orients
num_chars_gr_ii = fprintf('%05d/%05d', ii, num_orients);
self.lattice(o_ii).gr{ii} = gtCalculateGrain( ...
self.lattice(o_ii).gr{ii}, self.parameters, ...
'ref_omind', ref_omind, 'included', ref_included );
fprintf(repmat('\b', [1 num_chars_gr_ii]));
end
fprintf(repmat(' ', [1 num_chars_gr_ii]));
fprintf(repmat('\b', [1 num_chars_gr_ii]));
end
fprintf(repmat('\b', [1 num_chars_gr]));
fprintf(repmat(' ', [1 num_chars_gr]));
fprintf(repmat('\b', [1 num_chars_gr]));
fprintf(repmat('\b', [1 num_chars_o]));
end
total_orient = sum(arrayfun(@(x)numel(x.gr), self.lattice));
fprintf(' (orders: %d, total orientations: %04d, in %3.1f s) Done.\n', ...
numel(self.lattice), total_orient, toc(c));
end
function produce_geometry_supersampling(self)
ref_omind = self.ref_gr.allblobs.omind;
ref_included = self.ondet(self.included(self.selected));
fprintf('Computing super-sampling geometries: ')
c = tic();
for o_ii = 1:numel(self.lattice_ss)
num_orients = numel(self.lattice_ss(o_ii).gr);
num_chars_o = fprintf('order %d/%d: (%d) orientation: ', o_ii-1, numel(self.lattice_ss), num_orients);
for g_ii = 1:num_orients
num_chars_gr_ii = fprintf('%05d/%05d', g_ii, num_orients);
self.lattice_ss(o_ii).gr{g_ii} = gtCalculateGrain_p( ...
self.lattice_ss(o_ii).gr{g_ii}, self.parameters, ...
'ref_omind', ref_omind, 'included', ref_included );
fprintf(repmat('\b', [1 num_chars_gr_ii]));
end
fprintf(repmat('\b', [1 num_chars_o]));
end
total_orient = sum(arrayfun(@(x)numel(x.gr), self.lattice_ss)); fprintf(' (orders: %d, total orientations: %04d, in %3.1f s) Done.\n', ...
numel(self.lattice), total_orient, toc(c));
end
function produce_stats(self, type)
self.stats.sampling.type = type;
switch(type)
case 'cubic'
% self.compute_pixel_distances_simple_grid();
self.compute_blob_coverage_simple_grid();
self.compute_sampling_distances_cubic();
case 'fcc'
case 'bcc'
% self.compute_pixel_distances_bcc();
self.compute_blob_coverage_simple_grid();
self.compute_sampling_distances_bcc();
end
end
function compute_sampling_distances_cubic(self)
self.stats.sampling.gaps(1) = self.lattice.gr{2, 1, 1}.R_vector(1) - self.lattice.gr{1, 1, 1}.R_vector(1);
self.stats.sampling.gaps(2) = self.lattice.gr{1, 2, 1}.R_vector(2) - self.lattice.gr{1, 1, 1}.R_vector(2);
self.stats.sampling.gaps(3) = self.lattice.gr{1, 1, 2}.R_vector(3) - self.lattice.gr{1, 1, 1}.R_vector(3);
self.stats.sampling.bbox_size = self.lattice.gr{end}.R_vector - self.lattice.gr{1}.R_vector;
self.stats.sampling.bbox_shift = self.lattice.gr{1}.R_vector;
end
function compute_sampling_distances_bcc(self)
self.stats.sampling.gaps = self.lattice(2).gr{1}.R_vector - self.lattice(1).gr{1}.R_vector;
self.stats.sampling.bbox_size = self.lattice(1).gr{end}.R_vector - self.lattice(1).gr{1}.R_vector;
self.stats.sampling.bbox_shift = self.lattice(1).gr{1}.R_vector;
end
function compute_pixel_distances_simple_grid(self)
max_dist_u = 0;
max_dist_v = 0;
max_angle_u = 0;
max_angle_v = 0;
for z_ii = 1:size(self.lattice(1).gr, 3)
for y_ii = 1:size(self.lattice(1).gr, 2)
for x_ii = 1:size(self.lattice(1).gr, 1)-1
gr1 = self.lattice(1).gr{x_ii, y_ii, z_ii};
gr2 = self.lattice(1).gr{x_ii+1, y_ii, z_ii};
[gr12_dist, gr12_angle] = self.compute_max_proj_pixel_dist(gr1, gr2);
max_dist_u = max(max_dist_u, gr12_dist);
max_angle_u = max(max_angle_u, gr12_angle);
end
end
for y_ii = 1:size(self.lattice(1).gr, 2)-1
for x_ii = 1:size(self.lattice(1).gr, 1)
gr1 = self.lattice(1).gr{x_ii, y_ii, z_ii};
gr2 = self.lattice(1).gr{x_ii, y_ii+1, z_ii};
[gr12_dist, gr12_angle] = self.compute_max_proj_pixel_dist(gr1, gr2);
max_dist_v = max(max_dist_v, gr12_dist);
max_angle_v = max(max_angle_v, gr12_angle);
end
end
end
self.stats.proj.max_pixel_dist_u = max_dist_u;
self.stats.proj.max_pixel_dist_v = max_dist_v;
self.stats.proj.max_angle_u = max_angle_u;
self.stats.proj.max_angle_v = max_angle_v;
self.stats.sample.avg_dist = self.compute_average_sample_distance();
end
function compute_pixel_distances_bcc(self)
max_dist_x = 0;
max_dist_y = 0;
max_angle_x = 0;
max_angle_y = 0;
for z_ii = 1:size(self.lattice(2).gr, 3)
for y_ii = 1:size(self.lattice(2).gr, 2)
for x_ii = 1:size(self.lattice(2).gr, 1)
gr1 = self.lattice(1).gr{x_ii, y_ii, z_ii};
gr2 = self.lattice(2).gr{x_ii, y_ii, z_ii};
[gr12_dist, gr12_angle] = self.compute_max_proj_pixel_dist(gr1, gr2);
max_dist_x = max(max_dist_x, gr12_dist);
max_angle_x = max(max_angle_x, gr12_angle);
end
end
for y_ii = 1:size(self.lattice(2).gr, 2)
for x_ii = 1:size(self.lattice(2).gr, 1)
gr1 = self.lattice(1).gr{x_ii+1, y_ii, z_ii};
gr2 = self.lattice(2).gr{x_ii, y_ii, z_ii};
[gr12_dist, gr12_angle] = self.compute_max_proj_pixel_dist(gr1, gr2);
max_dist_x = max(max_dist_x, gr12_dist);
max_angle_x = max(max_angle_x, gr12_angle);
end
end
for y_ii = 1:size(self.lattice(2).gr, 2)
for x_ii = 1:size(self.lattice(2).gr, 1)
gr1 = self.lattice(1).gr{x_ii, y_ii+1, z_ii};
gr2 = self.lattice(2).gr{x_ii, y_ii, z_ii};
[gr12_dist, gr12_angle] = self.compute_max_proj_pixel_dist(gr1, gr2);
max_dist_x = max(max_dist_x, gr12_dist);
max_angle_x = max(max_angle_x, gr12_angle);
end
end
for y_ii = 1:size(self.lattice(2).gr, 2)
for x_ii = 1:size(self.lattice(2).gr, 1)
gr1 = self.lattice(1).gr{x_ii+1, y_ii+1, z_ii};
gr2 = self.lattice(2).gr{x_ii, y_ii, z_ii};
[gr12_dist, gr12_angle] = self.compute_max_proj_pixel_dist(gr1, gr2);
max_dist_x = max(max_dist_x, gr12_dist);
max_angle_x = max(max_angle_x, gr12_angle);
end
end
end
self.stats.proj.max_pixel_dist_u = max_dist_x;
self.stats.proj.max_pixel_dist_v = max_dist_y;
self.stats.proj.max_angle_u = max_angle_x;
self.stats.proj.max_angle_v = max_angle_y;
self.stats.sample.avg_dist = self.compute_average_sample_distance();
end
function dist = compute_average_sample_distance(self)
dist = 0;
tot_orient = get_total_num_orientations(self);
for o_ii = 1:numel(self.lattice)
for z_ii = 1:size(self.lattice(o_ii).gr, 3)
for y_ii = 1:size(self.lattice(o_ii).gr, 2)
for x_ii = 1:size(self.lattice(o_ii).gr, 1)
ab = self.lattice(o_ii).gr{x_ii, y_ii, z_ii}.allblobs;
if (isfield(ab, 'detector'))
positions = ab.detector(1).proj_geom(:, 4:6);
else
positions = ab.proj_geom(:, 4:6);
end dist = dist + sum(sqrt(sum(positions .^ 2, 2))) / size(positions, 1);
end
end
end
end
dist = dist / tot_orient;
end
function [dist, angle] = compute_max_proj_pixel_dist(self, gr1, gr2)
if (isfield(gr1.allblobs, 'detector'))
geom1 = gr1.allblobs.detector(1).proj_geom;
geom2 = gr2.allblobs.detector(1).proj_geom;
else
geom1 = gr1.allblobs.proj_geom;
geom2 = gr2.allblobs.proj_geom;
end
dirs1 = geom1(:, 1:3);
dirs2 = geom2(:, 1:3);
det_pos1 = geom1(:, 4:6);
det_pos2 = geom2(:, 4:6);
lengths1 = sqrt(sum(det_pos1 .^ 2, 2));
lengths2 = sqrt(sum(det_pos2 .^ 2, 2));
disps1 = dirs1 .* lengths1(:, [1 1 1]) - det_pos1;
disps2 = dirs2 .* lengths2(:, [1 1 1]) - det_pos2;
% Approximately
disps = disps1 - disps2;
[dist, ~] = max(sqrt(sum(disps .^ 2, 2)));
% rot_dirs = cross(det_pos2, det_pos1);
% norm_rot_dirs = sqrt(sum(rot_dirs .^ 2, 2));
% rot_dirs = rot_dirs ./ norm_rot_dirs(:, [1 1 1]);
%
% renorm_det_pos = lengths1 .* lengths2;
% dot_det_pos = det_pos1 .* det_pos2;
% rot_angles = abs(acos(sum(dot_det_pos, 2) ./ renorm_det_pos));
%
% rotated_dirs2 = zeros(size(dirs2));
%
% for r_ii = 1:numel(rot_angles)
% rot_comp = gtMathsRotationMatrixComp(rot_dirs(r_ii, :), 'row');
% rot_tensor = gtMathsRotationTensor(rot_angles(r_ii), rot_comp);
%
% rotated_dirs2(r_ii, :) = dirs2(r_ii, :) * rot_tensor;
% end
%
% angles = acos(sum(dirs1 .* rotated_dirs2, 2));
angles = abs(acos(sum(dirs1 .* dirs2, 2)));
[angle, ~] = max(angles);
end
function compute_blob_coverage_simple_grid(self)
omega_step = self.get_omega_step_deg();
bls = self.bl(self.selected);
num_blobs = numel(bls);
max_depth_blobs = 0;
for b_ii = 1:num_blobs
max_depth_blobs = max(max_depth_blobs, size(bls(b_ii).intm, 3));
end
switch(self.stats.sampling.type)
case 'cubic'
num_z_comps = size(self.lattice.gr, 3);
% Coverage of the slices in blobs (by Z layer of the sampling)
self.stats.pics.blob_coverage_map = cell(num_z_comps);
end
% Total intensity of each slice in a blob
self.stats.pics.blob_intensity = zeros(num_blobs, max_depth_blobs);
% Coverage of the slices in blobs (generic)
self.stats.pics.blob_coverage = zeros(num_blobs, max_depth_blobs);
min_blob_pos = zeros(num_blobs, 1);
max_blob_pos = zeros(num_blobs, 1);
for b_ii = 1:num_blobs
depth_blob = size(bls(b_ii).intm, 3);
min_blob_pos(b_ii) = floor((max_depth_blobs - depth_blob)/2) + 1;
max_blob_pos(b_ii) = min_blob_pos(b_ii) + depth_blob - 1;
self.stats.pics.blob_intensity(b_ii, min_blob_pos(b_ii):max_blob_pos(b_ii)) ...
= squeeze(sum(sum(bls(b_ii).intm, 2), 1));
end
limits = cat(1, bls(:).bbwim);
for o_ii = 1:numel(self.lattice)
for z_ii = 1:size(self.lattice(o_ii).gr, 3)
switch(self.stats.sampling.type)
case 'cubic'
self.stats.pics.blob_coverage_map{z_ii} = zeros(num_blobs, max_depth_blobs);
end
for y_ii = 1:size(self.lattice(o_ii).gr, 2)
for x_ii = 1:size(self.lattice(o_ii).gr, 1)
ab = self.lattice(o_ii).gr{x_ii, y_ii, z_ii}.allblobs;
omegas = ab.omega / omega_step;
within_lims = find(omegas >= limits(:, 1) & omegas <= limits(:, 2));
pos_in_blobs = round(omegas(within_lims)) - limits(within_lims, 1) + 1;
pos_in_coverage = pos_in_blobs + min_blob_pos(within_lims) - 1;
indx = sub2ind(size(self.stats.pics.blob_coverage), within_lims, pos_in_coverage);
switch(self.stats.sampling.type)
case 'cubic'
self.stats.pics.blob_coverage_map{z_ii}(indx) = ...
self.stats.pics.blob_coverage_map{z_ii}(indx) + 1;
end
self.stats.pics.blob_coverage(indx) = 1;
end
end
end
end
% Intensity reached by at least one orientation
self.stats.pics.blob_covered_intensity = ...
self.stats.pics.blob_intensity .* self.stats.pics.blob_coverage;
% Intensity NOT reached by any orientation
self.stats.pics.blob_uncovered_intensity = ...
self.stats.pics.blob_intensity .* (1 - self.stats.pics.blob_coverage);
end
function plot_blob_coverage(self)
f1 = figure();
lims = [min(self.stats.pics.blob_intensity(:)) max(self.stats.pics.blob_intensity(:))];
ax1 = subplot(1, 4, 1, 'parent', f1);
imagesc(self.stats.pics.blob_covered_intensity, 'parent', ax1, lims);
title(ax1, 'Covered Intensity');
ax2 = subplot(1, 4, 2, 'parent', f1);
imagesc(self.stats.pics.blob_uncovered_intensity, 'parent', ax2, lims);
title(ax2, 'Uncovered Intensity');
ax2 = subplot(1, 4, 3, 'parent', f1); imagesc(self.stats.pics.blob_coverage, 'parent', ax2);
title(ax2, 'Coverage');
ax2 = subplot(1, 4, 4, 'parent', f1);
imagesc(self.stats.pics.blob_intensity, 'parent', ax2);
title(ax2, 'Intensity');
switch(self.stats.sampling.type)
case 'cubic'
num_z_comps = size(self.lattice.gr, 3);
% Let's visualize the maps separately
grid_rows = floor(sqrt(num_z_comps));
grid_cols = ceil(num_z_comps / grid_rows);
f2 = figure();
for ii = 1:num_z_comps
ax = subplot(grid_rows, grid_cols, ii, 'parent', f2);
imagesc(self.stats.pics.blob_coverage_map{ii}, 'parent', ax);
title(ax, sprintf('Coverage of Z Layer %d', ii));
end
otherwise
end
drawnow();
end
function plot_sampling(self, gvdm, dmean)
points = self.R_vectors;
points = cat(1, points{:});
figure %, scatter3(points(:, 1), points(:, 2), points(:, 3))
hold on
scatter3(gvdm(1, :)', gvdm(2, :)', gvdm(3, :)', 20, 'c')
scatter3(points(:, 1), points(:, 2), points(:, 3), 20, 'r')
scatter3(mean(points(:, 1)), mean(points(:, 2)), mean(points(:, 3)), 30, 'm', 'filled')
scatter3(dmean(1), dmean(2), dmean(3), 30, 'g', 'filled')
hold off
drawnow();
end
function plot_supersampling(self)
points = [self.lattice.gr{:}];
points = cat(1, points.R_vector);
points_os = [self.lattice_ss.gr{:}];
points_os = [points_os{:}];
points_os = cat(1, points_os.R_vector);
figure %, scatter3(points(:, 1), points(:, 2), points(:, 3))
hold on
scatter3(points(:, 1), points(:, 2), points(:, 3), 20, 'c')
scatter3(points_os(:, 1), points_os(:, 2), points_os(:, 3), 20, 'r')
hold off
drawnow();
end
function print_sampling_details(self, dranges, oversize, edge_points)
if (numel(edge_points) == 1)
edge_points = edge_points([1 1 1]);
end
size_bb = 2 * atand(dranges * oversize);
max_resolution = self.estimate_maximum_resolution();
max_res_sampling_points = ceil(size_bb ./ max_resolution) + 1;
current_resolution = size_bb ./ (edge_points-1);
num_interp_range = current_resolution(3) / max_resolution(3);
num_interp_range = [num_interp_range, num_interp_range * 2];
detected_num_interp = self.parameters.rec.grains.options.num_interp; if (detected_num_interp < 0)
detected_num_interp = abs(detected_num_interp) * self.minimum_num_interp(current_resolution(3));
elseif (~detected_num_interp)
detected_num_interp = self.suggest_num_interp(current_resolution(3));
end
fprintf('\nOrientation-sapce sampling details:\n')
fprintf('- Bounding-box size (deg): [%g, %g, %g]\n', size_bb)
fprintf('- Maximum estimated resolution (deg): [%g, %g, %g] (Sampling points: [%d, %d, %d] => %d)\n', ...
max_resolution, max_res_sampling_points, prod(max_res_sampling_points));
fprintf('- Resolution (deg): [%g, %g, %g] (Sampling points: [%d, %d, %d] => %d)\n', ...
current_resolution, edge_points, prod(edge_points));
fprintf('- Suggested num_interp interval: [%g, %g], detected: %g\n\n', ...
num_interp_range, detected_num_interp);
end
function print_stats(self)
fprintf('Sampling:\n - lattice type: "%s" (size: %d, %d, %d)\n', ...
self.stats.sampling.type, size(self.lattice(1).gr));
fprintf(' - Gaps: x %e, y %e, z %e\n', ...
self.stats.sampling.gaps(1), ...
self.stats.sampling.gaps(2), ...
self.stats.sampling.gaps(3) );
fprintf(' - BBox: size (%e, %e, %e), shift (%e, %e, %e)\n', ...
self.stats.sampling.bbox_size, ...
self.stats.sampling.bbox_shift )
fprintf('Average smple distance: %g (pixels)\n', ...
self.stats.sample.avg_dist)
fprintf('Max projection pixel distances: x %f, y %f\n', ...
self.stats.proj.max_pixel_dist_u, ...
self.stats.proj.max_pixel_dist_v );
fprintf('Max projection angle distances: x %e (deg %f), y %e (deg %f)\n', ...
self.stats.proj.max_angle_u, ...
(180/pi*self.stats.proj.max_angle_u), ...
self.stats.proj.max_angle_v, ...
(180/pi*self.stats.proj.max_angle_v) );
end
function omega_step = get_omega_step_deg(self)
omega_step = 180 / self.parameters.acq.nproj;
end
function omega_step = get_omega_step_rad(self)
omega_step = pi / self.parameters.acq.nproj;
end
function gc = get_grain_center(self)
if (~isempty(self.fedpars))
gc = self.fedpars.vox000sam' + ...
(self.fedpars.grainsize' ./ 2 + [0.5; 0.5; 0.5]) .* self.parameters.recgeo.voxsize';
gc = gc';
else
gc = self.ref_gr.center; % .* self.parameters.recgeo.voxsize;
end
end
function [raw_depths, weighted_depths] = get_blob_depths(self)
% Try to account for Forentz factor, at least to avoid that max
% blows up. Computation is only approximated, so doesn't really
% count
indx = self.ondet(self.included(self.selected));
etas = self.ref_gr.allblobs.eta(indx);
factors = abs(sind(etas));
raw_depths = cellfun(@(x){sum(sum(abs(x), 1), 2)}, {self.bl(self.selected).intm}');
raw_depths = cellfun(@(x)(find(x >= 0, 1, 'last') - find(x >= 0, 1, 'first') + 1), raw_depths);
% raw_depths = cellfun(@(x)size(x, 3), {self.bl(self.selected).intm}') - 2;
weighted_depths = raw_depths .* factors;
end
function [max_blob_depth, min_blob_depth] = get_range_blob_depth(self)
max_blob_depth = 0;
min_blob_depth = +Inf;
num_blobs = numel(self.bl);
etas = self.ref_gr.allblobs.eta(self.ondet(self.included(self.selected)));
for b_ii = 1:num_blobs
% 2 is the pading we apply in forward simulation
current_blob_depth = (size(self.bl(b_ii).intm, 3) - 2);
lorentz_factored_depth = current_blob_depth * abs(sind(etas(b_ii)));
max_blob_depth = max(max_blob_depth, lorentz_factored_depth);
min_blob_depth = min(min_blob_depth, current_blob_depth);
end
end
end
methods (Access = protected, Static)
function range = get_cubic_range(edge_num, oversize)
if (~exist('oversize', 'var'))
oversize = 1.1;
end
limit = (edge_num - 1) / 2;
range = (-limit:limit) / (limit * 2) * oversize;
end
function range = get_bcc_center_range(edge_num, oversize)
if (~exist('oversize', 'var'))
oversize = 1.1;
end
limit_bcc = (edge_num - 1) / 2;
limit = (edge_num - 2) / 2;
range = (-limit:limit) / (limit_bcc * 2) * oversize;
end
end
methods (Access = public, Static)
function [avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = computeAverageOrientations(ors, or_vols, selected)
if (~exist('selected', 'var') || isempty(selected))
selected = true(size(or_vols));
end
final_vol_size = [size(or_vols{1}, 1), size(or_vols{1}, 2), size(or_vols{1}, 3)];
avg_R_vecs = zeros([final_vol_size, 3]);
avg_R_vecs_int = zeros(final_vol_size);
stddev_R_vecs = zeros(final_vol_size);
if (iscell(ors))
ors_sel = [ors{selected}];
% Converting to gvdm
gvdm_tmp = cat(1, ors_sel(:).R_vector)';
else
gvdm_tmp = ors(:, selected);
end
% computing weighted average orientation
for ii = 1:numel(or_vols)
weights = or_vols{ii};
r_vec = cat(4, ...
gvdm_tmp(1, ii) * ones(final_vol_size), ...
gvdm_tmp(2, ii) * ones(final_vol_size), ...
gvdm_tmp(3, ii) * ones(final_vol_size) );
avg_R_vecs_int = avg_R_vecs_int + weights;
avg_R_vecs = avg_R_vecs + weights(:, :, :, [1 1 1]) .* r_vec;
end
avg_R_vecs_int_norm = avg_R_vecs_int + (avg_R_vecs_int == 0); avg_R_vecs = avg_R_vecs ./ avg_R_vecs_int_norm(:, :, :, [1 1 1]);
% computing weighted standard deviations
for ii = 1:numel(or_vols)
r_vec = cat(4, ...
gvdm_tmp(1, ii) * ones(final_vol_size), ...
gvdm_tmp(2, ii) * ones(final_vol_size), ...
gvdm_tmp(3, ii) * ones(final_vol_size) );
square_dists = sum((r_vec - avg_R_vecs) .^ 2, 4);
weights = or_vols{ii};
stddev_R_vecs = stddev_R_vecs + weights .* square_dists;
end
stddev_R_vecs = stddev_R_vecs ./ avg_R_vecs_int_norm;
end
function avg_R_vec = computeAverageOrientation(gr, vols)
gr_sel = [gr{:}];
R_vecs = cat(1, gr_sel(:).R_vector);
vol_ints = GtOrientationSampling.computeODF(gr, vols);
vol_ints = reshape(vol_ints, [], 1);
vol_ints = vol_ints(:, [1 1 1]);
avg_R_vec = sum(R_vecs .* vol_ints, 1) ./ sum(vol_ints, 1);
end
function odf = computeODF(ors, vols)
odf = zeros(size(ors));
odf(:) = cellfun(@(x)gtMathsSumNDVol(x), vols);
end
function [dmean, dcenters, dranges] = get_deformation_stats(gvdm)
dmean = mean(gvdm, 2)';
dcenters = (max(gvdm, [], 2)' + min(gvdm, [], 2)') / 2;
dranges = max(gvdm, [], 2)' - min(gvdm, [], 2)';
end
end
end