From d2fb04cfe71653d54e1562d9c02e502ac14b8f33 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <vigano@yoda.esrf.fr>
Date: Thu, 19 May 2016 18:47:38 +0200
Subject: [PATCH] 6D-Reconstruction/multi-region: hopefuly fixed new num_interp
 behavior

Signed-off-by: Nicola Vigano <vigano@yoda.esrf.fr>
---
 .../Gt6DReconstructionAlgorithmFactory.m      | 206 +++++++++++++++---
 ...t6DCreateProjDataFromTwinnedGrainFwdProj.m |  64 +++---
 2 files changed, 206 insertions(+), 64 deletions(-)

diff --git a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
index 834ff7cb..f4d05fce 100644
--- a/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
+++ b/zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m
@@ -64,13 +64,13 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
         end
 
         function [algo, blobs] = get6DAlgorithmMultiRegion(self, sampler, num_interp, varargin)
-            num_grains = numel(sampler);
+            num_regions = numel(sampler);
 
             % Build Empty volumes
             ref_gr = sampler(1).get_reference_grain();
             volume_size = self.get_volume_size(ref_gr);
 
-            for ii_o_reg = 1:num_grains
+            for ii_o_reg = 1:num_regions
                 gr = sampler(ii_o_reg).get_reference_grain();
                 fprintf('%d) Grainid %d:\n', ii_o_reg, gr.id)
                 shape_funcs = self.get_shape_functions(sampler);
@@ -86,52 +86,192 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             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;
+
             % 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(proj.included));
-            included_pos_in_selected(proj.selected) = 1:numel(find(proj.selected));
+            included_pos_in_selected = zeros(size(ref_gr.proj(self.det_index).included));
+            included_pos_in_selected(proj.selected) = 1:numel(find(ref_gr.proj(self.det_index).selected));
+
+            num_parent_blobs = numel(algo_params(1).blobs);
+
+            % We should first go through all the shared blobs in each
+            % orientation and extend them if needed, and then extend the
+            % positions pointed by the orientations in the shared blobs
+            all_shared_blobs = find(ref_gr.proj(self.det_index).shared_bls_twins);
+
+            % First entry (in the third dimension) contains the final size
+            shared_blobs_w_lims = zeros(num_parent_blobs, 2, num_regions+1);
+            shared_blobs_w_lims(all_shared_blobs, :, 1) = algo_params(1).blobs_w_lims(all_shared_blobs, :);
+            shared_blobs_w_lims(all_shared_blobs, :, 2) = algo_params(1).blobs_w_lims(all_shared_blobs, :);
+
+            % The reference num_interps
+            shared_blobs_w_interp = zeros(num_parent_blobs, 1);
+            shared_blobs_w_interp(all_shared_blobs) = algo_params(1).blobs_w_interp(all_shared_blobs);
+
+            % Same as above, but rescaled with num_interp
+            rescaled_shared_blobs_w_lims(all_shared_blobs, :, 2) = zeros(num_parent_blobs, 2, num_regions+1);
+            rescaled_shared_blobs_w_lims(all_shared_blobs, :, 2) = shared_blobs_w_lims(all_shared_blobs, :, 2) ./ shared_blobs_w_interp(all_shared_blobs, [1 1]);
+
+            for ii_o_reg = 2:num_regions
+                curr_or = sampler(ii_o_reg).get_reference_grain();
+                curr_sel = sampler(ii_o_reg).selected;
+                % This contains the position of the shared blobs in
+                % parent's numbering
+                shared_pos = curr_or.proj(self.det_index).shared_bls_parent(curr_sel);
+                % This indicates the ones that are actually shared among
+                % the twin's blobs
+                shared = shared_pos > 0;
+
+                if (any(shared_blobs_w_interp(shared_pos(shared)) ~= algo_params(ii_o_reg).blobs_w_interp(shared)))
+                    warning('Gt6DReconstructionAlgorithmFactory:wrong_interpolation', ...
+                        ['For some reason, the num_interp used for the ' ...
+                        'main orientation-space region and the region ' ...
+                        'n.%d, are different.\nThe ' ...
+                        'reconstruction will not be correct!!'], ii_o_reg)
+                end
 
-            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) > 0;
-                shared_pos = temp_gr.proj.shared_parent(sampler(ii_g).selected);
+                shared_blobs_w_lims(all_shared_blobs, :, ii_o_reg+1) = algo_params(ii_o_reg).blobs_w_lims(shared, :);
+                rescaled_shared_blobs_w_lims(all_shared_blobs, :, ii_o_reg+1) = algo_params(ii_o_reg).blobs_w_lims(shared, :) ./ shared_blobs_w_interp(shared_pos(shared), [1 1]);
 
