diff --git a/zUtil_ForwardSim/gtForwardSimulate_v2.m b/zUtil_ForwardSim/gtForwardSimulate_v2.m
index 906b6fd087205dc01703ca00459d39bb29111a93..04babf8dc8affae1a7a3c239ecd96523374d0a96 100644
--- a/zUtil_ForwardSim/gtForwardSimulate_v2.m
+++ b/zUtil_ForwardSim/gtForwardSimulate_v2.m
@@ -234,7 +234,7 @@ for n = first : last
     fprintf(' (Done).\n - Producing blob masks..')
 
     proj_bls = predict_spot_masks_from_indexed( ...
-        gr, grain, segmentedSpots, spotsCommProps, onDet, parameters);
+        gr, grain, segmentedSpots, spotsCommProps, onDet, parameters, det_index);
 
     fprintf(' (Done).\n - Preallocating structs..')
 
@@ -687,7 +687,7 @@ function g_twins = forward_simulate_twin(gr, fwd_sim, proj, twin_params, ...
         uvw = twin.allblobs.detector(det_index).uvw;
         ondet = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
 
-        proj_bls = predict_spot_masks_from_fwdsimgrain(twin, proj, ondet, parameters);
+        proj_bls = predict_spot_masks_from_fwdsimgrain(twin, proj, ondet, parameters, det_index);
 
         [twin.fwd_sim, bl] = forward_simulate_grain(...
             twin, ondet, spotsCommProps, segmentedSpots, proj_bls, parameters, det_index);
@@ -770,7 +770,12 @@ function displayShowFsim(phaseID, grain, varargin)
 end
 
 function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
-    segmentedSpots, spotsCommProps, onDet, parameters)
+    segmentedSpots, spotsCommProps, onDet, parameters, det_ind)
+
+    samgeo = parameters.samgeo;
+    labgeo = parameters.labgeo;
+    detgeo = parameters.detgeo(det_ind);
+    omstep = gtGetOmegaStepDeg(parameters, det_ind);
 
     % First argument: included_INDXTR, might be needed in the future for
     % some reason...
@@ -785,9 +790,15 @@ function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
     BBsizes = cat(1, bls_INDXTR(:).mbbsize);
     BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
 
-    Ws = segmentedSpots.CentroidImage(difspotID_INDXTR) .* parameters.labgeo.omstep;
+    Ws = segmentedSpots.CentroidImage(difspotID_INDXTR) .* omstep;
+
+    gr_center = gr.center(ones(numel(difspotID_INDXTR), 1), :);
+
+    projvecs = [segmentedSpots.CentroidX(difspotID_INDXTR), segmentedSpots.CentroidY(difspotID_INDXTR)];
+    projvecs = gtGeoDet2Lab(projvecs, detgeo, false);
+    projvecs = gtGeoLab2Sam(projvecs, Ws, labgeo, samgeo, 0, 0) - gr_center;
 
-    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);
+    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
 
     % Producing original blobs' boundingboxes (as by original segmentation
     % info)
@@ -797,17 +808,19 @@ function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
         segmentedSpots.BoundingBoxXsize(difspotID_INDXTR), ...
         segmentedSpots.BoundingBoxYsize(difspotID_INDXTR) ];
 
-    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);
+    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
 
     for ii_b = 1:numel(proj_bls)
         proj_bls(ii_b).mbbsize = proj_bls_original(ii_b).bbsize;
         proj_bls(ii_b).mbbu = proj_bls_original(ii_b).bbuim;
         proj_bls(ii_b).mbbv = proj_bls_original(ii_b).bbvim;
         proj_bls(ii_b).mbbw = proj_bls_original(ii_b).bbwim;
+
+        proj_bls(ii_b).mask = proj_bls_original(ii_b).mask;
     end
 end
 
-function proj_bls = predict_spot_masks_from_fwdsimgrain(gr, proj, onDet, parameters)
+function proj_bls = predict_spot_masks_from_fwdsimgrain(gr, proj, onDet, parameters, det_ind)
 
     centroids = zeros(numel(proj.bl), 1);
     for ii_b = 1:numel(proj.bl)
@@ -824,7 +837,9 @@ function proj_bls = predict_spot_masks_from_fwdsimgrain(gr, proj, onDet, paramet
     BBsizes = cat(1, proj.bl(:).bbsize);
     BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
 
-    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);
+    projvecs = gr.allblobs.dvecsam(onDet, :);
+
+    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
 
     % Producing original blobs' boundingboxes (as by original segmentation
     % info)
