From 2e1640722daa3cddfcd74b5eda67ef18aa32d493 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Wed, 13 May 2015 19:39:30 +0200
Subject: [PATCH] 6D-Reconstruction: added reconstruction support for parent +
twin
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
.../Gt6DReconstructionAlgorithmFactory.m | 78 ++-
.../gt6DCreateProjDataFromTwinnedGrain.m | 92 +--
...t6DCreateProjDataFromTwinnedGrainFwdProj.m | 535 ++++++++++++++++++
3 files changed, 654 insertions(+), 51 deletions(-)
create mode 100644 zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
index fc8a06e6..5153113f 100644
--- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
+++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
@@ -9,6 +9,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
volume_downscaling = 1;
rspace_oversize = 1.2;
+ use_predicted_scatter_ints = true;
+
parameters;
end
@@ -59,7 +61,10 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
function [algo, good_orients, blobs] = getTwinReconstructionAlgo(self, sampler, num_interp, varargin)
num_grains = numel(sampler);
for ii_g = 1:num_grains
+ gr = sampler(ii_g).get_reference_grain();
+ fprintf('%d) Grainid %d:\n', ii_g, gr.id)
algo_params(ii_g) = self.get_algorithm_params(sampler(ii_g), num_interp); %#ok<AGROW>
+ fprintf('\n')
end
% Build Empty volumes
@@ -72,6 +77,12 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
end
blobs = cat(1, algo_params(:).blobs);
+ % This vector will tell us where we will find in the blobs, the
+ % ones that were selected among the included, because the
+ % information we have in the twins about the spots that
+ % coincide are about the included and not the selected
+ included_pos_in_selected = zeros(size(ref_gr.proj.included));
+ included_pos_in_selected(ref_gr.proj.selected) = 1:numel(find(ref_gr.proj.selected));
psf = cat(1, algo_params(:).psf);
good_orients = cat(1, algo_params(:).good_orients);
@@ -84,18 +95,27 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
base_offset_shift = numel(algo_params(1).blobs);
for ii_g = 2:num_grains
temp_gr = sampler(ii_g).get_reference_grain();
- shared = temp_gr.proj.shared_parent(sampler(ii_g).selected);
+ shared = temp_gr.proj.shared_parent(sampler(ii_g).selected) > 0;
+ shared_pos = temp_gr.proj.shared_parent(sampler(ii_g).selected);
orient_offsets = algo_params(ii_g).offsets;
for ii_o = 1:numel(orient_offsets)
to_be_summed = ~shared(orient_offsets{ii_o}.blob_offsets);
+
orient_offsets{ii_o}.blob_offsets(to_be_summed) ...
= orient_offsets{ii_o}.blob_offsets(to_be_summed) + base_offset_shift;
+
+ orient_offsets{ii_o}.blob_offsets(~to_be_summed) ...
+ = included_pos_in_selected(shared_pos(orient_offsets{ii_o}.blob_offsets(~to_be_summed)));
+
for ii_p = 1:numel(orient_offsets{ii_o}.proj)
ii_b = orient_offsets{ii_o}.proj(ii_p).blob_offset;
- if (~temp_gr.proj.shared_parent(ii_b))
+ if (~shared(ii_b))
orient_offsets{ii_o}.proj(ii_p).blob_offset ...
= ii_b + base_offset_shift;
+ else
+ orient_offsets{ii_o}.proj(ii_p).blob_offset ...
+ = included_pos_in_selected(shared_pos(ii_b));
end
end
end
@@ -106,8 +126,23 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
for ii_o = 1:numel(orient_offsets_ss)
for ii_ss = 1:numel(orient_offsets_ss{ii_o})
to_be_summed = ~shared(orient_offsets{ii_o}{ii_ss}.blob_offsets);
+
orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(to_be_summed) ...
= orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(to_be_summed) + base_offset_shift;
+
+ orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(~to_be_summed) ...
+ = included_pos_in_selected(shared_pos(orient_offsets_ss{ii_o}{ii_ss}.blob_offsets(~to_be_summed)));
+
+ for ii_p = 1:numel(orient_offsets{ii_o}{ii_ss}.proj)
+ ii_b = orient_offsets{ii_o}{ii_ss}.proj(ii_p).blob_offset;
+ if (~shared(ii_b))
+ orient_offsets{ii_o}{ii_ss}.proj(ii_p).blob_offset ...
+ = ii_b + base_offset_shift;
+ else
+ orient_offsets{ii_o}{ii_ss}.proj(ii_p).blob_offset ...
+ = included_pos_in_selected(shared_pos(ii_b));
+ end
+ end
end
end
offsets_ss = cat(1, offsets_ss, orient_offsets_ss);
@@ -162,6 +197,9 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
= self.get_sub_blob_lims(w_tab, bl, ref_ws, slices_interp, padding);
[sub_blob_slices, proj_coeffs] = self.get_sub_blobs(bl, slices_interp, padding);
+ if (self.use_predicted_scatter_ints)
+ sub_blob_slices = self.renormalize_blobs(sampler, sub_blob_slices);
+ end
[overflow_b, overflow_o] = self.compute_overflow(w_tab, extreemes_blobs_w);
@@ -181,11 +219,16 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
fprintf(['WARNING! Some selected blobs were filtered out, ' ...
'due to no orientation projecting to the said blob:\n'])
disp(find(~good_blobs)')
-
- titles = {'blobs img lims:', 'projs img lims:'};
- fprintf(' %15s %15s\n', titles{:})
- fprintf(' %4d %4d %4d %4d\n', ...
- [extreemes_blobs_w(~good_blobs, :), round(extreemes_projs_w(~good_blobs, :))]');
+ fprintf('Position in the included vector:\n')
+ incl_pos = false(size(sel_bls));
+ incl_pos(sel_incl_indx(~good_blobs)) = true;
+ disp(find(incl_pos)')
+
+ titles = {'blobs img lims:', 'blobs img lims:', 'projs img lims:'};
+ fprintf(' %15s %15s %15s\n', titles{:})
+ fprintf(' %4d %4d %4d %4d %4d %4d\n', ...
+ [extreemes_blobs_w(~good_blobs, :), cat(1, bl(~good_blobs).bbwim), ...
+ round(extreemes_projs_w(~good_blobs, :))]');
end
bl = bl(good_blobs);
sub_blob_slices = sub_blob_slices(good_blobs);
@@ -371,7 +414,9 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
% Let's build the interpolated slices
blob = blobs(ii).intm;
blob(blob < 0) = 0;
-% blob = blob ./ sum(blob(:));
+ if (self.use_predicted_scatter_ints)
+ blob = blob .* blobs(ii).intensity;
+ end
sub_blob_slices{ii} = zeros(blob_sizes(ii, 1), blob_sizes(ii, 2), nums_slices(ii), self.data_type);
for n = nums_slices(ii):-1:1
@@ -600,6 +645,23 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
geometries{ii} = double(geometries{ii}(valid_proj_dirs, :));
end
end
+
+ function blobs = renormalize_blobs(self, sampler, blobs)
+ ref_gr = sampler.get_reference_grain();
+ if (isfield(self.parameters.cryst(ref_gr.phaseid), 'int_exp'))
+ cryst_scatter_ints = self.parameters.cryst(ref_gr.phaseid).int_exp;
+ else
+ cryst_scatter_ints = self.parameters.cryst(ref_gr.phaseid).int;
+ end
+
+ sel = sampler.ondet(sampler.included(sampler.selected));
+ thetatypes = ref_gr.allblobs.thetatype(sel);
+ scatter_ints = cryst_scatter_ints(thetatypes);
+ lorentz_fact = ref_gr.allblobs.lorentzfac(sel);
+ for ii_b = 1:numel(blobs)
+ blobs{ii_b} = blobs{ii_b} ./ scatter_ints(ii_b) ./ lorentz_fact(ii_b);
+ end
+ end
end
methods (Access = protected, Static)
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m
index c7f51d01..f9a30f7b 100644
--- a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m
+++ b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m
@@ -22,12 +22,16 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
'save', false );
conf = parse_pv_pairs(conf, varargin);
+ p = gtLoadParameters();
+ symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
+
num_grains = numel(grs_list);
if (~isstruct(grs_list))
fprintf('Loading grains: ')
for ii_g = num_grains:-1:1
num_chars = fprintf('%02d/%02d (%d)', num_grains-ii_g+1, num_grains, grs_list(ii_g));
grs(ii_g) = gtLoadGrain(phase_id, grs_list(ii_g));
+ grs(ii_g) = gtCalculateGrain(grs(ii_g), p);
fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars]));
@@ -37,8 +41,6 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
grs = grs_list;
grs_list = cat(2, grs(:).id);
end
- p = gtLoadParameters();
- symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
if (conf.verbose)
fprintf('R_vectors:\n')
@@ -177,16 +179,19 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
if (conf.verbose)
vis_sel = fwdsim_inc_refl;
% vis_sel = ref_ondet;
- fprintf('%02d) w1 %6.2f, w2 %6.2f, w_inds (%d<->%d), hkls [%2d %2d %2d]<->[%2d %2d %2d] <- diff: %7.2f (!! %d) [twin: i %d, s %d]\n', ...
+ fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_inds (%d<->%d), hkls [%2d %2d %2d]<->[%2d %2d %2d] <- diff: %7.2f (!! %d, L %g) [twin: i %d, s %d]\n', ...
[(1:numel(vis_sel))', ...
grs(1).allblobs.omega(vis_sel), ...
grs(2).allblobs.omega(vis_sel), ...
+ grs(1).allblobs.eta(vis_sel), ...
+ grs(2).allblobs.eta(vis_sel), ...
grs(1).allblobs.omind(vis_sel), ...
grs(2).allblobs.omind(vis_sel), ...
grs(1).allblobs.hklsp(vis_sel, [1 2 4]), ...
grs(2).allblobs.hklsp(vis_sel, [1 2 4]), ...
- grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel), ...
- abs(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel)) < 1, ...
+ abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180), ...
+ abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180) < grs(1).allblobs.lorentzfac(vis_sel), ...
+ grs(1).allblobs.lorentzfac(vis_sel), ...
bool_twin_inc(vis_sel), bool_twin_sel(vis_sel)]');
end
@@ -292,6 +297,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
inc_ref_ondet = ref_ondet(ref_included(inc_reflections));
inc_ref_omegas = samp_ors(1).allblobs.omega(inc_ref_ondet);
+ inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(inc_ref_ondet);
for ii_g = 1:num_grains
fprintf('Loading raw images for grain %d: ', ii_g)
@@ -299,7 +305,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
if (ii_g > 1)
w_diff = inc_ref_omegas - samp_ors(ii_g).allblobs.omega(inc_ref_ondet);
% Could be improved
- shared_parent = abs(w_diff) < 1;
+ shared_parent = abs(mod(w_diff + 180, 360) - 180) < inc_ref_Lfac;
end
temp_gr = grs(ii_g);
@@ -361,43 +367,43 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
end
end
- partial_ints = zeros(numel(inc_reflections), num_grains);
- phys_renorm = ones(numel(inc_reflections), num_grains);
-
- partial_ints(:, 1) = cat(1, samp_ors(1).proj.bl.intensity);
- for ii_g = 2:num_grains
- to_be_summed = ~samp_ors(ii_g).proj.shared_parent;
- partial_ints(to_be_summed, ii_g) = cat(1, samp_ors(ii_g).proj.bl(to_be_summed).intensity);
- end
-
- if (~conf.flat_normalization)
- for ii_g = 1:num_grains
- try
- temp_or = gtComputeIncomingBeamIntensity(samp_ors(ii_g), p);
- beam_ints = temp_or.beam_intensity;
- catch mexc
- gtPrintException(mexc, 'Turning beam inensity correction off, because interlaced acquisition is not supported')
- beam_ints = 1;
- end
- phys_renorm(:, ii_g) = 1 ...
- ./ samp_ors(ii_g).allblobs.lorentzfac(inc_ref_ondet) ...
- ./ (beam_ints / mean(beam_ints));
- end
- end
- total_ints = sum(partial_ints .* phys_renorm, 2);
-
- % Let's renormalize
- for ii_g = 1:num_grains
- for ii_b = 1:num_blobs
- ints_renorm = phys_renorm(:, ii_g) ./ total_ints(:);
-
- if (ii_g > 1 && samp_ors(ii_g).proj.shared_parent(ii_b))
- fprintf('Shared!\n')
- else
- samp_ors(ii_g).proj.bl(ii_b).intm = samp_ors(ii_g).proj.bl(ii_b).intm * ints_renorm(ii_b);
- end
- end
- end
+% partial_ints = zeros(numel(inc_reflections), num_grains);
+% phys_renorm = ones(numel(inc_reflections), num_grains);
+%
+% partial_ints(:, 1) = cat(1, samp_ors(1).proj.bl.intensity);
+% for ii_g = 2:num_grains
+% to_be_summed = ~samp_ors(ii_g).proj.shared_parent;
+% partial_ints(to_be_summed, ii_g) = cat(1, samp_ors(ii_g).proj.bl(to_be_summed).intensity);
+% end
+%
+% if (~conf.flat_normalization)
+% for ii_g = 1:num_grains
+% try
+% temp_or = gtComputeIncomingBeamIntensity(samp_ors(ii_g), p);
+% beam_ints = temp_or.beam_intensity;
+% catch mexc
+% gtPrintException(mexc, 'Turning beam inensity correction off, because interlaced acquisition is not supported')
+% beam_ints = 1;
+% end
+% phys_renorm(:, ii_g) = 1 ...
+% ./ samp_ors(ii_g).allblobs.lorentzfac(inc_ref_ondet) ...
+% ./ (beam_ints / mean(beam_ints));
+% end
+% end
+% total_ints = sum(partial_ints .* phys_renorm, 2);
+%
+% % Let's renormalize
+% for ii_g = 1:num_grains
+% for ii_b = 1:num_blobs
+% ints_renorm = phys_renorm(:, ii_g) ./ total_ints(:);
+%
+% if (ii_g > 1 && samp_ors(ii_g).proj.shared_parent(ii_b))
+% fprintf('Shared!\n')
+% else
+% samp_ors(ii_g).proj.bl(ii_b).intm = samp_ors(ii_g).proj.bl(ii_b).intm * ints_renorm(ii_b);
+% end
+% end
+% end
if (conf.verbose)
f = figure();
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
new file mode 100644
index 00000000..a9ce90a8
--- /dev/null
+++ b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
@@ -0,0 +1,535 @@
+function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFromTwinnedGrainFwdProj(grs_list, phase_id, varargin)
+% FUNCTION [refor] = gt6DCreateProjDataFromGrainCluster(grs_list, phase_id, varargin)
+% proj: is a grain.proj structure
+% refor: is a grain structure for the average orientation in the average
+% point
+% or: is a set of grain structures for the extreeme corners of the
+% orientation space that should be considered
+
+ if (~exist('phase_id', 'var'))
+ phase_id = 1;
+ end
+
+ conf = struct( ...
+ 'verbose', false, ...
+ 'min_eta', 15, ...
+ 'ospace_oversize', 1.1, ...
+ 'rspace_oversize', 1.1, ...
+ 'flat_normalization', true, ...
+ 'check_spots', false, ...
+ 'mask_spots', true, ...
+ 'oversize', 1.4, ...
+ 'save', false );
+ conf = parse_pv_pairs(conf, varargin);
+
+ p = gtLoadParameters();
+ symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
+
+ num_grains = numel(grs_list);
+ if (~isstruct(grs_list))
+ fprintf('Loading grains: ')
+ for ii_g = num_grains:-1:1
+ num_chars = fprintf('%02d/%02d (%d)', num_grains-ii_g+1, num_grains, grs_list(ii_g));
+ grs(ii_g) = gtLoadGrain(phase_id, grs_list(ii_g));
+ grs(ii_g) = gtCalculateGrain(grs(ii_g), p);
+ fprintf(repmat('\b', [1 num_chars]));
+ fprintf(repmat(' ', [1 num_chars]));
+ fprintf(repmat('\b', [1 num_chars]));
+ end
+ fprintf('Done.\n')
+ else
+ grs = grs_list;
+ grs_list = cat(2, grs(:).id);
+ end
+
+ if (conf.verbose)
+ fprintf('R_vectors:\n')
+ cat(1, grs(:).R_vector)
+ fprintf('Reciprocal disorientations:\n')
+ for ii_g_1 = 1:num_grains
+ grs_to_check = (ii_g_1+1):num_grains;
+ if (numel(grs_to_check) > 0)
+ fprintf(' - Disorientations from grain: %d\n', grs(ii_g_1).id)
+ for ii_g_2 = grs_to_check
+ dis_angle = gtDisorientation(grs(ii_g_1).R_vector', grs(ii_g_2).R_vector', symm);
+ fprintf(' + Grain %d: %f\n', grs(ii_g_2).id, dis_angle);
+ end
+ end
+ end
+ end
+
+ refgr = grs(1);
+ grs(2:end) = reorder_reflections(grs(2:end), refgr);
+
+ if (isfield(refgr.proj, 'ondet'))
+ ref_ondet = refgr.proj.ondet;
+ ref_included = refgr.proj.included;
+ ref_selected = refgr.proj.selected;
+ else
+ ref_ondet = refgr.ondet;
+ ref_included = refgr.included;
+ ref_selected = refgr.selected;
+ end
+
+ if (conf.verbose)
+ try
+ refgr = gtComputeIncomingBeamIntensity(refgr, p);
+ beam_ints = refgr.beam_intensity(ref_selected);
+ beam_ints = beam_ints / mean(beam_ints);
+ catch mexc
+ gtPrintException(mexc, 'Skipping beam intensity renormalization')
+ beam_ints = 1;
+ end
+ blob_ints = refgr.intensity(ref_selected);
+
+ vis_spots = ref_ondet(ref_included(ref_selected));
+ thetatypes = refgr.allblobs.thetatype(vis_spots);
+
+ if (isfield(p.cryst(phase_id), 'int_exp'))
+ predicted_ints = p.cryst(phase_id).int_exp(thetatypes)';
+
+ f = figure();
+ ax = axes('parent', f);
+ plot(ax, blob_ints)
+ hold(ax, 'on')
+ plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
+ predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
+ plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
+ predicted_ints = predicted_ints .* beam_ints;
+ plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
+ end
+
+ predicted_ints = p.cryst(phase_id).int(thetatypes)';
+
+ f = figure();
+ ax = axes('parent', f);
+ plot(ax, blob_ints)
+ hold(ax, 'on')
+ plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
+ predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
+ plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
+ predicted_ints = predicted_ints .* beam_ints;
+ plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
+
+ fprintf('%d) %d, %g\n', ...
+ [(1:numel(blob_ints))', refgr.allblobs.thetatype(vis_spots), ...
+ blob_ints - predicted_ints / mean(predicted_ints) * mean(blob_ints)]')
+ hold(ax, 'off')
+ end
+
+ if (conf.verbose)
+ vis_sel = ref_ondet(ref_included);
+% vis_sel = ref_ondet;
+ fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_inds (%d<->%d), hkls [%2d %2d %2d]<->[%2d %2d %2d] <- diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
+ [(1:numel(vis_sel))', ...
+ grs(1).allblobs.omega(vis_sel), ...
+ grs(2).allblobs.omega(vis_sel), ...
+ grs(1).allblobs.eta(vis_sel), ...
+ grs(2).allblobs.eta(vis_sel), ...
+ grs(1).allblobs.omind(vis_sel), ...
+ grs(2).allblobs.omind(vis_sel), ...
+ grs(1).allblobs.hklsp(vis_sel, [1 2 4]), ...
+ grs(2).allblobs.hklsp(vis_sel, [1 2 4]), ...
+ abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180), ...
+ abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180) < grs(1).allblobs.lorentzfac(vis_sel), ...
+ grs(1).allblobs.lorentzfac(vis_sel), ...
+ acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
+ acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
+ ]');
+ end
+
+ vol_half_size = [refgr.proj.vol_size_y, refgr.proj.vol_size_x, refgr.proj.vol_size_z] / 2;
+ space_bbox = [ ...
+ refgr.proj.centerpix - vol_half_size * conf.rspace_oversize, ...
+ refgr.proj.centerpix + vol_half_size * conf.rspace_oversize ];
+ estim_space_bbox_pix = [floor(space_bbox(1:3)), ceil(space_bbox(4:6))];
+ bbox_size_pix = estim_space_bbox_pix(4:6) - estim_space_bbox_pix(1:3);
+
+ estim_space_bbox_mm = [ ...
+ gtGeoSam2Sam(estim_space_bbox_pix(1:3), p.recgeo, p.samgeo, false, false), ...
+ gtGeoSam2Sam(estim_space_bbox_pix(4:6), p.recgeo, p.samgeo, false, false) ];
+ bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3);
+
+ img_sizes = zeros(num_grains, 2);
+
+ for ii_g = num_grains:-1:1
+ sampler = GtOrientationSampling(p, grs(ii_g), 'verbose', conf.verbose);
+ r_vecs = sampler.guess_ODF_BB()';
+
+ estim_orient_bbox = [min(r_vecs, [], 1), max(r_vecs, [], 1)];
+ bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
+
+ % oversizing the orienation a bit
+ delta_bbox_size_rod = bbox_size_rod * 0.05;
+ estim_orient_bbox = estim_orient_bbox + [-delta_bbox_size_rod, delta_bbox_size_rod];
+ bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
+
+ bbox_size_deg = 2 * atand(bbox_size_rod);
+
+ if (conf.verbose)
+ fprintf('\n');
+ fprintf('Estimated spatial voxel BBox: [%3d, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix);
+ fprintf(' BBox size: %3d, %3d, %3d (%f, %f, %f mm)\n', bbox_size_pix, bbox_size_mm);
+
+ fprintf(' Estimated orientation BBox: [%3.3f, %3.3f, %3.3f] -> [%3.3f, %3.3f, %3.3f]\n', estim_orient_bbox);
+ fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
+ fprintf('\n');
+ end
+
+ refor = struct( ...
+ 'id', grs(ii_g).id, 'phaseid', grs(ii_g).phaseid, ...
+ 'center', estim_space_bbox_mm(1:3) + bbox_size_mm / 2, ...
+ 'R_vector', grs(ii_g).R_vector );
+ refor = gtCalculateGrain(refor, p);
+
+ refor.bb_ors = get_6D_space_corners(estim_space_bbox_mm, ...
+ bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p);
+
+ img_sizes(ii_g, :) = [size(grs(ii_g).proj.stack, 1), size(grs(ii_g).proj.stack, 3), ];
+
+ samp_ors(ii_g) = refor;
+ end
+ samp_ors = reorder_reflections(samp_ors, refgr);
+
+ [stackUSize, stackVSize] = get_stack_UV_size(cat(1, img_sizes), p, conf);
+
+ bool_ref_inc = false(size(refgr.allblobs.omega));
+ bool_ref_inc(ref_ondet(ref_included)) = true;
+ bool_ref_sel = false(size(refgr.allblobs.omega));
+ bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;
+
+ for ii_g = 1:num_grains
+ fprintf('Loading raw images for grain %d: ', ii_g)
+
+ if (ii_g > 1)
+ twingr = grs(ii_g);
+
+ if (isfield(twingr.proj, 'ondet'))
+ twin_ondet = twingr.proj.ondet;
+ twin_inc = twingr.proj.included;
+ twin_sel = twingr.proj.selected;
+ else
+ twin_ondet = twingr.ondet;
+ twin_inc = twingr.included;
+ twin_sel = twingr.selected;
+ end
+
+ bool_twin_inc = false(size(refgr.allblobs.omega));
+ bool_twin_inc(twin_ondet(twin_inc)) = true;
+ bool_twin_sel = false(size(refgr.allblobs.omega));
+ bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;
+
+ bool_test_inc = bool_ref_inc | bool_twin_inc;
+
+ inc_ref_omegas = samp_ors(1).allblobs.omega(bool_test_inc);
+ inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(bool_test_inc);
+ inc_ref_etas = samp_ors(1).allblobs.eta(bool_test_inc);
+
+ inc_twin_omegas = samp_ors(ii_g).allblobs.omega(bool_test_inc);
+ inc_twin_etas = samp_ors(ii_g).allblobs.eta(bool_test_inc);
+
+ w_diff = abs(mod(inc_ref_omegas - inc_twin_omegas + 180, 360) - 180);
+ n_diff = abs(mod(inc_ref_etas - inc_twin_etas + 180, 360) - 180);
+
+ test_shared = (w_diff < inc_ref_Lfac) & (n_diff < 5);
+ % let's add them as included
+ bool_test_inc_indx = find(bool_test_inc);
+
+ shared_parent = false(size(refgr.allblobs.omega));
+ shared_parent(bool_test_inc_indx(test_shared)) = true;
+
+ % So if it is shared, we use the parent's blob
+ bool_twin_inc = bool_twin_inc | shared_parent;
+ % If the parent enables it, we enable it as well
+ bool_twin_sel = bool_twin_sel | (shared_parent & bool_ref_sel);
+
+ twin_inc = find(bool_twin_inc(twin_ondet));
+ twin_sel = bool_twin_sel(twin_ondet(twin_inc));
+
+ shared_parent_pos = bool_ref_inc & shared_parent;
+ shared_parent_pos = shared_parent_pos(ref_ondet(ref_included));
+
+ % We cut it down to the size of the included
+ shared_parent = shared_parent(twin_ondet(twin_inc));
+
+ local_ondet = twin_ondet;
+ local_included = twin_inc;
+ local_selected = twin_sel;
+ else
+ local_ondet = ref_ondet;
+ local_included = ref_included;
+ local_selected = ref_selected;
+ end
+
+ num_blobs = numel(local_included);
+
+ if (ii_g == 1)
+ blobs = grs(ii_g).proj.bl;
+ else
+ blobs = get_default_blob();
+ blobs(~shared_parent) = grs(ii_g).proj.bl;
+ % Using the same because it gives the same size/uvw-limits
+ blobs(shared_parent) = samp_ors(1).proj.bl(shared_parent_pos);
+ end
+
+ for ii_b = 1:num_blobs
+ num_chars = fprintf('%02d/%02d', ii_b, num_blobs);
+
+ blobs(ii_b) = include_blob(blobs(ii_b), stackUSize, stackVSize);
+
+ fprintf(repmat('\b', [1 num_chars]));
+ fprintf(repmat(' ', [1 num_chars]));
+ fprintf(repmat('\b', [1 num_chars]));
+ end
+ fprintf('Done.\n')
+
+ spots = arrayfun(@(x){sum(x.intm, 3)}, blobs);
+ spots = permute(cat(3, spots{:}), [1 3 2]);
+
+ proj.centerpix = estim_space_bbox_pix(1:3) + bbox_size_pix / 2;
+ proj.bl = blobs;
+ proj.stack = spots;
+ vol_size = bbox_size_pix + mean(bbox_size_pix) * 0.3;
+ proj.vol_size_x = vol_size(2);
+ proj.vol_size_y = vol_size(1);
+ proj.vol_size_z = vol_size(3);
+
+ proj.ondet = local_ondet;
+ proj.included = local_included;
+ proj.selected = local_selected;
+
+ if (ii_g > 1)
+ proj.shared_parent = zeros(size(shared_parent));
+ proj.shared_parent(shared_parent) = find(shared_parent_pos);
+ end
+
+ samp_ors(ii_g).proj = proj;
+
+ if (conf.check_spots)
+ samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true);
+ end
+ end
+
+ if (conf.verbose)
+ f = figure();
+ ax = axes('parent', f);
+ hold(ax, 'on');
+ for ii_g = 1:num_grains
+ gt6DPlotOrientationBBox(ax, cat(1, samp_ors(ii_g).bb_ors(:).R_vector));
+ scatter3(ax, grs(ii_g).R_vector(1), grs(ii_g).R_vector(2), grs(ii_g).R_vector(3), 30);
+ end
+ hold(ax, 'off');
+
+ drawnow();
+ end
+
+ if (conf.save)
+ str_ids = sprintf('_%04d', grs_list);
+ grain_filename = fullfile(p.acq.dir, '4_grains', ...
+ sprintf('phase_%02d', phase_id), ...
+ sprintf('grain_twin_cluster%s.mat', str_ids));
+ fprintf('Saving the cluster file "%s"..', grain_filename)
+ save(grain_filename, 'samp_ors', '-v7.3');
+ fprintf('\b\b: Done.\n')
+
+% fprintf('Saving to sample.mat..')
+% sample = GtSample.loadFromFile();
+% sample.phases{phase_id}.clusters(end+1) = ...
+% GtPhase.makeCluster(grs_list, refor.R_vector, estim_space_bbox_pix, estim_orient_bbox);
+% sample.saveToFile();
+% fprintf('\b\b: Done.\n')
+ end
+end
+
+function orients = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p)
+ orients = {};
+ for ii_x1 = 0:1
+ base_x1 = estim_space_bbox_mm(1:3) + [ii_x1 * bbox_size_mm(1), 0, 0];
+ for ii_x2 = 0:1
+ base_x2 = base_x1 + [0, ii_x2 * bbox_size_mm(2), 0];
+ for ii_x3 = 0:1
+ base_x3 = base_x2 + [0, 0, ii_x3 * bbox_size_mm(3)];
+ for ii_r1 = 0:1
+ base_r1 = estim_orient_bbox(1:3) + [ii_r1 * bbox_size_rod(1), 0, 0];
+ for ii_r2 = 0:1
+ base_r2 = base_r1 + [0, ii_r2 * bbox_size_rod(2), 0];
+ for ii_r3 = 0:1
+ base_r3 = base_r2 + [0, 0, ii_r3 * bbox_size_rod(3)];
+
+% fprintf(' center: %f, %f, %f, r_vec: %f, %f, %f\n', base_x3, base_r3)
+
+ orients{end+1} = struct( ...
+ 'id', refgr.id, 'phaseid', refgr.phaseid, ...
+ 'center', base_x3, 'R_vector', base_r3); %#ok<AGROW>
+ end
+ end
+ end
+ end
+ end
+ end
+ orients = gtCalculateGrain_p(orients, p);
+ orients = [orients{:}];
+
+ orients = reorder_reflections(orients, refgr);
+end
+
+function orients = reorder_reflections(orients, refgr)
+ refgr_omind = refgr.allblobs.omind;
+ refgr_omind_ind = [ ...
+ find(refgr_omind == 1); find(refgr_omind == 2); ...
+ find(refgr_omind == 3); find(refgr_omind == 4) ];
+
+ for ii_g = 1:numel(orients)
+ orients(ii_g) = gtGrainAllblobsFilterOrder(orients(ii_g), refgr_omind_ind);
+ end
+end
+
+function uvw_tab = recenter_w_tab(uvw_tab, refor_ws, p)
+ num_images = gtGetTotNumberOfImages(p);
+
+ % Let's treat those blobs at the w edge 360->0
+ % (from the sampled orientations perspective)
+ opposite_ws = mod(refor_ws + num_images / 2, num_images);
+ gt_ref_int = opposite_ws > num_images / 2;
+ lt_ref_int = ~gt_ref_int;
+ opposite_ws_repmat = opposite_ws(:, ones(1, size(uvw_tab(:, :, 3), 2)));
+
+ uvw_tab(gt_ref_int, :, 3) = uvw_tab(gt_ref_int, :, 3) ...
+ - num_images .* (uvw_tab(gt_ref_int, :, 3) > opposite_ws_repmat(gt_ref_int, :));
+ uvw_tab(lt_ref_int, :, 3) = uvw_tab(lt_ref_int, :, 3) ...
+ + num_images .* (uvw_tab(lt_ref_int, :, 3) < opposite_ws_repmat(lt_ref_int, :));
+end
+
+function uvw_tab = get_uvw_tab(orientations, refl_sel)
+ num_ors = numel(orientations);
+ or_abs = cat(1, orientations(:).allblobs);
+ uvw_tab = zeros(numel(refl_sel), num_ors, 3);
+ for ii_o = 1:num_ors
+ if (isfield(or_abs(1), 'detector'))
+ uvw_tab(:, ii_o, :) = or_abs(ii_o).detector(1).uvw(refl_sel, :);
+ else
+ uvw_tab(:, ii_o, :) = or_abs(ii_o).uvw(refl_sel, :);
+ end
+ end
+end
+
+function [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf)
+
+ max_img_sizes = [max(img_sizes(:, 1)), max(img_sizes(:, 2))];
+ stack_imgs_oversize = min(p.fsim.oversize, conf.oversize);
+ u_size = round(max_img_sizes(1) * stack_imgs_oversize);
+ v_size = round(max_img_sizes(2) * stack_imgs_oversize);
+
+ if (conf.verbose)
+ fprintf('\n');
+ fprintf(' Maximum images size: [%3d, %3d]\n', max_img_sizes);
+ fprintf('Stack images size (oversize: %1.2f): [%3d, %3d]\n', stack_imgs_oversize, u_size, v_size);
+ fprintf('\n');
+ end
+end
+
+function blob = load_blob(img_bboxes, img_sizes, stackUSize, stackVSize, gbl_norm, gbl_cut, p)
+ blob = get_default_blob();
+
+ bb = [img_bboxes(1, 1:2), img_sizes(1, 1:2)];
+ blob_vol = gtGetRawRoi(img_bboxes(3), img_bboxes(6), p.acq, bb);
+ blob_vol(blob_vol < 0) = 0;
+ blob_vol(isnan(blob_vol)) = 0;
+ blob_bb = [img_bboxes(1, 1:3), img_sizes(1, :)];
+
+ % Transposing to keep the same convention as spots
+ blob_vol = permute(blob_vol, [2 1 3]);
+
+ blob.mbbsize = blob_bb(4:6);
+ blob.mbbu = [blob_bb(1), blob_bb(1) + blob_bb(4) + 1];
+ blob.mbbv = [blob_bb(2), blob_bb(2) + blob_bb(5) + 1];
+ blob.mbbw = [blob_bb(3), blob_bb(3) + blob_bb(6) + 1];
+
+ blob_size_im = [stackUSize, stackVSize, blob_bb(6)+2];
+
+ shifts_blob = gtFwdSimGetStackShifts(stackUSize, stackVSize, blob_bb, false);
+ shifts_blob = [shifts_blob.u, shifts_blob.v, 1];
+
+ blob.intm = gtPlaceSubVolume( ...
+ zeros(blob_size_im, 'single'), single(blob_vol), shifts_blob);
+
+ blob.bbsize = blob_size_im;
+
+ blob_bb_im = [blob_bb(1, 1:3) - shifts_blob, blob_size_im];
+
+ blob.bbuim = [blob_bb_im(1), blob_bb_im(1) + blob_bb_im(4) - 1];
+ blob.bbvim = [blob_bb_im(2), blob_bb_im(2) + blob_bb_im(5) - 1];
+ blob.bbwim = [blob_bb_im(3), blob_bb_im(3) + blob_bb_im(6) - 1];
+
+ % We should build the mask from the segmented reflections of the
+ % grains
+ blob.mask = false(blob_size_im);
+
+ % spreading the flat mask over the full expected BBox in UVW
+ grain_mask = gbl_cut.mask(:, :, ones(size(blob.mask, 3)));
+
+ gbl_shifts_im = [gbl_cut.bbuim(1), gbl_cut.bbvim(1), blob.bbwim(1)]; % <- note w component: blob.bbwim(1)
+ blob_shifts_im = [blob.bbuim(1), blob.bbvim(1), blob.bbwim(1)];
+ mask_shifts = gbl_shifts_im - blob_shifts_im;
+
+ blob.mask = gtPlaceSubVolume(blob.mask, grain_mask, mask_shifts);
+ blob.intm(~blob.mask) = 0;
+
+ % We should build the mask from the segmented reflections of the
+ % grains
+ % for the moment we only use one (the reference one), to
+ % renormalize the different reflections
+ blob.mask = false(blob_size_im);
+
+ % spreading the flat mask over the full expected BBox in UVW
+ grain_mask = gbl_norm.mask(:, :, ones(size(blob.mask, 3)));
+
+ gbl_shifts_im = [gbl_norm.bbuim(1), gbl_norm.bbvim(1), blob.bbwim(1)]; % <- note w component: blob.bbwim(1)
+ blob_shifts_im = [blob.bbuim(1), blob.bbvim(1), blob.bbwim(1)];
+ mask_shifts = gbl_shifts_im - blob_shifts_im;
+
+ blob.mask = gtPlaceSubVolume(blob.mask, grain_mask, mask_shifts);
+
+ blob_int = sum(blob.intm(blob.mask));
+% blob.intm = blob.intm / blob_int;
+
+ blob.intensity = blob_int;
+end
+
+function blob = include_blob(blob, stackUSize, stackVSize)
+ new_bbsize = [stackUSize, stackVSize, blob.bbsize(3)];
+ shift = floor((new_bbsize - blob.bbsize) / 2);
+
+ new_bbuim = blob.bbuim - shift(1);
+ new_bbuim(2) = new_bbuim(1) + new_bbsize(1) - 1;
+
+ new_bbvim = blob.bbvim - shift(2);
+ new_bbvim(2) = new_bbvim(1) + new_bbsize(2) - 1;
+
+ blob.bbuim = new_bbuim;
+ blob.bbvim = new_bbvim;
+ blob.bbsize = new_bbsize;
+
+ blob.intm = gtPlaceSubVolume(zeros(new_bbsize, 'like', blob.intm), blob.intm, shift);
+ blob.mask = logical(gtPlaceSubVolume(zeros(new_bbsize, 'uint8'), blob.mask, shift));
+end
+
+function blob = get_default_blob()
+ blob = struct( ...
+ 'intm', [], ... % The normalized blob
+ 'mask', [], ... % The segmentation mask
+ 'bbuim', [], ... % Image BB on U
+ 'bbvim', [], ... % Image BB on V
+ 'bbwim', [], ... % Image BB on W <- Not strictly w
+ 'bbsize', [], ... % Image BoundingBox (includes margins)
+ 'mbbsize', [], ... % Real (Measured) Blob Bounding Box
+ 'mbbu', [], ... % Real (Measured) Blob BB on U
+ 'mbbv', [], ... % Real (Measured) Blob BB on V
+ 'mbbw', [], ... % Real (Measured) Blob BB on W <- Not strictly w
+ 'intensity', [] ); % Blob original intensity
+end
+
+
+
+
+
--
GitLab