diff --git a/zUtil_Deformation/gt6DCalculateFullSpaceCorners.m b/zUtil_Deformation/gt6DCalculateFullSpaceCorners.m
new file mode 100644
index 0000000000000000000000000000000000000000..f986ac7bb38b8438d437b552b007f938e8220b92
--- /dev/null
+++ b/zUtil_Deformation/gt6DCalculateFullSpaceCorners.m
@@ -0,0 +1,30 @@
+function orients = gt6DCalculateFullSpaceCorners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p)
+    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
+            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{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);
+                        end
+                    end
+                end
+            end
+        end
+    end
+
+    orients = gtCalculateGrain_p(orients, p, 'ref_omind', refgr.allblobs.omind);
+    orients = [orients{:}];
+end
\ No newline at end of file
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromExtendedGrain.m b/zUtil_Deformation/gt6DCreateProjDataFromExtendedGrain.m
index fa3ed314a6e88d7edd1043e07633154397ac6377..f815a01273ec63aff10bee8ac4e81762bcf35e90 100644
--- a/zUtil_Deformation/gt6DCreateProjDataFromExtendedGrain.m
+++ b/zUtil_Deformation/gt6DCreateProjDataFromExtendedGrain.m
@@ -13,6 +13,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     conf = struct( ...
         'verbose', false, ...
         'min_eta', 20, ...
+        'det_index', 1, ...
         'ospace_oversize', 1.1, ...
         'rspace_oversize', 1.1, ...
         'ospace_lims', [], ...
@@ -31,20 +32,10 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     p = gtLoadParameters();
 
     refgr = gr;
-    refgr_omind = refgr.allblobs.omind;
-    refgr_omind_ind = [ ...
-        find(refgr_omind == 1); find(refgr_omind == 2); ...
-        find(refgr_omind == 3); find(refgr_omind == 4) ];
-
-    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
+
+    ref_ondet = refgr.proj(conf.det_index).ondet;
+    ref_included = refgr.proj(conf.det_index).included;
+    ref_selected = refgr.proj(conf.det_index).selected;
 
     fwdsim_sel_refl = ref_ondet(ref_included(ref_selected));
     fwdsim_sel_refl_logic = false(size(ref_ondet));
@@ -58,12 +49,8 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
         used_refl = fwdsim_sel_refl;
     end
 
-    vol_half_size = [gr.proj.vol_size_y, gr.proj.vol_size_x, gr.proj.vol_size_z] / 2;
-    space_bbox = [ ...
-            gr.proj.centerpix - vol_half_size * conf.rspace_oversize, ...
-            gr.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);
+    [gr_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(gr, conf.det_index);
+    gr_center_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo, p.samgeo, false, false);
 
     estim_space_bbox_mm = [ ...
         gtGeoSam2Sam(estim_space_bbox_pix(1:3), p.recgeo, p.samgeo, false, false), ...
@@ -71,7 +58,8 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3);
 
     if (isempty(conf.ospace_lims))
-        sampler = GtOrientationSampling(p, gr, 'verbose', conf.verbose);
+        sampler = GtOrientationSampling(p, gr, ...
+            'detector_index', conf.det_index, '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);
@@ -85,44 +73,17 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     end
     bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
 
+    gr_center_rod = (estim_orient_bbox(4:6) + estim_orient_bbox(1:3)) / 2;
+
     bbox_size_deg = 2 * atand(bbox_size_rod);
 
     % 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 = {};
-    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)
-
-                            or{end+1} = struct( ...
-                                'id', refgr.id, 'phaseid', refgr.phaseid, ...
-                                'center', base_x3, 'R_vector', base_r3); %#ok<AGROW>
-                        end
-                    end
-                end
-            end
-        end
-    end
-    or = gtCalculateGrain_p(or, p);
-    or = [or{:}];
-
-    num_ors = numel(or);
-
-    for ii_g = num_ors:-1:1
-        or(ii_g) = gtGrainAllblobsFilterOrder(or(ii_g), refgr_omind_ind);
-    end
+    ors = gt6DCalculateFullSpaceCorners( ...
+        estim_space_bbox_mm, bbox_size_mm, ...
+        estim_orient_bbox, bbox_size_rod, ...
+        refgr, p);
 
     if (conf.verbose)
         fprintf('\n');
@@ -134,11 +95,11 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
         fprintf('\n');
     end
 
-    or_abs = cat(1, or(:).allblobs);
+    or_abs = cat(1, ors(:).allblobs);
     uvw_tab = zeros(numel(used_refl), num_ors, 3);
     for ii_g = 1:num_ors
         if (isfield(or_abs(1), 'detector'))
-            uvw_tab(:, ii_g, :) = or_abs(ii_g).detector(1).uvw(used_refl, :);
+            uvw_tab(:, ii_g, :) = or_abs(ii_g).detector(conf.det_index).uvw(used_refl, :);
         else
             uvw_tab(:, ii_g, :) = or_abs(ii_g).uvw(used_refl, :);
         end
@@ -146,25 +107,13 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
 
     refor = struct( ...
         'id', refgr.id, 'phaseid', refgr.phaseid, ...
-        '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 = gtGrainAllblobsFilterOrder(refor, refgr_omind_ind);
-
-    num_images = gtGetTotNumberOfImages(p);
-    refor_ws = refor.allblobs.omega(used_refl) / p.labgeo.omstep;
+        'center', gr_center_mm, 'R_vector', gr_center_rod );
+    refor = gtCalculateGrain(refor, p, 'ref_omind', refgr.allblobs.omind);
 
     % 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, :));