@@ -833,13 +848,15 @@ function proj_bls = predict_spot_masks_from_fwdsimgrain(gr, proj, onDet, paramet
     BBsizes = cat(1, proj.bl(:).mbbsize);
     BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
 
-    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);
+    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
 
     for ii_b = 1:numel(proj_bls)
         proj_bls(ii_b).mbbsize = proj_bls_original(ii_b).bbsize;
         proj_bls(ii_b).mbbu = proj_bls_original(ii_b).bbuim;
         proj_bls(ii_b).mbbv = proj_bls_original(ii_b).bbvim;
         proj_bls(ii_b).mbbw = proj_bls_original(ii_b).bbwim;
+
+        proj_bls(ii_b).mask = proj_bls_original(ii_b).mask;
     end
 end
 
@@ -950,8 +967,9 @@ function ls = compute_proj_bbs_partial_likelihoods(bls, proj_bls)
         mismatch_less(ii) = sum(proj_bl_mask(:));
     end
 
-    % 3 and 2 are temporary
-    ls = (1 .* mismatch_more + 1 .* mismatch_less) ./ proj_count;
+    % 3/4 makes up for the recent change to the use of a bigger estimate of
+    % the grain volume
+    ls = (1 .* mismatch_more + 3 / 4 .* mismatch_less) ./ proj_count;
 
     ls = exp(- (ls .^ 2));
 end
diff --git a/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m b/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
index 21e9f94def2778b477d64e199f4a3e9f6b1a8580..a15ab266851de75fb01ebcfc40384b083b4a2ea7 100644
--- a/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
+++ b/zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m
@@ -6,6 +6,7 @@ function proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb,
     end
 
     acq = parameters.acq;
+    fsim = parameters.fsim;
     labgeo = parameters.labgeo;
     samgeo = parameters.samgeo;
     recgeo = parameters.recgeo(det_ind);
@@ -54,11 +55,24 @@ function proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb,
     proj.num_rows  = stackVSize;
     proj.num_cols  = stackUSize;
 
-    proj.centerpix = gtGeoLab2Sam(gr_center, eye(3), labgeo, recgeo, true);
+    proj.centerpix = gtGeoSam2Sam(gr_center, samgeo, recgeo, true);
 
     use_polyhedron = numel(find(selected)) >= 3;
-    vol_size = gtFwdSimPredictVolumeSize(gr_center, omegas(selected), ...
-        bb(selected, :), parameters, stackUSize, stackVSize, use_polyhedron);
+    if (use_polyhedron)
+        % This should behave better with vertical detector
+        % (if no bug is introduced with it :D)
+        verts = gtFwdSimComputeCircumscribingPolyhedron(...
+            gr_center, diff_beam_dirs(selected, :), omegas(selected), ...
+            bb(selected, :), parameters, det_ind);
+
+        vol_size = 2 * max(abs(max(verts, [], 1)), abs(min(verts, [], 1)));
+    else
+        % We should in the future handle properly vertical detector
+        % (general geometry) maybe determining a convex shape of the grain!
+        vol_size = [stackUSize, stackUSize, stackVSize] / fsim.oversize;
+    end
+
+    vol_size = round(vol_size * fsim.oversizeVol);
 
     proj.vol_size_x = vol_size(2);
     proj.vol_size_y = vol_size(1);
diff --git a/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m b/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m
index c1fa7e1b7c9518f697a1cc95408f61cb6a375641..5ef68c3f6c5b4dea0ae4d015106947753da5f729 100644
--- a/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m
+++ b/zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m
@@ -1,19 +1,31 @@
-function verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, omegas, bb, parameters, det_index, display_figure)
-% verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, omegas, bb, parameters, det_index, display_figure)
+function verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, projvecs, omegas, bb, parameters, det_index, criterion, display_figure)
+% verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, d_beam_dirs, omegas, bb, parameters, det_index, display_figure)
 %
 % bb is [x_origin, y_origin, x_size, y_size]
 %
 % We can think of spot U/V edges like planes that cut the space in a
 % polyhedron when they intersect, having the grain volume on the inside.