-                % Shifting projection to the concatenated blobs, unless that
-                % blob was shared between parent and twin
-                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);
+                shared_blobs_w_lims(all_shared_blobs, :, 1) = [ ...
+                    min(shared_blobs_w_lims(all_shared_blobs, 1), algo_params(ii_o_reg).blobs_w_lims(shared, 1)), ...
+                    max(shared_blobs_w_lims(all_shared_blobs, 2), algo_params(ii_o_reg).blobs_w_lims(shared, 2)) ];
+            end
 
-                            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;
+            rescaled_shared_blobs_w_lims(all_shared_blobs, :, 1) = shared_blobs_w_lims(all_shared_blobs, :, 1) ./ shared_blobs_w_interp(all_shared_blobs, [1 1]);
+            to_be_changed = any(rescaled_shared_blobs_w_lims(all_shared_blobs, :, 1) ~= rescaled_shared_blobs_w_lims(all_shared_blobs, :, 2), 2);
+            rescaled_shared_blobs_sizes = rescaled_shared_blobs_w_lims(all_shared_blobs, 2, 1) - rescaled_shared_blobs_w_lims(all_shared_blobs, 1, 1) + 1;
+
+            % And now let's extend those blobs that need it!
+            for ii_b = 1:numel(all_shared_blobs)
+                % if they changed, let's enlarge them
+                if (to_be_changed(ii_b))
+                    proj_shift = rescaled_shared_blobs_w_lims(all_shared_blobs(ii_b), 1, 2) ...
+                        - rescaled_shared_blobs_w_lims(all_shared_blobs(ii_b), 1, 1);
+
+                    curr_blob_size = size(blobs{all_shared_blobs(ii_b)});
+                    new_size = [curr_blob_size(1), rescaled_shared_blobs_sizes(ii_b), curr_blob_size(3)];
+
+                    % Creating the new blob
+                    new_blob = zeros(new_size, 'like', blobs{all_shared_blobs(ii_b)});
+                    new_blob = gtPlaceSubVolume(new_blob, blobs{all_shared_blobs(ii_b)}, [0, proj_shift, 0]);
+                    blobs{all_shared_blobs(ii_b)} = new_blob;
+
+                    % Fixing the shift in parent's projection
+                    if (~isempty(geometries_ss))
+                        for ii_o = 1:numel(offsets_ss)
+                            for ii_ss = 1:numel(offsets_ss{ii_o})
+                                inds_of_shared_blobs = offsets_ss{ii_o}{ii_ss}.blob_offsets == all_shared_blobs(ii_b);
+
+                                offsets_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) ...
+                                    = offsets_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) + proj_shift;
+                            end
+                        end
+                    else
+                        for ii_o = 1:numel(offsets)
+                            inds_of_shared_blobs = offsets{ii_o}.blob_offsets == all_shared_blobs(ii_b);
 
-                            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)));
+                            offsets{ii_o}.blob_offsets(inds_of_shared_blobs) ...
+                                = offsets{ii_o}.blob_offsets(inds_of_shared_blobs) + proj_shift;
                         end
                     end
-                    offsets_ss = cat(1, offsets_ss, orient_offsets_ss);
-                else
-                    orient_offsets = algo_params(ii_g).offsets;
-                    for ii_o = 1:numel(orient_offsets)
-                        to_be_summed = ~shared(orient_offsets{ii_o}.blob_offsets);
+                end
+            end
 
-                        orient_offsets{ii_o}.blob_offsets(to_be_summed) ...
-                            = orient_offsets{ii_o}.blob_offsets(to_be_summed) + base_offset_shift;
+            base_offset_shift = num_parent_blobs;
+            for ii_o_reg = 2:num_regions
+                curr_or = sampler(ii_o_reg).get_reference_grain();
+                curr_sel = sampler(ii_o_reg).selected;
+                % This contains the position of the shared blobs in
+                % parent's numbering
+                shared_pos = curr_or.proj(self.det_index).shared_bls_parent(curr_sel);
+                % This indicates the ones that are actually shared among
+                % the twin's blobs
+                shared = shared_pos > 0;
+
+                proj_shifts = rescaled_shared_blobs_w_lims(shared_pos(shared), 1, ii_o_reg+1) ...
+                    - rescaled_shared_blobs_w_lims(shared_pos(shared), 1, 1);
+                blobs_with_proj_to_be_shifted = find(proj_shifts > 0);
 