+    refor_ws = refor.allblobs.omega(used_refl) / gtGetOmegaStepDeg(p, conf.det_index);
+    uvw_tab(:, :, 3) = gtGrainAnglesTabularFix360deg(uvw_tab(:, :, 3), refor_ws, p);
 
     img_bboxes_initial = [
         squeeze(floor(min(uvw_tab, [], 2))), ...
@@ -214,7 +163,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     end
     sel_reflections = sel_reflections(~inconvenient_etas);
 
-    blobs(1:num_blobs) = gtFwdSimBlobDefinition();
+    blobs = gtFwdSimBlobDefinition('blob', num_blobs);
 
     num_sel_refl = numel(sel_reflections);
     for ii_i = 1:num_sel_refl
@@ -223,7 +172,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
         ii_b = sel_reflections(ii_i);
 
         bb = [img_bboxes(ii_i, 1:2), img_sizes(ii_i, 1:2)];
-        blob_vol = gtGetRawRoi(img_bboxes(ii_i, 3), img_bboxes(ii_i, 6), p.acq, bb);
+        blob_vol = gtGetRawRoi(img_bboxes(ii_i, 3), img_bboxes(ii_i, 6), p.acq, bb, conf.det_index);
         blob_vol(blob_vol < 0) = 0;
         blob_vol(isnan(blob_vol)) = 0;
         blob_bb = [img_bboxes(ii_i, 1:3), img_sizes(ii_i, :)];
@@ -298,6 +247,8 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     end
     fprintf('Done.\n')
 
+    proj = gtFwdSimProjDefinition();
+
     proj.ondet = ref_ondet;
     if (~conf.include_all)
         proj.included = ref_included;
@@ -338,8 +289,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     spots = permute(cat(3, spots{:}), [1 3 2]);
 
     % These reconstructions usually need some padding
-    proj.centerpix = estim_space_bbox_pix(1:3) + bbox_size_pix / 2;
-%     proj.centerpix = gtGeoSam2Sam(refor.center, p.samgeo, p.recgeo, false, false);
+    proj.centerpix = gr_center_pix;
     proj.bl = blobs;
     proj.stack = spots;
     % We apply a simple oversize
@@ -348,8 +298,8 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     proj.vol_size_y = vol_size(1);
     proj.vol_size_z = vol_size(3);
 
-    refor.proj = proj;
-    refor.bb_ors = or;
+    refor.proj(conf.det_index) = proj;
+    refor.bb_ors = ors;
 
     if (conf.save)
         fprintf('Saving the extended grain file..')
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m b/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m
index d261734b3c05f93674d5846fa69d0df4eb2fa08d..1c1054b8b44d15b68489eec5875202e1f7d9afed 100644
--- a/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m
+++ b/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m
@@ -13,6 +13,8 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     conf = struct( ...
         'verbose', false, ...
         'min_eta', 15, ...
+        'det_index', 1, ...
+        'ospace_oversize', 1.1, ...
         'save', false );
     conf = parse_pv_pairs(conf, varargin);
 
@@ -51,96 +53,49 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
 
     refgr = grs(1);
     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 = num_grains:-1:2
-        grs(ii_g) = gtGrainAllblobsFilterOrder(grs(ii_g), refgr_omind_ind);
+        grs(ii_g) = gtGrainAllblobsFilterOrder(grs(ii_g), refgr_omind);
     end
 
-    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
+    ref_ondet = refgr.proj(conf.det_index).ondet;
+    ref_included = refgr.proj(conf.det_index).included;
+    ref_selected = refgr.proj(conf.det_index).selected;
+
     fwdsim_sel_refl = ref_ondet(ref_included(ref_selected));
 
-    space_bboxes = zeros(num_grains, 6);
-    for ii_g = 1:num_grains
-        vol_half_size = [grs(ii_g).proj.vol_size_y, grs(ii_g).proj.vol_size_x, grs(ii_g).proj.vol_size_z] / 2;
-        space_bboxes(ii_g, :) = [ ...
-            grs(ii_g).proj.centerpix - vol_half_size, ...
-            grs(ii_g).proj.centerpix + vol_half_size ];
-    end
-    estim_space_bbox_pix = [floor(min(space_bboxes(:, 1:3), [], 1)), ceil(max(space_bboxes(:, 4:6), [], 1))];
-    bbox_size_pix = estim_space_bbox_pix(4:6) - estim_space_bbox_pix(1:3);
+    % Real-space volume size
+    [gr_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(grs, conf.det_index);
+    gr_center_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo, p.samgeo, false, false);
 
     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);
 