+%
+% Input:
+%   gr_center : <1 x 3> vector in sample coordinates
+%   projvecs : <n x 3> vector in sample coordinates
+%   omegas : <n> vector
+%   bb : <n x 4> vector in detector coordinates
+%   parameters : dataset parameters
+%
 
     if (~exist('det_index', 'var'))
         det_index = 1;
     end
+    if (~exist('criterion', 'var'))
+        criterion = 'convhull';
+    end
     if (~exist('display_figure', 'var'))
         display_figure = false;
     end
 
     labgeo = parameters.labgeo;
+    samgeo = parameters.samgeo;
     recgeo = parameters.recgeo(det_index);
     detgeo = parameters.detgeo(det_index);
 
@@ -21,56 +33,125 @@ function verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, omegas, bb,
     rot_comp_w = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
     rot_omega = gtMathsRotationTensor(omegas, rot_comp_w);
 
-    gc_pos = gr_center(ones(numel(omegas), 1), :);
+    gc_pos_rec = gtGeoSam2Sam(gr_center, samgeo, recgeo, true);
+    gc_pos_rec = gc_pos_rec(ones(numel(omegas), 1), :);
 
-    % vectors bu and bv, in positions left, right, bottom, top
-    % bu are the plane normals in the 'u' directions, while bv are the
-    % plane normals in 'v' directions.
+    % Taking the four corners: ctl = corner top left, cbr = corner bottom
+    % right, ...
     % Here we compute their positions in the lab coordinates
-    bu_l_pos_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)), (bb(:, 2) + bb(:, 4)/2)], detgeo, false);
-    bu_r_pos_lab = gtGeoDet2Lab([bb(:, 1), (bb(:, 2) + bb(:, 4)/2)], detgeo, false);
-    bv_b_pos_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)/2), (bb(:, 2) + bb(:, 4))], detgeo, false);
-    bv_t_pos_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)/2), bb(:, 2)], detgeo, false);
+    ctl_pos_lab = gtGeoDet2Lab([bb(:, 1), bb(:, 2)], detgeo, false);
+    ctr_pos_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)), bb(:, 2)], detgeo, false);
+    cbl_pos_lab = gtGeoDet2Lab([bb(:, 1), (bb(:, 2) + bb(:, 4))], detgeo, false);
+    cbr_pos_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)), (bb(:, 2) + bb(:, 4))], detgeo, false);
+
+    % positions in the reconstruction coordinates
+    ctl_pos_rec = gtGeoLab2Sam(ctl_pos_lab, rot_omega, labgeo, recgeo, 0, 0);
+    ctr_pos_rec = gtGeoLab2Sam(ctr_pos_lab, rot_omega, labgeo, recgeo, 0, 0);
+    cbl_pos_rec = gtGeoLab2Sam(cbl_pos_lab, rot_omega, labgeo, recgeo, 0, 0);
+    cbr_pos_rec = gtGeoLab2Sam(cbr_pos_lab, rot_omega, labgeo, recgeo, 0, 0);
+
+    projvecs_rec = gtGeoSam2Sam(projvecs, samgeo, recgeo, true);
+    projvecs_rec = gtMathsNormalizeVectorsList(projvecs_rec);
 
-    bb_cent_lab = gtGeoDet2Lab([(bb(:, 1) + bb(:, 3)/2), (bb(:, 2) + bb(:, 4)/2)], detgeo, false);
+    % computing now the closest points on the projvec lines from the four
+    % corners to the grain center
+    dot_ctl_gc = sum((gc_pos_rec - ctl_pos_rec) .* projvecs_rec, 2);
+    dot_ctr_gc = sum((gc_pos_rec - ctr_pos_rec) .* projvecs_rec, 2);
+    dot_cbl_gc = sum((gc_pos_rec - cbl_pos_rec) .* projvecs_rec, 2);
+    dot_cbr_gc = sum((gc_pos_rec - cbr_pos_rec) .* projvecs_rec, 2);
 