-                        orient_offsets{ii_o}.blob_offsets(~to_be_summed) ...
-                            = included_pos_in_selected(shared_pos(orient_offsets{ii_o}.blob_offsets(~to_be_summed)));
+                % Shifting projection to the concatenated blobs, unless that
+                % blob was shared between parent and twin
+                if (~isempty(geometries_ss))
+                    proj_coeffs_ss = algo_params(ii_o_reg).offsets_ss;
+                    for ii_o = 1:numel(proj_coeffs_ss)
+                        for ii_ss = 1:numel(proj_coeffs_ss{ii_o})
+                            % Indices of the blobs, that we project to, in
+                            % curr_or
+                            inds_of_blobs = proj_coeffs_ss{ii_o}{ii_ss}.blob_offsets;
+
+                            % to_be_summed contains the indices to the
+                            % blobs belonging only to the current region
+                            % identified by temp_gr
+                            to_be_summed = ~shared(inds_of_blobs);
+
+                            % The indices of the given blobs belonging to
+                            % the current orientation will be shifted by
+                            % the number of blobs accumulated already
+                            inds_of_blobs(to_be_summed) ...
+                                = inds_of_blobs(to_be_summed) + base_offset_shift;
+
+                            % Treating the shared blobs
+                            inds_of_blobs(~to_be_summed) ...
+                                = included_pos_in_selected(shared_pos(inds_of_blobs(~to_be_summed)));
+
+                            proj_coeffs_ss{ii_o}{ii_ss}.blob_offsets = inds_of_blobs;
+
+                            % And now taking care of the positions inside
+                            % those blobs
+                            for ii_b = 1:numel(blobs_with_proj_to_be_shifted)
+                                inds_of_shared_blobs = inds_of_blobs == shared_pos(shared(blobs_with_proj_to_be_shifted(ii_b)));
+
+                                proj_coeffs_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) = ...
+                                    proj_coeffs_ss{ii_o}{ii_ss}.proj_offsets(inds_of_shared_blobs) + proj_shifts(blobs_with_proj_to_be_shifted(ii_b));
+                            end
+                        end
+                    end
+                    offsets_ss = cat(1, offsets_ss, proj_coeffs_ss);
+                else
+                    proj_coeffs = algo_params(ii_o_reg).offsets;
+                    for ii_o = 1:numel(proj_coeffs)
+                        % Indices of the blobs, that we project to, in
+                        % curr_or
+                        inds_of_blobs = proj_coeffs{ii_o}.blob_offsets;
+
+                        % to_be_summed contains the indices to the
+                        % blobs belonging only to the current region
+                        % identified by temp_gr
+                        to_be_summed = ~shared(inds_of_blobs);
+
+                        % The indices of the given blobs belonging to
+                        % the current orientation will be shifted by
+                        % the number of blobs accumulated already
+                        inds_of_blobs(to_be_summed) ...
+                            = inds_of_blobs(to_be_summed) + base_offset_shift;
+
+                        % Treating the shared blobs
+                        inds_of_blobs(~to_be_summed) ...
+                            = included_pos_in_selected(shared_pos(inds_of_blobs(~to_be_summed)));
+
+                        proj_coeffs{ii_o}.blob_offsets = inds_of_blobs;
+
+                        % And now taking care of the positions inside
+                        % those blobs
+                        for ii_b = 1:numel(blobs_with_proj_to_be_shifted)
+                            inds_of_shared_blobs = inds_of_blobs == shared_pos(shared(blobs_with_proj_to_be_shifted(ii_b)));
+
+                            proj_coeffs{ii_o}.proj_offsets(inds_of_shared_blobs) = ...
+                                proj_coeffs{ii_o}.proj_offsets(inds_of_shared_blobs) + proj_shifts(blobs_with_proj_to_be_shifted(ii_b));
+                        end
                     end
-                    offsets = cat(1, offsets, orient_offsets);
+                    offsets = cat(1, offsets, proj_coeffs);
                 end
 
