From d2422b5a71b5ff9a54287c60432f5c267467882f Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Thu, 30 Apr 2015 14:18:46 +0200
Subject: [PATCH] 6D-Reconstruction: added initial twin reconstruction
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
.../gtReconstruct6DLaunchAlgorithm.m | 2 +-
.../gtReconstructGrainTwinCluster.m | 111 ++++
.../Gt6DReconstructionAlgorithmFactory.m | 116 +++-
.../gt6DCreateProjDataFromTwinnedGrain.m | 606 ++++++++++++++++++
zUtil_Deformation/test_generation/do6D.m | 4 +-
5 files changed, 820 insertions(+), 19 deletions(-)
create mode 100644 5_reconstruction/gtReconstructGrainTwinCluster.m
create mode 100644 zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m
diff --git a/5_reconstruction/gtReconstruct6DLaunchAlgorithm.m b/5_reconstruction/gtReconstruct6DLaunchAlgorithm.m
index 2ac1d85b..890031b5 100644
--- a/5_reconstruction/gtReconstruct6DLaunchAlgorithm.m
+++ b/5_reconstruction/gtReconstruct6DLaunchAlgorithm.m
@@ -5,7 +5,7 @@ function [algo, good_or] = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, par
rec_factory = Gt6DReconstructionAlgorithmFactory(parameters, ...
'volume_downscaling', rec_opts.volume_downscaling, ...
'rspace_oversize', rec_opts.rspace_oversize );
- [algo, good_or] = rec_factory.getSubBlobReconstructionAlgo(sampler, rec_opts.num_interp);
+ [algo, good_or] = rec_factory.getReconstructionAlgo(sampler, rec_opts.num_interp);
mem_cons = [ ...
algo.get_peak_memory_consumption(rec_opts.algorithm), ...
diff --git a/5_reconstruction/gtReconstructGrainTwinCluster.m b/5_reconstruction/gtReconstructGrainTwinCluster.m
new file mode 100644
index 00000000..e7bf3881
--- /dev/null
+++ b/5_reconstruction/gtReconstructGrainTwinCluster.m
@@ -0,0 +1,111 @@
+function gtReconstructGrainTwinCluster(grain_ids, phase_id, parameters, varargin)
+% gtReconstructGrainOrientation 6D reconstructions on a GPU machine
+% gtAstraReconstructGrain(grainIDs, phaseID, [parameters])
+% -------------------------------------------------------
+ if (~exist('parameters', 'var') || isempty(parameters))
+ parameters = gtLoadParameters();
+ end
+ parameters.fsim.check_spot = false;
+
+ rec_opts = gtReconstruct6DGetParamenters(parameters);
+
+ phase_dir = fullfile(parameters.acq.dir, '4_grains', ...
+ sprintf('phase_%02d', phase_id));
+ grains_str_ids = sprintf('_%04d', grain_ids);
+ grain_cluster_file = fullfile(phase_dir, ...
+ sprintf('grain_twin_cluster%s.mat', grains_str_ids));
+
+ fprintf('Loading the cluster file..')
+ grs = load(grain_cluster_file);
+ grs = grs.samp_ors;
+ fprintf('\b\b: Done.\n')
+ num_grains = numel(grs);
+
+ for ii_g = 1:num_grains
+ bb_gvdm = cat(1, grs(ii_g).bb_ors(:).R_vector)';
+
+ sampler(ii_g) = GtOrientationSampling(parameters, grs(ii_g)); %#ok<AGROW>
+ sampler(ii_g).make_even_simple_grid('cubic', rec_opts.grid_edge, bb_gvdm, 1)
+ if (rec_opts.super_sampling > 1)
+ sampler(ii_g).make_supersampling_simple_grid([1 2 3], rec_opts.super_sampling);
+ end
+ end
+
+ [algo, good_or] = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters);
+
+ vols = algo.getCurrentSolution();
+
+ % Restoring initial volume size (depending on the rounding)
+ if (rec_opts.volume_downscaling > 1)
+ fprintf('Expanding volumes..')
+ c = tic();
+ vols = gtMathsUpsampleVolume(vols, rec_opts.volume_downscaling);
+ fprintf('\b\b: Done (%f seconds).\n', toc(c))
+ end
+
+ ranges = zeros(num_grains, 2);
+ ranges(1, 1) = 1;
+ for ii_g = 2:num_grains
+ ranges(ii_g, 1) = ranges(ii_g-1) + sampler(ii_g-1).get_total_num_orientations();
+ end
+ ranges(1:end-1, 2) = ranges(2:end, 1) - 1;
+ ranges(end, 2) = numel(vols);
+
+ for ii_g = num_grains:-1:1
+ or_sizes = sampler(ii_g).get_orientation_sampling_size();
+ gr_vols = vols(ranges(ii_g, 1):ranges(ii_g, 2));
+ gr_good_or = good_or(ranges(ii_g, 1):ranges(ii_g, 2));
+
+ [avg_R_vecs, avg_R_vecs_int] = sampler(ii_g).getAverageOrientations(gr_vols, gr_good_or);
+ avg_R_vec = sampler(ii_g).getAverageOrientation(gr_vols, gr_good_or);
+ s_g_odf = reshape(sampler(ii_g).getODF(gr_vols, gr_good_or), or_sizes{1});
+
+ vol_size = size(avg_R_vecs_int);
+ shift = gtFwdSimComputeVolumeShifts(grs(ii_g).proj, parameters, vol_size);
+
+ [kam, gam(ii_g)] = gtDefComputeKernelAverageMisorientation(avg_R_vecs, avg_R_vecs_int);
+ [igm, gos(ii_g)] = gtDefComputeIntraGranularMisorientation(avg_R_vecs, avg_R_vecs_int, 'R_vector', grs(ii_g).R_vector);
+
+ ODF6D(ii_g) = struct( ...
+ 'voxels_avg_R_vectors', {avg_R_vecs}, ...
+ 'intensity', {avg_R_vecs_int}, ...
+ 'shift', {shift}, ...
+ 'R_vectors', {sampler(ii_g).get_R_vectors()}, ...
+ 'single_grain_ODF', {s_g_odf}, ...
+ 'single_grain_avg_R_vector', {avg_R_vec}, ...
+ 'kernel_average_misorientation', {kam}, ...
+ 'intra_granular_misorientation', {igm} );
+ end
+
+ cluster_rec_file = fullfile(phase_dir, ...
+ sprintf('grain_twin_cluster_details%s.mat', grains_str_ids));
+ % Saving and cleaning at the same time
+ if (exist(cluster_rec_file, 'file'))
+ cl_rec = load(cluster_rec_file);
+ cl_rec.ODF6D = ODF6D;
+ else
+ cl_rec = struct('ODF6D', ODF6D); %#ok<NASGU>
+ end
+ save(cluster_rec_file, '-struct', 'cl_rec', '-v7.3');
+
+ if (algo.verbose)
+ clear('ODF6D')
+ [proj_blobs, proj_spots] = algo.getProjectionOfCurrentSolution();
+
+ for ii_g = num_grains:-1:1
+ gr_vols = vols(ranges(ii_g, 1):ranges(ii_g, 2));
+
+ ODF6D(ii_g) = struct( ...
+ 'orientation_volumes', {gr_vols}, ...
+ 'fwd_projected_blobs', {proj_blobs}, ...
+ 'fwd_projected_spots', {proj_spots}, ...
+ 'grain_average_misorientation', {gam(ii_g)}, ...
+ 'grain_orientation_spread', {gos(ii_g)}, ...
+ 'compute_statistics', {algo.get_statistics()} );
+ end
+ grain_full_details_file = fullfile(phase_dir, ...
+ sprintf('grain_twin_cluster_full_details%s.mat', grains_str_ids));
+ save(grain_full_details_file, 'ODF6D', '-v7.3');
+ end
+end
+
diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
index 03743b0b..fc8a06e6 100644
--- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
+++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
@@ -21,7 +21,15 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
self.parameters = parameters;
end
- function [algo, good_orients, blobs] = getSubBlobReconstructionAlgo(self, sampler, num_interp, varargin)
+ function [algo, good_orients, blobs] = getReconstructionAlgo(self, sampler, num_interp, varargin)
+ if (numel(sampler) == 1)
+ [algo, good_orients, blobs] = self.getGrainReconstructionAlgo(sampler, num_interp, varargin{:});
+ else
+ [algo, good_orients, blobs] = self.getTwinReconstructionAlgo(sampler, num_interp, varargin{:});
+ end
+ end
+
+ function [algo, good_orients, blobs] = getGrainReconstructionAlgo(self, sampler, num_interp, varargin)
algo_params = self.get_algorithm_params(sampler, num_interp);
% Build Empty volumes
@@ -31,19 +39,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
if (self.volume_downscaling > 1)
volume_size = ceil(volume_size / self.volume_downscaling);
- num_geoms = numel(algo_params.geometries);
- for ii = 1:num_geoms
- algo_params.geometries{ii}(:, 4:12) = algo_params.geometries{ii}(:, 4:12) / self.volume_downscaling;
- end
- if (~isempty(algo_params.geometries_ss))
- for ii = 1:num_geoms
- geom_ss = algo_params.geometries_ss{ii};
- for ii_ss = 1:numel(geom_ss)
- geom_ss{ii_ss}(:, 4:12) = geom_ss{ii_ss}(:, 4:12) / self.volume_downscaling;
- end
- algo_params.geometries_ss{ii} = geom_ss;
- end
- end
end
good_orients = algo_params.good_orients;
@@ -60,6 +55,78 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
'psf', algo_params.psf, ...
varargin{:} );
end
+
+ function [algo, good_orients, blobs] = getTwinReconstructionAlgo(self, sampler, num_interp, varargin)
+ num_grains = numel(sampler);
+ for ii_g = 1:num_grains
+ algo_params(ii_g) = self.get_algorithm_params(sampler(ii_g), num_interp); %#ok<AGROW>
+ end
+
+ % Build Empty volumes
+ ref_gr = sampler(1).get_reference_grain();
+ spacing = mean([ref_gr.proj.vol_size_y, ref_gr.proj.vol_size_x, ref_gr.proj.vol_size_z]) * (self.rspace_oversize - 1);
+ volume_size = round([ref_gr.proj.vol_size_y, ref_gr.proj.vol_size_x, ref_gr.proj.vol_size_z] + spacing);
+
+ if (self.volume_downscaling > 1)
+ volume_size = ceil(volume_size / self.volume_downscaling);
+ end
+
+ blobs = cat(1, algo_params(:).blobs);
+ psf = cat(1, algo_params(:).psf);
+
+ good_orients = cat(1, algo_params(:).good_orients);
+
+ geometries = cat(1, algo_params(:).geometries);
+ geometries_ss = cat(1, algo_params(:).geometries_ss);
+
+ offsets = algo_params(1).offsets;
+ offsets_ss = algo_params(1).offsets_ss;
+ 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);
+
+ 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;
+ 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))
+ orient_offsets{ii_o}.proj(ii_p).blob_offset ...
+ = ii_b + base_offset_shift;
+ end
+ end
+ end
+ offsets = cat(1, offsets, orient_offsets);
+
+ if (~isempty(geometries_ss))
+ orient_offsets_ss = algo_params(ii_g).offsets_ss;
+ 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;
+ end
+ end
+ offsets_ss = cat(1, offsets_ss, orient_offsets_ss);
+ end
+
+ base_offset_shift = base_offset_shift + numel(algo_params(ii_g).blobs);
+ end
+
+ algo = Gt6DBlobReconstructor(volume_size, blobs, ...
+ 'data_type', self.data_type, ...
+ 'geometries', geometries, ...
+ 'offsets', offsets, ...
+ 'ss_geometries', geometries_ss, ...
+ 'ss_offsets', offsets_ss, ...
+ 'use_astra_projectors', true, ...
+ 'volume_ss', self.volume_downscaling, ...
+ 'psf', psf, ...
+ varargin{:} );
+ end
end
methods (Access = protected)
@@ -246,6 +313,21 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
offsets_ss = {};
end
+ if (self.volume_downscaling > 1)
+ for ii = 1:num_geoms
+ geometries{ii}(:, 4:12) = geometries{ii}(:, 4:12) / self.volume_downscaling;
+ end
+ if (~isempty(geometries_ss))
+ for ii = 1:num_geoms
+ geom_ss = geometries_ss{ii};
+ for ii_ss = 1:numel(geom_ss)
+ geom_ss{ii_ss}(:, 4:12) = geom_ss{ii_ss}(:, 4:12) / self.volume_downscaling;
+ end
+ geometries_ss{ii} = geom_ss;
+ end
+ end
+ end
+
if (isfield(bl, 'psf'))
psf = arrayfun(@(x){ permute(x.psf, [1 3 2]) }, bl);
elseif (isfield(sampler.ref_gr.proj, 'psf'))
@@ -261,7 +343,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
'geometries_ss', {geometries_ss}, ...
'offsets_ss', {offsets_ss}, ...
'psf', {psf}, ...
- 'good_orients', {good_orients} );
+ 'good_orients', {good_orients'} );
end
function [sub_blob_slices, proj_coeffs] = get_sub_blobs(self, blobs, slices_interp, padding)
@@ -441,6 +523,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
else
geometries = {grains_props(:).proj_geom};
end
+ geometries = reshape(geometries, [], 1);
+
num_geoms = numel(geometries);
w_int_tab(:, :, 1) = floor(w_tab);
@@ -451,7 +535,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
slice_pos_in_blob = w_int_tab - repmat(extreemes_blobs_w(:, 1), [1 num_geoms 2]) + 1;
% Assign offsets
- offsets = cell(1, num_geoms);
+ offsets = cell(num_geoms, 1);
for ii = 1:num_geoms
valid_proj_dirs = ...
(w_tab(:, ii) >= abs_lims(:, 1)) ...
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m
new file mode 100644
index 00000000..c7f51d01
--- /dev/null
+++ b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m
@@ -0,0 +1,606 @@
+function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFromTwinnedGrain(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_twin_spots', false, ...
+ 'mask_spots', true, ...
+ 'oversize', 1.4, ...
+ 'save', false );
+ conf = parse_pv_pairs(conf, varargin);
+
+ 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));
+ 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
+ p = gtLoadParameters();
+ symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
+
+ 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);
+ twingr = grs(2);
+
+ if (isfield(refgr.proj, 'ondet'))
+ ref_ondet = refgr.proj.ondet;
+ ref_included = refgr.proj.included;
+ ref_selected = refgr.proj.selected;
+
+ twin_ondet = twingr.proj.ondet;
+ twin_inc = twingr.proj.included;
+ twin_sel = twingr.proj.selected;
+ else
+ ref_ondet = refgr.ondet;
+ ref_included = refgr.included;
+ ref_selected = refgr.selected;
+
+ twin_ondet = twingr.ondet;
+ twin_inc = twingr.included;
+ twin_sel = twingr.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
+
+ % let's now filter the reflections that are not available for the twin
+ bool_ref_ondet = false(size(refgr.allblobs.omega));
+ bool_ref_ondet(ref_ondet) = true;
+ 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;
+
+ bool_twin_ondet = false(size(refgr.allblobs.omega));
+ bool_twin_ondet(twin_ondet) = true;
+ 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_deleted_ref_ondet = bool_ref_ondet & ~bool_twin_ondet;
+% bool_deleted_ref_inc = bool_ref_inc & ~bool_twin_ondet;
+%
+% ref_deleted_ondet = find(bool_deleted_ref_ondet);
+% ref_deleted_included = find(bool_deleted_ref_inc(ref_deleted_ondet));
+
+ bool_ref_ondet = bool_ref_ondet & bool_twin_ondet;
+ bool_ref_inc = bool_ref_inc & bool_twin_ondet;
+ bool_ref_sel = bool_ref_sel & bool_twin_ondet;
+
+% bool_twin_ondet = bool_ref_ondet & bool_twin_ondet;
+% bool_twin_inc = bool_ref_inc & bool_twin_inc;
+% bool_twin_sel = bool_ref_sel & bool_twin_sel;
+
+ ref_ondet = find(bool_ref_ondet);
+ ref_included = find(bool_ref_inc(ref_ondet));
+ ref_selected = bool_ref_sel(ref_ondet(ref_included));
+
+ if (isfield(refgr.proj, 'ondet'))
+ refgr.proj.ondet = ref_ondet;
+ refgr.proj.included = ref_included;
+ refgr.proj.selected = ref_selected;
+ else
+ refgr.ondet = ref_ondet;
+ refgr.included = ref_included;
+ refgr.selected = ref_selected;
+ end
+
+% ref_twin_ondet_bool = bool_twin_ondet(ref_ondet);
+% ref_twin_included_bool = bool_twin_inc(ref_ondet(ref_included));
+% ref_twin_selected_bool = bool_twin_sel(ref_ondet(ref_included(ref_selected)));
+
+ % And now..
+ fwdsim_inc_refl = ref_ondet(ref_included);
+
+ 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', ...
+ [(1:numel(vis_sel))', ...
+ grs(1).allblobs.omega(vis_sel), ...
+ grs(2).allblobs.omega(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, ...
+ bool_twin_inc(vis_sel), bool_twin_sel(vis_sel)]');
+ 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_bboxes = cell(num_grains, 1);
+ img_sizes = cell(num_grains, 1);
+
+ 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
+
+ % Let's now compute the bb on the images, by computing for each corner
+ % of the space bb, the position on the detector of each corner of the
+ % orientation bb.
+ or = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p);
+
+ refor = struct( ...
+ 'id', refgr.id, 'phaseid', refgr.phaseid, 'gids', grs_list, ...
+ 'center', estim_space_bbox_mm(1:3) + bbox_size_mm / 2, ...
+ 'R_vector', estim_orient_bbox(1:3) + bbox_size_rod / 2 );
+ refor = gtCalculateGrain(refor, p);
+ refor = reorder_reflections(refor, refgr);
+
+ uvw_tab = get_uvw_tab(or, fwdsim_inc_refl);
+
+ refor_ws = refor.allblobs.omega(fwdsim_inc_refl) / p.labgeo.omstep;
+ uvw_tab = recenter_w_tab(uvw_tab, refor_ws, p);
+
+ img_bboxes{ii_g} = [
+ squeeze(floor(min(uvw_tab, [], 2))), ...
+ squeeze( ceil(max(uvw_tab, [], 2))) ];
+ img_sizes{ii_g} = img_bboxes{ii_g}(:, 4:6) - img_bboxes{ii_g}(:, 1:3) + 1;
+
+ refor.bb_ors = or;
+
+ samp_ors(ii_g) = refor;
+ end
+
+ % We avoid the vertical spots for convenience
+ bad_etas = false(size(fwdsim_inc_refl));
+ for ii_g = 1:num_grains
+ refor_ns = samp_ors(ii_g).allblobs.eta(fwdsim_inc_refl);
+ bad_etas = bad_etas | acosd(abs(cosd(refor_ns))) < conf.min_eta;
+ if (conf.verbose)
+ fprintf('Grain: %d\n', ii_g)
+ fprintf('%02d) du %8d, dv %8d, dw %8d, eta: %7.3f <- used: %d, u: [%4d %4d], v: [%4d %4d], w: [%4d %4d]\n', ...
+ [(1:numel(bad_etas))', img_sizes{1}, refor_ns, ~bad_etas, img_bboxes{1}(:, [1 4 2 5 3 6]) ]');
+ fprintf('\n')
+ end
+ end
+
+ for ii_g = 1:num_grains
+ img_bboxes{ii_g} = img_bboxes{ii_g}(~bad_etas, :);
+ img_sizes{ii_g} = img_sizes{ii_g}(~bad_etas, :);
+ end
+
+ [stackUSize, stackVSize] = get_stack_UV_size(cat(1, img_sizes{:}), p, conf);
+
+ inc_reflections = find(1:numel(fwdsim_inc_refl));
+ inc_reflections = inc_reflections(~bad_etas);
+
+ ref_BBus = cat(1, refgr.proj.bl(inc_reflections).mbbu);
+ ref_BBvs = cat(1, refgr.proj.bl(inc_reflections).mbbv);
+ ref_BBsizes = cat(1, refgr.proj.bl(inc_reflections).mbbsize);
+ ref_norm_BBs = [ref_BBus(:, 1), ref_BBvs(:, 1), ref_BBsizes(:, 1:2)];
+
+ ref_BBus = (ref_BBus + cat(1, refgr.proj.bl(inc_reflections).bbuim)) / 2;
+ ref_BBvs = (ref_BBvs + cat(1, refgr.proj.bl(inc_reflections).bbvim)) / 2;
+ ref_BBsizes = (ref_BBsizes + cat(1, refgr.proj.bl(inc_reflections).bbsize)) / 2;
+ ref_cut_BBs = [ref_BBus(:, 1), ref_BBvs(:, 1), ref_BBsizes(:, 1:2)];
+
+ inc_ref_ondet = ref_ondet(ref_included(inc_reflections));
+
+ inc_ref_omegas = samp_ors(1).allblobs.omega(inc_ref_ondet);
+
+ for ii_g = 1:num_grains
+ fprintf('Loading raw images for grain %d: ', ii_g)
+
+ 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;
+ end
+
+ temp_gr = grs(ii_g);
+ temp_gr.center = refgr.center;
+ proj_bls = gtFwdSimPredictProjectedGrainBBox(temp_gr, ref_norm_BBs, inc_ref_omegas, inc_ref_ondet, p);
+ proj_cut_bls = gtFwdSimPredictProjectedGrainBBox(temp_gr, ref_cut_BBs, inc_ref_omegas, inc_ref_ondet, p);
+
+ num_blobs = numel(inc_reflections);
+ blobs(1:num_blobs) = get_default_blob();
+
+ for ii_b = 1:num_blobs
+ num_chars = fprintf('%02d/%02d', ii_b, num_blobs);
+
+ if (ii_g == 1 || ~shared_parent(ii_b))
+ gbl_norm = proj_bls(ii_b);
+ gbl_cut = proj_cut_bls(ii_b);
+ blobs(ii_b) = load_blob(img_bboxes{ii_g}(ii_b, :), ...
+ img_sizes{ii_g}(ii_b, :), stackUSize, stackVSize, gbl_norm, gbl_cut, p); %#ok<AGROW>
+ else
+ % Using the same because it gives the same size/uvw-limits
+ blobs(ii_b) = samp_ors(1).proj.bl(ii_b); %#ok<AGROW>
+ end
+
+ 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 = ref_ondet;
+ proj.included = ref_included(inc_reflections);
+ proj.selected = true(size(proj.included));
+
+ if (ii_g > 1)
+ proj.selected = samp_ors(1).proj.selected;
+ proj.shared_parent = shared_parent; % <- size: proj.included
+ else
+ proj.selected = true(size(proj.included));
+ end
+
+ samp_ors(ii_g).proj = proj;
+
+ if (conf.check_twin_spots)
+ samp_ors(ii_g).proj.selected = samp_ors(1).proj.selected & gtGuiGrainMontage(samp_ors(ii_g), [], true);
+ for ii_g_o = 1:ii_g-1
+ samp_ors(ii_g_o).proj.selected = samp_ors(ii_g_o).proj.selected & samp_ors(ii_g).proj.selected;
+ 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();
+ 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 = numel(orients):-1:1
+ 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 = 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
+
+
+
+
+
diff --git a/zUtil_Deformation/test_generation/do6D.m b/zUtil_Deformation/test_generation/do6D.m
index 8cc63fee..f279642d 100755
--- a/zUtil_Deformation/test_generation/do6D.m
+++ b/zUtil_Deformation/test_generation/do6D.m
@@ -186,7 +186,7 @@ switch(mode)
recon_algo_factory = Gt6DReconstructionAlgorithmFactory();
% [reconstructor, blobs] = recon_algo_factory.getBlobReconstructionAlgo( ...
% test_data.bl, blobs_on_det, sampler, '2D');
- [reconstructor, blobs] = recon_algo_factory.getSubBlobReconstructionAlgo( ...
+ [reconstructor, blobs] = recon_algo_factory.getReconstructionAlgo( ...
test_data.bl, blobs_on_det, sampler, '2D', num_interp);
case 'testPeterSubBlob3D'
cd('~/data/testdata/')
@@ -230,7 +230,7 @@ switch(mode)
end
recon_algo_factory = Gt6DReconstructionAlgorithmFactory();
- [reconstructor, blobs] = recon_algo_factory.getSubBlobReconstructionAlgo( ...
+ [reconstructor, blobs] = recon_algo_factory.getReconstructionAlgo( ...
test_data.bl, blobs_on_det, sampler, '3D', num_interp);
case 'peterDataset'
cd('/mntdirect/_data_id19_degisher/mi1026/id18/AlLi2b_vert_ ')
--
GitLab