-    % bus and bvs positions in the reconstruction coordinates
-    bu_l_pos_rec  = gtGeoLab2Sam(bu_l_pos_lab, rot_omega, labgeo, recgeo, 0, 0) - gc_pos;
-    bu_r_pos_rec  = gtGeoLab2Sam(bu_r_pos_lab, rot_omega, labgeo, recgeo, 0, 0) - gc_pos;
-    bv_b_pos_rec  = gtGeoLab2Sam(bv_b_pos_lab, rot_omega, labgeo, recgeo, 0, 0) - gc_pos;
-    bv_t_pos_rec  = gtGeoLab2Sam(bv_t_pos_lab, rot_omega, labgeo, recgeo, 0, 0) - gc_pos;
+    ctl_pos_rec = ctl_pos_rec + projvecs_rec .* dot_ctl_gc(:, [1 1 1]);
+    ctr_pos_rec = ctr_pos_rec + projvecs_rec .* dot_ctr_gc(:, [1 1 1]);
+    cbl_pos_rec = cbl_pos_rec + projvecs_rec .* dot_cbl_gc(:, [1 1 1]);
+    cbr_pos_rec = cbr_pos_rec + projvecs_rec .* dot_cbr_gc(:, [1 1 1]);
 
-    bb_cent_rec = gtGeoLab2Sam(bb_cent_lab, rot_omega, labgeo, recgeo, 0, 0) - gc_pos;
+    % Recentering to the grain center position
+    ctl_pos_rec = ctl_pos_rec - gc_pos_rec;
+    ctr_pos_rec = ctr_pos_rec - gc_pos_rec;
+    cbl_pos_rec = cbl_pos_rec - gc_pos_rec;
+    cbr_pos_rec = cbr_pos_rec - gc_pos_rec;
+
+    % computing plane normals
+    pv_t_rec = gtMathsCross(projvecs_rec, ctr_pos_rec - ctl_pos_rec);
+    pv_b_rec = gtMathsCross(projvecs_rec, cbr_pos_rec - cbl_pos_rec);
+    pu_l_rec = gtMathsCross(projvecs_rec, ctr_pos_rec - cbr_pos_rec);
+    pu_r_rec = gtMathsCross(projvecs_rec, ctl_pos_rec - cbl_pos_rec);
 
-    % Let's finally find the normalized plane normals for the boundaries
-    pv_t_rec = gtMathsCross(bb_cent_rec, bu_l_pos_rec);
     pv_t_rec = gtMathsNormalizeVectorsList(pv_t_rec);
-    pu_l_rec = gtMathsCross(bb_cent_rec, bv_b_pos_rec);
-    pu_l_rec = gtMathsNormalizeVectorsList(pu_l_rec);
-    pv_b_rec = gtMathsCross(bb_cent_rec, bu_r_pos_rec);
     pv_b_rec = gtMathsNormalizeVectorsList(pv_b_rec);
-    pu_r_rec = gtMathsCross(bb_cent_rec, bv_t_pos_rec);
+    pu_l_rec = gtMathsNormalizeVectorsList(pu_l_rec);
     pu_r_rec = gtMathsNormalizeVectorsList(pu_r_rec);
 
+    % Computing displacements
+    pv_t_rec_norm = sum(ctr_pos_rec .* pv_t_rec, 2);
+    pv_b_rec_norm = sum(cbl_pos_rec .* pv_b_rec, 2);
+    pu_l_rec_norm = sum(ctl_pos_rec .* pu_l_rec, 2);
+    pu_r_rec_norm = sum(cbr_pos_rec .* pu_r_rec, 2);
+
     % We now compute the norm of the vector from the grain_center to
     % the closest point on the plane normals
-    pv_t_rec_norm = sum(pv_t_rec .* bv_t_pos_rec, 2);
     pv_t_rec = pv_t_rec_norm(:, [1 1 1]) .* pv_t_rec;
-    pv_b_rec_norm = sum(pv_b_rec .* bv_b_pos_rec, 2);
     pv_b_rec = pv_b_rec_norm(:, [1 1 1]) .* pv_b_rec;
-    pu_l_rec_norm = sum(pu_l_rec .* bu_l_pos_rec, 2);
     pu_l_rec = pu_l_rec_norm(:, [1 1 1]) .* pu_l_rec;
-    pu_r_rec_norm = sum(pu_r_rec .* bu_r_pos_rec, 2);
     pu_r_rec = pu_r_rec_norm(:, [1 1 1]) .* pu_r_rec;
 