+    % Orientation-space volume size
     for ii_g = num_grains:-1:1
-        sampler = GtOrientationSampling(p, grs(ii_g), 'verbose', conf.verbose);
+        sampler = GtOrientationSampling(p, grs(ii_g), ...
+            'detector_index', conf.det_index, 'verbose', conf.verbose);
         odfw_R_vectors{ii_g} = sampler.guess_ODF_BB()';
     end
     r_vecs = cat(1, odfw_R_vectors{:});
     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;
+    % oversizing the orienation BB a bit
+    delta_bbox_size_rod = bbox_size_rod * (conf.ospace_oversize - 1) / 2;
     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);
-
-    % 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 = {};
-    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)
-
-                            or{end+1} = struct( ...
-                                'id', refgr.id, 'phaseid', refgr.phaseid, ...
-                                'center', base_x3, 'R_vector', base_r3); %#ok<AGROW>
-                        end
-                    end
-                end
-            end
-        end
-    end
-    or = gtCalculateGrain_p(or, p);
-    or = [or{:}];
-
-    num_ors = numel(or);
+    gr_center_rod = (estim_orient_bbox(4:6) + estim_orient_bbox(1:3)) / 2;
 
-    for ii_g = num_ors:-1:1
-        or(ii_g) = gtGrainAllblobsFilterOrder(or(ii_g), refgr_omind_ind);
-    end
+    bbox_size_deg = 2 * atand(bbox_size_rod);
 
     if (conf.verbose)
         f = figure();
         ax = axes('parent', f);
         hold(ax, 'on');
-        gt6DPlotOrientationBBox(ax, cat(1, or(:).R_vector));
+        gt6DPlotOrientationBBox(ax, [estim_orient_bbox(1:3); estim_orient_bbox(4:6)]);
         for ii_g = 1:num_grains
             scatter3(ax, grs(ii_g).R_vector(1), grs(ii_g).R_vector(2), grs(ii_g).R_vector(3), 30);
         end
@@ -157,37 +112,30 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
         fprintf('\n');
     end
 