-                base_offset_shift = base_offset_shift + numel(algo_params(ii_g).blobs);
+                base_offset_shift = base_offset_shift + numel(algo_params(ii_o_reg).blobs);
             end
 
             algo = Gt6DBlobReconstructor(volume_size, blobs, ...
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
index 17ce5df7..c15470ce 100644
--- a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
+++ b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
@@ -14,6 +14,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
 
     conf = struct( ...
         'verbose', false, ...
+        'det_ind', 1, ...
         'min_eta', 15, ...
         'rspace_oversize', 1.1, ...
         'check_spots', false, ...
@@ -25,8 +26,8 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
     symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
 
     if (numel(parent_id) > 1)
-        error('gt6DCreateProjDataFromTwinnedGrainTheoFwdProj:wrong_argument', ...
-            'Only one grain should be passed')
+        error('gt6DCreateProjDataFromTwinnedGrainFwdProj:wrong_argument', ...
+            'Only one grain should be passed as parent')
     end
     if (~isstruct(parent_id))
         grs(1) = gtLoadGrain(phase_id, parent_id);
@@ -107,15 +108,11 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
     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
+    refgr_proj = refgr.proj(conf.det_ind);
+
+    ref_ondet = refgr_proj.ondet;
+    ref_included = refgr_proj.included;
+    ref_selected = refgr_proj.selected;
 
     if (conf.verbose && false)
         try
@@ -184,10 +181,14 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
             ]');
     end
 
-    vol_half_size = [refgr.proj.vol_size_y, refgr.proj.vol_size_x, refgr.proj.vol_size_z] / 2;
+    vol_size = [ ...
+        refgr_proj.vol_size_y, ...
+        refgr_proj.vol_size_x, ...
+        refgr_proj.vol_size_z ];
+    vol_half_size = vol_size / 2;
     space_bbox = [ ...
-            refgr.proj.centerpix - vol_half_size * conf.rspace_oversize, ...
-            refgr.proj.centerpix + vol_half_size * conf.rspace_oversize ];
+        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);
 
@@ -250,15 +251,11 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
         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
+            twingr_proj = twingr.proj(conf.det_ind);
+
+            twin_ondet = twingr_proj.ondet;
+            twin_inc = twingr_proj.included;
+            twin_sel = twingr_proj.selected;
 
             bool_twin_ondet = false(size(refgr.allblobs.omega));
             bool_twin_ondet(twin_ondet) = true;
@@ -316,16 +313,17 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
             local_selected = ref_selected;
         end
 
+        num_blobs = numel(local_included);
+
         if (ii_g == 1)
             blobs = grs(ii_g).proj.bl;
         else
-            blobs = gtFwdSimBlobDefinition();
+            blobs = gtFwdSimBlobDefinition('blob', num_blobs);
             blobs(~shared_parent) = grs(ii_g).proj.bl(~shared_parent_inc_redundant);
             % Using the same because it gives the same size/uvw-limits
             blobs(shared_parent) = samp_ors(1).proj.bl(shared_parent_pos);
         end
 
-        num_blobs = numel(local_included);
         for ii_b = 1:num_blobs
             num_chars = fprintf('%02d/%02d', ii_b, num_blobs);
 
@@ -353,11 +351,15 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
         proj.selected = local_selected;
 
         if (ii_g > 1)
-            proj.shared_parent = zeros(size(shared_parent));
-            proj.shared_parent(shared_parent) = find(shared_parent_pos);
+            proj.shared_bls_parent = zeros(size(shared_parent));
+            proj.shared_bls_parent(shared_parent) = find(shared_parent_pos);
+
+            samp_ors(ii_g).proj(conf.det_ind).shared_bls_twins(shared_parent_pos) = true;
+        else
+            proj.shared_bls_twins = false(num_blobs, 1);
         end
 
-        samp_ors(ii_g).proj = proj;
+        samp_ors(ii_g).proj(conf.det_ind) = proj;
 
         if (conf.check_spots && ii_g > 1)
             samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true);
@@ -396,7 +398,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
 end
 
 function orients = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p)
-    orients = {};
+    orients = cell(2, 2, 2, 2, 2, 2);
     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
@@ -412,9 +414,9 @@ function orients = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim
 
 %                             fprintf(' center: %f, %f, %f, r_vec: %f, %f, %f\n', base_x3, base_r3)
 
-                            orients{end+1} = struct( ...
+                            orients{ii_r3+1, ii_r2+1, ii_r1+1, ii_x3+1, ii_x2+1, ii_x1+1} = struct( ...
                                 'id', refgr.id, 'phaseid', refgr.phaseid, ...
-                                'center', base_x3, 'R_vector', base_r3); %#ok<AGROW>
+                                'center', base_x3, 'R_vector', base_r3);
                         end
                     end
                 end
-- 
GitLab