-    verts = gtMathsGetPolyhedronVerticesFromPlaneNormals( ...
-        [pv_t_rec; pv_b_rec; pu_l_rec; pu_r_rec] );
+    switch (criterion)
+        case 'inscribed_polyhedron'
+            verts = gtMathsGetPolyhedronVerticesFromPlaneNormals( ...
+                [pv_t_rec; pv_b_rec; pu_l_rec; pu_r_rec] );
+        case 'circumscribed_polyhedron'
+            tot_vecs = [pv_t_rec; pv_b_rec; pu_l_rec; pu_r_rec];
+            k = convhull(tot_vecs);
+            verts = gtMathsGetPolyhedronVerticesFromPlaneNormals(tot_vecs(k, :));
+        case 'convhull'
+            tot_points = [ctl_pos_rec; ctr_pos_rec; cbl_pos_rec; cbr_pos_rec];
+            k = convhull(tot_points);
+            verts = tot_points(k, :);
+        case 'box'
+            tot_points = [ctl_pos_rec; ctr_pos_rec; cbl_pos_rec; cbr_pos_rec];
+            min_points = min(tot_points, [], 1);
+            max_points = max(tot_points, [], 1);
+            verts = [min_points; ...
+                [min_points(1:2), max_points(3)]; ...
+                [min_points(1), max_points(2), min_points(3)]; ...
+                [min_points(1), max_points(2:3)]; ...
+                [max_points(1), min_points(2:3)]; ...
+                [max_points(1), min_points(2), max_points(3)]; ...
+                [max_points(1:2), min_points(3)]; ...
+                max_points ];
+    end
 
     if (display_figure)
         f = figure();
         ax = axes('parent', f);
-        hold(ax, 'on')
+        hold(ax, 'on');
+        grid(ax, 'on');
         scatter3(ax, verts(:, 1), verts(:, 2), verts(:, 3), 30, 'b', 'filled')
-        hold(ax, 'off')
+        hold(ax, 'off');
+
+        f = figure();
+        ax = axes('parent', f);
+        hold(ax, 'on');
+        grid(ax, 'on');
+        scatter3(ax, 0, 0, 0, 30, 'r', 'filled');
+        for ii = 1:numel(omegas)
+            quiver3(ax, 0, 0, 0, projvecs_rec(ii, 1), projvecs_rec(ii, 2), projvecs_rec(ii, 3));
+
+            scatter3(ax, ctl_pos_rec(ii, 1), ctl_pos_rec(ii, 2), ctl_pos_rec(ii, 3), 20, 'b');
+            scatter3(ax, ctr_pos_rec(ii, 1), ctr_pos_rec(ii, 2), ctr_pos_rec(ii, 3), 20, 'b');
+
+            quiver3(ax, 0, 0, 0, pv_t_rec(ii, 1), pv_t_rec(ii, 2), pv_t_rec(ii, 3));
+%             pause
+
+            scatter3(ax, cbl_pos_rec(ii, 1), cbl_pos_rec(ii, 2), cbl_pos_rec(ii, 3), 20, 'b');
+            scatter3(ax, cbr_pos_rec(ii, 1), cbr_pos_rec(ii, 2), cbr_pos_rec(ii, 3), 20, 'b');
+
+            quiver3(ax, 0, 0, 0, pv_b_rec(ii, 1), pv_b_rec(ii, 2), pv_b_rec(ii, 3));
+%             pause
+
+            quiver3(ax, 0, 0, 0, pu_l_rec(ii, 1), pu_l_rec(ii, 2), pu_l_rec(ii, 3));
+            quiver3(ax, 0, 0, 0, pu_r_rec(ii, 1), pu_r_rec(ii, 2), pu_r_rec(ii, 3));
+        end
+        hold(ax, 'off');
     end
 end
\ No newline at end of file
diff --git a/zUtil_ForwardSim/gtFwdSimPredictProjectedGrainBBox.m b/zUtil_ForwardSim/gtFwdSimPredictProjectedGrainBBox.m
index 8dcaeaee11b049976a8abedaccb02bfffc1e2b03..81c11dfcd746890cdaf86ba037c92bf31ea51c45 100644
--- a/zUtil_ForwardSim/gtFwdSimPredictProjectedGrainBBox.m
+++ b/zUtil_ForwardSim/gtFwdSimPredictProjectedGrainBBox.m
@@ -1,4 +1,4 @@
-function proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, sel_bb, sel_omegas, ondet, parameters, det_ind)
+function proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, sel_projvecs, sel_bb, sel_omegas, ondet, parameters, det_ind)
     if (~exist('det_ind', 'var') || isempty(det_ind))
         det_ind = 1;
     end