-    or_abs = cat(1, or(:).allblobs);
+    % 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.
+    ors = gt6DCalculateFullSpaceCorners( ...
+        estim_space_bbox_mm, bbox_size_mm, ...
+        estim_orient_bbox, bbox_size_rod, ...
+        refgr, p);
+
+    or_abs = cat(1, ors(:).allblobs);
+    num_ors = numel(ors);
     uvw_tab = zeros(numel(fwdsim_sel_refl), num_ors, 3);
     for ii_g = 1:num_ors
-        if (isfield(or_abs(1), 'detector'))
-            uvw_tab(:, ii_g, :) = or_abs(ii_g).detector(1).uvw(fwdsim_sel_refl, :);
-        else
-            uvw_tab(:, ii_g, :) = or_abs(ii_g).uvw(fwdsim_sel_refl, :);
-        end
+        uvw_tab(:, ii_g, :) = or_abs(ii_g).detector(conf.det_index).uvw(fwdsim_sel_refl, :);
     end
 
     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 = gtGrainAllblobsFilterOrder(refor, refgr_omind_ind);
-
-    num_images = gtGetTotNumberOfImages(p);
-    refor_ws = refor.allblobs.omega(fwdsim_sel_refl) / p.labgeo.omstep;
+        'center', gr_center_mm, 'R_vector', gr_center_rod );
+    refor = gtCalculateGrain(refor, p, 'ref_omind', refgr_omind);
 
     % 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, :));
+    refor_ws = refor.allblobs.omega(fwdsim_sel_refl) / gtGetOmegaStepDeg(p, conf.det_index);
+    uvw_tab(:, :, 3) = gtGrainAnglesTabularFix360deg(uvw_tab(:, :, 3), refor_ws, p);
 
     img_bboxes = [
         squeeze(floor(min(uvw_tab, [], 2))), ...
@@ -222,7 +170,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     fprintf('Loading raw images: ')
     num_blobs = numel(ref_included);
 
-    blobs(1:num_blobs) = gtFwdSimBlobDefinition();
+    blobs = gtFwdSimBlobDefinition('blob', num_blobs);
 
     sel_reflections = find(ref_selected);
     sel_reflections = sel_reflections(~inconvenient_etas);
@@ -234,7 +182,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
         ii_b = sel_reflections(ii_i);
 
         bb = [img_bboxes(ii_i, 1:2), img_sizes(ii_i, 1:2)];
-        blob_vol = gtGetRawRoi(img_bboxes(ii_i, 3), img_bboxes(ii_i, 6), p.acq, bb);
+        blob_vol = gtGetRawRoi(img_bboxes(ii_i, 3), img_bboxes(ii_i, 6), p.acq, bb, conf.det_index);
         blob_vol(blob_vol < 0) = 0;
         blob_bb = [img_bboxes(ii_i, 1:3), img_sizes(ii_i, :)];
 
@@ -318,7 +266,9 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     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 = gtFwdSimProjDefinition();
+
+    proj.centerpix = gr_center_pix;
     proj.bl = blobs;
     proj.stack = spots;
     vol_size = bbox_size_pix + mean(bbox_size_pix) * 0.3;
@@ -333,8 +283,8 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     finally_selected(sel_reflections) = true;
     proj.selected = finally_selected;
 
-    refor.proj = proj;
-    refor.bb_ors = or;
+    refor.proj(conf.det_index) = proj;
+    refor.bb_ors = ors;
 
     if (conf.save)
         fprintf('Saving the cluster file..')
diff --git a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
index c15470ce73ec6dd0eafeb54a8ad073d2dc2c4377..4d447226f998815c6f272149e9464e7695d448d2 100644
--- a/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
+++ b/zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
@@ -14,9 +14,10 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
 
     conf = struct( ...
         'verbose', false, ...
-        'det_ind', 1, ...
+        'det_index', 1, ...
         'min_eta', 15, ...
         'rspace_oversize', 1.1, ...
+        'ospace_oversize', 1.1, ...
         'check_spots', false, ...
         'stack_oversize', 1.4, ...
         'save', false );
@@ -106,9 +107,13 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
     end
 
     refgr = grs(1);
-    grs(2:end) = reorder_reflections(grs(2:end), refgr);
 
-    refgr_proj = refgr.proj(conf.det_ind);
+    refgr_omind = refgr.allblobs.omind;
+    for ii_g = 2:numel(grs)
+        grs(ii_g) = gtGrainAllblobsFilterOrder(grs(ii_g), refgr_omind); %#ok<AGROW>
+    end
+
+    refgr_proj = refgr.proj(conf.det_index);
 
     ref_ondet = refgr_proj.ondet;
     ref_included = refgr_proj.included;
@@ -162,7 +167,6 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
 
     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), ...