@@ -7,8 +7,10 @@ function proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, sel_bb, sel_omegas, on
     samgeo = parameters.samgeo;
     recgeo = parameters.recgeo(det_ind);
 
+    omstep = gtGetOmegaStepDeg(parameters, det_ind);
+
     % We first compute the bbox of the grain volume
-    verts_rec = gtFwdSimComputeCircumscribingPolyhedron(gr.center, sel_omegas, sel_bb, parameters);
+    verts_rec = gtFwdSimComputeCircumscribingPolyhedron(gr.center, sel_projvecs, sel_omegas, sel_bb, parameters, det_ind, 'convhull');
 
     % Now we project the vertices to the spots and take the convex hull
     dvecs_lab = gr.allblobs.dvec(ondet, :);
@@ -30,15 +32,16 @@ function proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, sel_bb, sel_omegas, on
     rot_omega = gtMathsRotationTensor(-omegas, rot_comp_w);
 
     ones_verts = ones(size(verts_rec, 1), 1);
-    gc_sam_pos = gr.center(ones_verts, :);
+    gc_sam = gr.center(ones_verts, :);
+
+    verts_sam = gtGeoSam2Sam(verts_rec, recgeo, samgeo, false) + gc_sam;
 
     for ii = 1:num_omegas
         rotw = rot_omega(:, :, ii);
-        gc_lab_pos = gtGeoSam2Lab(gc_sam_pos, rotw, labgeo, samgeo, false);
-        verts_lab_pos = gtGeoSam2Lab(verts_rec, rotw, labgeo, recgeo, false) + gc_lab_pos;
+        verts_lab = gtGeoSam2Lab(verts_sam, rotw, labgeo, samgeo, false);
         dvec_vers = dvecs_lab(ii .* ones_verts, :);
 
-        verts_det_pos = gtFedPredictUVMultiple([], dvec_vers', verts_lab_pos', ...
+        verts_det_pos = gtFedPredictUVMultiple([], dvec_vers', verts_lab', ...
             detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
             [detgeo.detrefu, detgeo.detrefv]')';
 
@@ -67,7 +70,7 @@ function proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, sel_bb, sel_omegas, on
 
         proj_bls(ii).bbuim = [bb_2d(1), (bb_2d(1) + bb_2d(3) - 1)];
         proj_bls(ii).bbvim = [bb_2d(2), (bb_2d(2) + bb_2d(4) - 1)];
-        proj_bls(ii).bbwim = round([omegas(ii), omegas(ii)] / labgeo.omstep);
+        proj_bls(ii).bbwim = round([omegas(ii), omegas(ii)] / omstep);
         proj_bls(ii).bbsize = [bb_2d(3:4), 1];
         proj_bls(ii).mask = mask;
 
diff --git a/zUtil_ForwardSim/gtFwdSimPredictVolumeSize.m b/zUtil_ForwardSim/gtFwdSimPredictVolumeSize.m
deleted file mode 100644
index 18a7d36cdb983aa5b932caf07bf2b4144d8e68c2..0000000000000000000000000000000000000000
--- a/zUtil_ForwardSim/gtFwdSimPredictVolumeSize.m
+++ /dev/null
@@ -1,24 +0,0 @@
-function vol_size = gtFwdSimPredictVolumeSize(gr_center, omegas, bb, parameters, stackUSize, stackVSize, use_polyhedron)
-    if (~exist('use_polyhedron', 'var') || isempty(use_polyhedron))
-        use_polyhedron = false;
-    end
-
-    fsim = parameters.fsim;
-
-    if (use_polyhedron)
-        % This should behave better with vertical detector
-        % (if no bug is introduced with it :D)
-        verts = gtFwdSimComputeCircumscribingPolyhedron(...
-            gr_center, omegas, bb, parameters);
-
-        vol_size = max(verts, [], 1) - min(verts, [], 1);
-        oversize_vol = fsim.oversizeVol * 1.1;
-    else
-        % We should in the future handle properly vertical detector
-        % (general geometry) maybe determining a convex shape of the grain!
-        vol_size = [stackUSize, stackUSize, stackVSize];
-        oversize_vol = fsim.oversizeVol / fsim.oversize;
-    end
-
-    vol_size = round(vol_size * oversize_vol);
-end