@@ -181,62 +185,57 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
             ]');
     end
 
-    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 ];
-    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);
+    [gr_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(grs, conf.det_index, conf.rspace_oversize);
+    gr_center_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo, p.samgeo, false, false);
 
     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);
 
+    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);
+    end
+
     img_sizes = zeros(num_grains, 2);
 
     for ii_g = num_grains:-1:1
-        sampler = GtOrientationSampling(p, grs(ii_g), 'verbose', conf.verbose);
+        sampler = GtOrientationSampling(p, grs(ii_g), ...
+            'detector_index', conf.det_index, '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;
+        delta_bbox_size_rod = bbox_size_rod * (conf.ospace_oversize - 1) / 2;
         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('%d) Estimated orientation BBox: [%3.3f, %3.3f, %3.3f] -> [%3.3f, %3.3f, %3.3f]\n', ii_g, 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);
+            'center', gr_center_mm, 'R_vector', grs(ii_g).R_vector );
+        refor = gtCalculateGrain(refor, p, 'ref_omind', refgr_omind);
 
-        refor.bb_ors = get_6D_space_corners(estim_space_bbox_mm, ...
+        refor.bb_ors = gt6DCalculateFullSpaceCorners(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), ];
+        img_sizes(ii_g, :) = [ ...
+            size(grs(ii_g).proj(conf.det_index).stack, 1), ...
+            size(grs(ii_g).proj(conf.det_index).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);
 
@@ -245,13 +244,16 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
     bool_ref_sel = false(size(refgr.allblobs.omega));
     bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;
 
+    % At the moment we only support the configuration: parent + twins
+    refl_matches = gt6DGetMatchingReflections(refgr, grs(2:end), 2);
+
     for ii_g = 1:num_grains
         fprintf('Loading raw images for grain %d: ', ii_g)
 
         if (ii_g > 1)
             twingr = grs(ii_g);
 
-            twingr_proj = twingr.proj(conf.det_ind);
+            twingr_proj = twingr.proj(conf.det_index);
 
             twin_ondet = twingr_proj.ondet;
             twin_inc = twingr_proj.included;
@@ -268,17 +270,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
             % twin, and were included for the parent
             bool_test_inc = bool_ref_inc & bool_twin_ondet;
 
-            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);
+            test_shared = refl_matches(bool_test_inc, ii_g-1);
 
             % let's add the shared as included
             shared_parent = false(size(refgr.allblobs.omega));
@@ -338,7 +330,9 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
         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 = gtFwdSimProjDefinition();
+
+        proj.centerpix = gr_center_pix;
         proj.bl = blobs;
         proj.stack = spots;
         vol_size = bbox_size_pix + mean(bbox_size_pix) * 0.3;
@@ -354,15 +348,15 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
             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;
+            samp_ors(1).proj(conf.det_index).shared_bls_twins(shared_parent_pos) = true;
         else
             proj.shared_bls_twins = false(num_blobs, 1);
         end
 
-        samp_ors(ii_g).proj(conf.det_ind) = proj;
+        samp_ors(ii_g).proj(conf.det_index) = proj;
 
         if (conf.check_spots && ii_g > 1)
-            samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true);
+            samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true, conf.det_index);
         end
     end
 
@@ -397,49 +391,6 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
     end
 end
 
-function orients = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p)
-    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
-            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{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);
-                        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 [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf)
 
     max_img_sizes = [max(img_sizes(:, 1)), max(img_sizes(:, 2))];
@@ -496,4 +447,3 @@ end
 
 
 
-
diff --git a/zUtil_Deformation/gt6DGetMatchingReflections.m b/zUtil_Deformation/gt6DGetMatchingReflections.m
new file mode 100644
index 0000000000000000000000000000000000000000..f756d2deb9c7d4a039ee7f102a9b9fbef3dc32f3
--- /dev/null
+++ b/zUtil_Deformation/gt6DGetMatchingReflections.m
@@ -0,0 +1,28 @@
+function match = gt6DGetMatchingReflections(ref_or, ors, angular_toll)
+
+    if (~exist('angular_toll', 'var') || isempty(angular_toll))
+        angular_toll = 1; % degrees
+    end
+
+    ref_ws = ref_or.allblobs.omega;
+    ref_ns = ref_or.allblobs.eta;
+    ref_Lfactors = ref_or.allblobs.lorentzfac;
+
+    w_toll = angular_toll * ref_Lfactors;
+    n_toll = angular_toll;
+
+    ors_allblobs = [ors(:).allblobs];
+    ors_ws = cat(2, ors_allblobs(:).omega);
+    ors_ns = cat(2, ors_allblobs(:).eta);
+
+    ones_num_ors = ones(numel(ors), 1);
+    ref_ws = ref_ws(:, ones_num_ors);
+    ref_ns = ref_ns(:, ones_num_ors);
+
+    w_toll = w_toll(:, ones_num_ors);
+
+    w_diff = abs(mod(ref_ws - ors_ws + 180, 360) - 180);
+    n_diff = abs(mod(ref_ns - ors_ns + 180, 360) - 180);
+
+    match = (w_diff < w_toll) & (n_diff < n_toll);
+end
diff --git a/zUtil_Deformation/gt6DMergeRealSpaceVolumes.m b/zUtil_Deformation/gt6DMergeRealSpaceVolumes.m
new file mode 100644
index 0000000000000000000000000000000000000000..36e824630748d3ed57c20768e227ac006abb59b6
--- /dev/null
+++ b/zUtil_Deformation/gt6DMergeRealSpaceVolumes.m
@@ -0,0 +1,39 @@
+function [gr_center_pix, vol_size_pix, bbox_pix] = gt6DMergeRealSpaceVolumes(grs, det_index, oversize)
+%     num_ors = numel(grs);
+%     vol_sizes = zeros(num_ors, 3);
+%     vol_centers = zeros(num_ors, 3);
+%     for ii = 1:num_ors
+%         vol_sizes(ii, :) = [ ...
+%             grs(ii).proj(det_index).vol_size_y, ...
+%             grs(ii).proj(det_index).vol_size_x, ...
+%             grs(ii).proj(det_index).vol_size_z ];
+%         vol_centers(ii, :) = grs(ii).proj(det_index).centerpix;
+%     end
+
+    projs = arrayfun(@(x){x.proj(det_index)}, grs);
+    projs = [projs{:}];
+
+    vol_sizes = [cat(1, projs.vol_size_y), cat(1, projs.vol_size_x), cat(1, projs.vol_size_z)];
+    vol_centers = cat(1, projs.centerpix);
+
+    vol_half_sizes = vol_sizes / 2;
+
+    min_corners = vol_centers - vol_half_sizes;
+    max_corners = vol_centers + vol_half_sizes;
+
+    full_vol_min_corners = min(min_corners, [], 1);
+    full_vol_max_corners = max(max_corners, [], 1);
+
+    gr_center_pix = (full_vol_max_corners + full_vol_min_corners) / 2;
+
+    vol_size_pix = full_vol_max_corners - full_vol_min_corners;
+
+    if (exist('oversize', 'var'))
+        vol_size_pix = vol_size_pix * oversize;
+    end
+
+    vol_size_pix = ceil(vol_size_pix);
+
+    vol_half_size_pix = vol_size_pix / 2;
+    bbox_pix = [gr_center_pix - vol_half_size_pix, gr_center_pix + vol_half_size_pix];
+end