Forward Simulation: fixed twin variant search, and added initial multi-detector search support

(multi-detector not fully implemented yet!)

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_ForwardSim/gtForwardSimulate_v2.m

@@ -118,8 +118,6 @@ end
 labgeo = parameters.labgeo;
 samgeo = parameters.samgeo;
 
-det_index = 1;
-
 if (~isfield(parameters.seg, 'writehdf5') || ~parameters.seg.writehdf5)
     error('gtForwardSimulate_v2:wrong_parameter', ...
         [ 'Segmentation has to write HDF5 blobs in order for ' ...
@@ -143,13 +141,18 @@ elseif (~isempty(varargin))
     % in interactive mode we can modify parameters
     [parameters.fsim, rej_pars] = parse_pv_pairs(parameters.fsim, varargin);
 end
-twin_params = struct( ...
+conf = struct( ...
+    'det_index', [], ...
     'twin_angles', zeros(0, 1), ...
     'twin_axes', zeros(0, 3) );
 if (~exist('rej_pars', 'var'))
     rej_pars = {};
 end
-[twin_params, rej_pars] = parse_pv_pairs(twin_params, rej_pars);
+[conf, rej_pars] = parse_pv_pairs(conf, rej_pars);
+
+if (isempty(conf.det_index))
+    conf.det_index = 1:numel(detgeo);
+end
 
 indexgrainModel = fullfile(parameters.acq.dir, '4_grains', 'phase_%02d', 'index.mat');
 fprintf('\nLoading indexing results...\n');
@@ -218,168 +221,181 @@ for n = first : last
 
     %%% Predict spot positions
     gr = gtCalculateGrain(gr, parameters);
-    uvw = gr.allblobs.detector(det_index).uvw;
 
-    %%% Get:
-    % - average difspot X and Y sizes for the current spots in the grain
-    % - theshold values used during the database search
-    % - offset in possible lateral grain orientation
-    % - diffstack Vertical and Horizontal sizes
-    spotsCommProps = gtFwdSimGetDiffspotsCommonProperties(gr, parameters, det_index);
+    proj = gtFwdSimProjDefinition('fwd_sim', numel(conf.det_index));
+    fwd_sim = repmat(gtFwdSimDataStructureDefinition(0), [numel(conf.det_index), 1]);
 
-    %%% Get the indices of all spots which fully intersect the active
-    % region of the detector
-    onDet = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
+    for ii_d = conf.det_index
+        uvw = gr.allblobs.detector(ii_d).uvw;
 
-    fprintf(' (Done).\n - Producing blob masks..')
+        %%% Get:
+        % - average difspot X and Y sizes for the current spots in the grain
+        % - theshold values used during the database search
+        % - offset in possible lateral grain orientation
+        % - diffstack Vertical and Horizontal sizes
+        spotsCommProps = gtFwdSimGetDiffspotsCommonProperties(gr, parameters, ii_d);
 
-    proj_bls = predict_spot_masks_from_indexed( ...
-        gr, grain, segmentedSpots, spotsCommProps, onDet, parameters, det_index);
-
-    fprintf(' (Done).\n - Preallocating structs..')
-
-    if (parameters.fsim.assemble_figure)
-        full = zeros(detgeo(det_index).detsizev, detgeo(det_index).detsizeu);
-    else
-        full = [];
-    end
-
-    % number of exspected diffraction spots
-    numspots = numel(onDet);
+        %%% Get the indices of all spots which fully intersect the active
+        % region of the detector
+        ondet = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
 
-    % Diffraction vectors directions
-    diff_beam_dirs  = zeros(numspots, 3);
+        fprintf(' (Done).\n - Producing blob masks..')
 
-    fprintf(' (Done).\n - Matching reflections and spots..');
+        proj_bls = predict_spot_masks_from_indexed( ...
+            gr, grain, segmentedSpots, spotsCommProps, ondet, parameters, ii_d);
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %%% Forward simulation:
-    % We try to match the reflections with the segmented spots, so we
-    % examine all hklsp and search corresponding spots / intensity on the
-    % detector..
-    %%% But also build the diffraction blob stack
-    [fwd_sim, bl] = forward_simulate_grain(gr, onDet, spotsCommProps, ...
-        segmentedSpots, proj_bls, parameters, det_index);
+        fprintf(' (Done).\n - Preallocating structs..')
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %%% Re-check flag == 1, to see if there's intensity at the predicted
-    % position, if requested.
-    if (parameters.fsim.check_spot)
-        fprintf(' (Done).\n -  Checking spots in Raw images..');
-
-        toBeChecked = find(fwd_sim.flag == 1);
-        for ii = reshape(toBeChecked, 1, []);
-            [found_intensity, spotInfo] = gtFwdSimCheckSpotInRawImages(...
-                uvw(onDet(ii), :), spotsCommProps, parameters);
-
-            if (found_intensity)
-                % Found intensity in FULL images at predicted position
-                fwd_sim.flag(ii) = 2;
+        if (parameters.fsim.assemble_figure)
+            full = zeros(detgeo(ii_d).detsizev, detgeo(ii_d).detsizeu);
+        else
+            full = [];
+        end
 
-                if (parameters.fsim.assemble_figure)
-                    full = gtPlaceSubImage2(spotInfo.spot, full, spotInfo.bb(1), spotInfo.bb(2), 'summed');
+        % number of exspected diffraction spots
+        numspots = numel(ondet);
+
+        % Diffraction vectors directions
+        diff_beam_dirs  = zeros(numspots, 3);
+
+        fprintf(' (Done).\n - Matching reflections and spots..');
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%% Forward simulation:
+        % We try to match the reflections with the segmented spots, so we
+        % examine all hklsp and search corresponding spots / intensity on
+        % the detector..
+        %%% But also build the diffraction blob stack
+        [fwd_sim(ii_d), bl] = forward_simulate_grain(gr, ondet, ...
+            spotsCommProps, segmentedSpots, proj_bls, parameters, ii_d);
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%% Re-check flag == 1, to see if there's intensity at the
+        % predicted position, if requested.
+        if (parameters.fsim.check_spot)
+            fprintf(' (Done).\n -  Checking spots in Raw images..');
+
+            toBeChecked = find(fwd_sim(ii_d).flag == 1);
+            for ii = reshape(toBeChecked, 1, []);
+                [found_intensity, spotInfo] = gtFwdSimCheckSpotInRawImages(...
+                    uvw(ondet(ii), :), spotsCommProps, parameters);
+
+                if (found_intensity)
+                    % Found intensity in FULL images at predicted position
+                    fwd_sim(ii_d).flag(ii) = 2;
+
+                    if (parameters.fsim.assemble_figure)
+                        full = gtPlaceSubImage2(spotInfo.spot, full, ...
+                            spotInfo.bb(1), spotInfo.bb(2), 'summed');
+                    end
                 end
+                fwd_sim(ii_d).check_spots(ii).info = spotInfo;
             end
-            fwd_sim.check_spots(ii).info = spotInfo;
         end
-    end
-
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    % Reflections that have been matched to segmented spots
-    included = find(fwd_sim.flag >= 3);
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %%% We now build the diffraction stack
-    fprintf(' (Done).\n - Building difstack..');
-
-    difstack = gtFwdSimBuildDifstackSpots(fwd_sim.spotid, included, ...
-        parameters, spotsCommProps.stackUSize, spotsCommProps.stackVSize, ...
-        spotsCommProps.Xsize, spotsCommProps.Ysize, ...
-        fwd_sim.bb, info);
-
-    if (parameters.fsim.assemble_figure)
-        % no padded spots
-        full = sfBuildFull(fwd_sim.spotid(included), ...
-            parameters, spotsCommProps, fwd_sim.bb(included, :));
-    end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % Reflections that have been matched to segmented spots
+        included = find(fwd_sim(ii_d).flag >= 3);
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %%% We now build the Geometry
-    fprintf(' (Done).\n - Building geometry..');
-
-    %%% Diffraction geometry
-    % diff_beam_dirs: direction of diffracted beam used in ASTRA
-    identif_indexed = fwd_sim.flag(included) == 5;
-    incl_indx = included(identif_indexed);
-    incl_not_indx = included(~identif_indexed);
-
-    if (isempty(incl_indx))
-        warning('gtForwardSimulate_v2:wrong_grain', ...
-            ['Grain "%d" is probably a badly indexed grain.' ...
-            ' It will be deselected'], n)
-        sample.phases{phaseID}.setSelected(n, false);
-        continue;
-    end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%% We now build the diffraction stack
+        fprintf(' (Done).\n - Building difstack..');
 
-    % we found a segmented spot - can use the theoretically
-    % predicted direction or apparent direction (based on CoM of
-    % grain and spot) for backprojection
-
-    % If available, we want to use the angles as determined from
-    % the Friedel pairs - interrogate spotpair table
-    [diff_beam_dirs(incl_indx, :), fwd_sim.uvw(incl_indx, 3)] = sfGetFriedelPair( ...
-        pairsTable, fwd_sim.spotid(incl_indx), fwd_sim.spot_type(incl_indx) == 'a');
-
-    % Otherwise we have two choices
-    if (parameters.fsim.use_th)
-        % use theoretically predited directions (default)
-        fwd_sim.uvw(incl_not_indx, 3) = gr.allblobs.omega(onDet(incl_not_indx));
-        % theoretically predicted beam direction (in sample
-        % coordinate system)
-        diff_beam_dirs(incl_not_indx, :) = gr.allblobs.dvecsam(onDet(incl_not_indx), :);
-    else
-        % use experimentally determined omega angles (which we find in: 
-        % fwd_sim.uvw(incl_not_indx, 3))
-        % projection direction as inferred from grain / spot center
-        % positions
-        diff_beam_dirs(incl_not_indx, :) = gtGeoDiffVecInSample( ...
-            fwd_sim.uvw(incl_not_indx, 1:2), fwd_sim.uvw(incl_not_indx, 3), ...
-            gr.center, detgeo(det_index), labgeo, samgeo);
-    end
+        difstack = gtFwdSimBuildDifstackSpots(fwd_sim(ii_d).spotid, ...
+            included, parameters, ...
+            spotsCommProps.stackUSize, spotsCommProps.stackVSize, ...
+            spotsCommProps.Xsize, spotsCommProps.Ysize, ...
+            fwd_sim(ii_d).bb, info);
 
-    %%% Selection of the spots
-    [selected, goodness] = gtFwdSimSelectSpots(fwd_sim.likelihood, fwd_sim.intensity, ...
-        fwd_sim.avg_pixel_int, fwd_sim.flag, included, parameters);
-
-    %%% We now Build ASTRA's geometry
-    proj = gtFwdSimBuildProjGeometry(...
-        diff_beam_dirs(included, :), gr.center, fwd_sim.uvw(included, 3), ...
-        fwd_sim.bb(included, :), parameters, ...
-        spotsCommProps.stackUSize, spotsCommProps.stackVSize, ...
-        selected);
-    % Selections
-    proj.ondet    = onDet;
-    proj.included = included;
-    proj.selected = selected; % logicals with pre-selected active spots
-    % blob/spot stacks
-    proj.bl    = bl(included);
-    proj.stack = difstack(:, included, :);
+        if (parameters.fsim.assemble_figure)
+            % no padded spots
+            full = sfBuildFull(fwd_sim(ii_d).spotid(included), ...
+                parameters, spotsCommProps, fwd_sim(ii_d).bb(included, :));
+        end
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %%% We now look for possible twins if asked
-    if (~isempty(twin_params.twin_angles))
-        fprintf(' (Done).\n - Searching twins..');
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%% We now build the Geometry
+        fprintf(' (Done).\n - Building geometry..');
+
+        %%% Diffraction geometry
+        % diff_beam_dirs: direction of diffracted beam used in ASTRA
+        identif_indexed = fwd_sim(ii_d).flag(included) == 5;
+        incl_indx = included(identif_indexed);
+        incl_not_indx = included(~identif_indexed);
+
+        if (isempty(incl_indx))
+            warning('gtForwardSimulate_v2:wrong_grain', ...
+                ['Grain "%d" is probably a badly indexed grain.' ...
+                ' It will be deselected'], n)
+            sample.phases{phaseID}.setSelected(n, false);
+            continue;
+        end
 
-        fp_twins = forward_simulate_twin(gr, fwd_sim, proj, twin_params, ...
-            spotsCommProps, segmentedSpots, parameters, det_index);
+        % we found a segmented spot - can use the theoretically
+        % predicted direction or apparent direction (based on CoM of
+        % grain and spot) for backprojection
+
+        % If available, we want to use the angles as determined from
+        % the Friedel pairs - interrogate spotpair table
+        [diff_beam_dirs(incl_indx, :), fwd_sim(ii_d).uvw(incl_indx, 3)] = sfGetFriedelPair( ...
+            pairsTable, fwd_sim(ii_d).spotid(incl_indx), fwd_sim(ii_d).spot_type(incl_indx) == 'a');
+
+        % Otherwise we have two choices
+        if (parameters.fsim.use_th)
+            % use theoretically predited directions (default)
+            fwd_sim(ii_d).uvw(incl_not_indx, 3) = gr.allblobs.omega(ondet(incl_not_indx));
+            % theoretically predicted beam direction (in sample
+            % coordinate system)
+            diff_beam_dirs(incl_not_indx, :) = gr.allblobs.dvecsam(ondet(incl_not_indx), :);
+        else
+            % use experimentally determined omega angles (which we find in: 
+            % fwd_sim(ii_d).uvw(incl_not_indx, 3))
+            % projection direction as inferred from grain / spot center
+            % positions
+            diff_beam_dirs(incl_not_indx, :) = gtGeoDiffVecInSample( ...
+                fwd_sim(ii_d).uvw(incl_not_indx, 1:2), ...
+                fwd_sim(ii_d).uvw(incl_not_indx, 3), ...
+                gr.center, detgeo(ii_d), labgeo, samgeo);
+        end
 
-        fprintf(' (Done).\n - Num spots per twin:\n')
-        for ii_t = 1:numel(fp_twins)
-            fprintf('   %d) %d\n', ii_t, numel(fp_twins(ii_t).proj.included));
+        %%% Selection of the spots
+        [selected, goodness] = gtFwdSimSelectSpots( ...
+            fwd_sim(ii_d).likelihood, fwd_sim(ii_d).intensity, ...
+            fwd_sim(ii_d).avg_pixel_int, fwd_sim(ii_d).flag, ...
+            included, parameters);
+
+        %%% We now Build ASTRA's geometry
+        [proj(ii_d), fwd_sim(ii_d).gv_verts] = gtFwdSimBuildProjGeometry( ...
+            diff_beam_dirs(included, :), ...
+            gr.center, fwd_sim(ii_d).uvw(included, 3), ...
+            fwd_sim(ii_d).bb(included, :), parameters, ...
+            spotsCommProps.stackUSize, spotsCommProps.stackVSize, ...
+            selected);
+        % Selections
+        proj(ii_d).ondet    = ondet;
+        proj(ii_d).included = included;
+        % selected: logicals with pre-selected active spots
+        proj(ii_d).selected = selected;
+        % blob/spot stacks
+        proj(ii_d).bl    = bl(included);
+        proj(ii_d).stack = difstack(:, included, :);
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%% We now look for possible twins if asked
+        if (~isempty(conf.twin_angles))
+            fprintf(' (Done).\n - Searching twins..');
+
+            fp_twins = forward_simulate_twin(gr, fwd_sim(ii_d), proj, ...
+                conf, spotsCommProps, segmentedSpots, parameters, ii_d);
+
+            fprintf(' (Done).\n - Num spots per twin:\n')
+            for ii_t = 1:numel(fp_twins)
+                fprintf('   %d) %d\n', ii_t, numel(fp_twins(ii_t).proj.included));
+            end
+            fprintf(' - Assembling final data-structure..')
+        else
+            fprintf(' (Done).\n - Assembling final data-structure..')
         end
-        fprintf(' - Assembling final data-structure..')
-    else
-        fprintf(' (Done).\n - Assembling final data-structure..')
     end
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -391,39 +407,40 @@ for n = first : last
     out.R_vector  = gr.R_vector;
     out.center    = gr.center;
     out.allblobs  = gr.allblobs;
+    out.fwd_sim   = fwd_sim;
 
-    out.check     = fwd_sim.check;
-    out.flag      = fwd_sim.flag;
+    out.check     = fwd_sim(1).check;
+    out.flag      = fwd_sim(1).flag;
     out.pairid    = gr.pairid;
-    out.spotid    = fwd_sim.spotid;
+    out.spotid    = fwd_sim(1).spotid;
 
     % only included spots info
-    out.cands     = fwd_sim.cands(included);
-    out.likelihood = fwd_sim.likelihood(included);
-    out.difspotID = fwd_sim.spotid(included);
-    out.uv_exp    = fwd_sim.uvw(included, 1:2);
-    out.om_exp    = fwd_sim.uvw(included, 3);
-    out.bb_exp    = fwd_sim.bb(included, :);
-    out.intensity = fwd_sim.intensity(included);
+    out.cands     = fwd_sim(1).cands(included);
+    out.likelihood = fwd_sim(1).likelihood(included);
+    out.difspotID = fwd_sim(1).spotid(included);
+    out.uv_exp    = fwd_sim(1).uvw(included, 1:2);
+    out.om_exp    = fwd_sim(1).uvw(included, 3);
+    out.bb_exp    = fwd_sim(1).bb(included, :);
+    out.intensity = fwd_sim(1).intensity(included);
 
     % save full image
     out.full      = full;
 
     % check spot info
-    out.checkSpots = fwd_sim.check_spots;
+    out.checkSpots = fwd_sim(1).check_spots;
 
     % diffraction geometry
     out.proj = proj;
 
-    if (~isempty(twin_params.twin_angles))
+    if (~isempty(conf.twin_angles))
         out.twin_vars = fp_twins;
     else
         out.twin_vars = [];
     end
 
-    found_reflections = numel(find(fwd_sim.flag > 1));
+    found_reflections = numel(find(fwd_sim(1).flag > 1));
     out.completeness  = found_reflections / numspots;
-    out.completeness_incl = numel(find(fwd_sim.flag >= 3)) / numspots;
+    out.completeness_incl = numel(find(fwd_sim(1).flag >= 3)) / numspots;
     out.goodness      = goodness;
     out.goodness_sel  = mean(out.likelihood(selected));
     out.confidence    = out.goodness * out.completeness_incl;
@@ -476,7 +493,7 @@ for n = first : last
 
     fprintf('\n Grain number %d:\n', n);
 
-    newspots = length(find(fwd_sim.flag == 4));
+    newspots = length(find(fwd_sim(1).flag == 4));
     fprintf(['Forward simulation found intensity at %d out of %d expected ' ...
              'spot positions \n'], found_reflections, numspots);
     fprintf('Completeness estimate:  %f\n', out.completeness);
@@ -539,7 +556,8 @@ function fwd_sim = gtFwdSimDataStructureDefinition(numspots)
             'likelihood_dev_w', cell(numspots, 1) ), ...
         'uvw', zeros(numspots, 3), ... % CoM positions of difspots from difspottable
         'spot_type', char(zeros(numspots, 1)), ...
-        'check_spots', struct('info', cell(numspots, 1)) );
+        'check_spots', struct('info', cell(numspots, 1)), ...
+        'gv_verts', [] );
 end
 
 function [fwd_sim, bl] = forward_simulate_grain(gr, onDet, ...
@@ -685,12 +703,12 @@ function g_twins = forward_simulate_twin(gr, fwd_sim, proj, twin_params, ...
         twin.difspotidB = [];
 
         uvw = twin.allblobs.detector(det_index).uvw;
-        ondet = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
+        ondet_twin = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
 
-        proj_bls = predict_spot_masks_from_fwdsimgrain(twin, proj, ondet, parameters, det_index);
+        proj_bls = predict_spot_masks_from_fwdsimgrain(twin, gr, proj, fwd_sim, ondet_twin, parameters, det_index);
 
         [twin.fwd_sim, bl] = forward_simulate_grain(...
-            twin, ondet, spotsCommProps, segmentedSpots, proj_bls, parameters, det_index);
+            twin, ondet_twin, spotsCommProps, segmentedSpots, proj_bls, parameters, det_index);
 
         included = find(twin.fwd_sim.flag >= 3);
 
@@ -700,7 +718,7 @@ function g_twins = forward_simulate_twin(gr, fwd_sim, proj, twin_params, ...
             twin.proj(det_index).stack = permute(cat(3, spots{:}), [1 3 2]);
         end
 
-        twin.proj(det_index).ondet = ondet;
+        twin.proj(det_index).ondet = ondet_twin;
         twin.proj(det_index).included = included;
         twin.proj(det_index).selected = false(numel(included), 1);
 
@@ -770,7 +788,7 @@ function displayShowFsim(phaseID, grain, varargin)
 end
 
 function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
-    segmentedSpots, spotsCommProps, onDet, parameters, det_ind)
+    segmentedSpots, spotsCommProps, ondet, parameters, det_ind)
 
     samgeo = parameters.samgeo;
     labgeo = parameters.labgeo;
@@ -798,7 +816,9 @@ function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
     projvecs = gtGeoDet2Lab(projvecs, detgeo, false);
     projvecs = gtGeoLab2Sam(projvecs, Ws, labgeo, samgeo, 0, 0) - gr_center;
 
-    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
+    verts_rec = gtFwdSimComputeCircumscribingPolyhedron(gr.center, projvecs, Ws, BBs, parameters, det_ind, 'convhull');
+
+    proj_bls = gtFwdSimProjectVertices2Det(gr, verts_rec, ondet, parameters, det_ind);
 
     % Producing original blobs' boundingboxes (as by original segmentation
     % info)
@@ -808,7 +828,9 @@ function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
         segmentedSpots.BoundingBoxXsize(difspotID_INDXTR), ...
         segmentedSpots.BoundingBoxYsize(difspotID_INDXTR) ];
 
-    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
+    verts_rec = gtFwdSimComputeCircumscribingPolyhedron(gr.center, projvecs, Ws, BBs, parameters, det_ind, 'convhull');
+
+    proj_bls_original = gtFwdSimProjectVertices2Det(gr, verts_rec, ondet, parameters, det_ind);
 
     for ii_b = 1:numel(proj_bls)
         proj_bls(ii_b).mbbsize = proj_bls_original(ii_b).bbsize;
@@ -820,35 +842,48 @@ function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
     end
 end
 
-function proj_bls = predict_spot_masks_from_fwdsimgrain(gr, proj, onDet, parameters, det_ind)
+function proj_bls = predict_spot_masks_from_fwdsimgrain(new_gr, main_gr, proj, main_gr_fwd_sim, new_gr_ondet, parameters, det_ind)
+% We here produce the expected bounding box using the main grain's proj
+% datastructure, and then we produce the projected blobs at the predicted
+% new grain's positions
+
+    if (isempty(main_gr_fwd_sim.gv_verts))
+        centroids_w = zeros(numel(proj.bl), 1);
+        for ii_b = 1:numel(proj.bl)
+            int_profile = reshape(sum(sum(proj.bl(ii_b).intm, 1), 2), 1, []);
+            num_slices = numel(int_profile);
+            centroids_w(ii_b) = sum((0:num_slices-1) .* int_profile) / num_slices;
+        end
+        Ws = cat(1, proj.bl(:).bbwim);
+        Ws = (Ws(:, 1) + centroids_w) .* gtGetOmegaStepDeg(parameters, det_ind);
 
-    centroids = zeros(numel(proj.bl), 1);
-    for ii_b = 1:numel(proj.bl)
-        int_profile = reshape(sum(sum(proj.bl(ii_b).intm, 1), 2), 1, []);
-        num_slices = numel(int_profile);
-        centroids(ii_b) = sum((0:num_slices-1) .* int_profile) / num_slices;
-    end
-    Ws = cat(1, proj.bl(:).bbwim);
-    Ws = Ws(:, 1) + centroids;
+        % Producing image blobs' boundingboxes (as by blob size info)
+        BBus = cat(1, proj.bl(:).bbuim);
+        BBvs = cat(1, proj.bl(:).bbvim);
+        BBsizes = cat(1, proj.bl(:).bbsize);
+        BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
 
-    % Producing image blobs' boundingboxes (as by blob size info)
-    BBus = cat(1, proj.bl(:).bbuim);
-    BBvs = cat(1, proj.bl(:).bbvim);
-    BBsizes = cat(1, proj.bl(:).bbsize);
-    BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
+        projvecs = main_gr.allblobs.dvecsam(proj.ondet(proj.included), :);
 
-    projvecs = gr.allblobs.dvecsam(onDet, :);
+        verts_rec = gtFwdSimComputeCircumscribingPolyhedron(main_gr.center, projvecs, Ws, BBs, parameters, det_ind, 'convhull');
+    else
+        verts_rec = main_gr_fwd_sim.gv_verts;
+    end
 
-    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
+    proj_bls = gtFwdSimProjectVertices2Det(new_gr, verts_rec, new_gr_ondet, parameters, det_ind);
 
-    % Producing original blobs' boundingboxes (as by original segmentation
-    % info)
-    BBus = cat(1, proj.bl(:).mbbu);
-    BBvs = cat(1, proj.bl(:).mbbv);
-    BBsizes = cat(1, proj.bl(:).mbbsize);
-    BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
+    if (isempty(main_gr_fwd_sim.gv_verts))
+        % Producing original blobs' boundingboxes (as by original segmentation
+        % info)
+        BBus = cat(1, proj.bl(:).mbbu);
+        BBvs = cat(1, proj.bl(:).mbbv);
+        BBsizes = cat(1, proj.bl(:).mbbsize);
+        BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];
+
+        verts_rec = gtFwdSimComputeCircumscribingPolyhedron(main_gr.center, projvecs, Ws, BBs, parameters, det_ind, 'convhull');
+    end
 
-    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, projvecs, BBs, Ws, onDet, parameters, det_ind);
+    proj_bls_original = gtFwdSimProjectVertices2Det(new_gr, verts_rec, new_gr_ondet, parameters, det_ind);
 
     for ii_b = 1:numel(proj_bls)
         proj_bls(ii_b).mbbsize = proj_bls_original(ii_b).bbsize;

zUtil_ForwardSim/gtFwdSimBuildProjGeometry.m

-function proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb, parameters, stackUSize, stackVSize, selected, det_ind)
+function [proj, verts] = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb, parameters, stackUSize, stackVSize, selected, det_ind)
 % FUNCTION proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb, parameters, stackUSize, stackVSize, selected, det_ind)
 %
     if (~exist('det_ind', 'var') || isempty(det_ind))
@@ -47,6 +47,8 @@ function proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb,
         bbpos_det_abs, [], detgeo, labgeo, samgeo, recgeo,  ...
         'ASTRA_absorption');
 
+    proj = gtFwdSimProjDefinition('fwd_sim');
+
     % diffraction geometry
     proj.geom      = proj_geom;
     proj.geom_full = proj_geom_full;
@@ -67,6 +69,7 @@ function proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb,
 
         vol_size = 2 * max(abs(max(verts, [], 1)), abs(min(verts, [], 1)));
     else
+        verts = [];
         % 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;

zUtil_ForwardSim/gtFwdSimComputeCircumscribingPolyhedron.m

@@ -31,7 +31,7 @@ function verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, projvecs, om
 
     % All functions that use the same omegas, need the same rotation tensors:
     rot_comp_w = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
-    rot_omega = gtMathsRotationTensor(omegas, rot_comp_w);
+    rot_s2l_w = gtMathsRotationTensor(omegas, rot_comp_w);
 
     gc_pos_rec = gtGeoSam2Sam(gr_center, samgeo, recgeo, true);
     gc_pos_rec = gc_pos_rec(ones(numel(omegas), 1), :);
@@ -45,10 +45,10 @@ function verts = gtFwdSimComputeCircumscribingPolyhedron(gr_center, projvecs, om
     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);
+    ctl_pos_rec = gtGeoLab2Sam(ctl_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
+    ctr_pos_rec = gtGeoLab2Sam(ctr_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
+    cbl_pos_rec = gtGeoLab2Sam(cbl_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
+    cbr_pos_rec = gtGeoLab2Sam(cbr_pos_lab, rot_s2l_w, labgeo, recgeo, 0, 0);
 
     projvecs_rec = gtGeoSam2Sam(projvecs, samgeo, recgeo, true);
     projvecs_rec = gtMathsNormalizeVectorsList(projvecs_rec);

zUtil_ForwardSim/gtFwdSimPredictProjectedGrainBBox.m

@@ -2,78 +2,8 @@ function proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, sel_projvecs, sel_bb,
     if (~exist('det_ind', 'var') || isempty(det_ind))
         det_ind = 1;
     end
-    detgeo = parameters.detgeo(det_ind);
-    labgeo = parameters.labgeo;
-    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_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, :);
-    omegas = gr.allblobs.omega(ondet);
-%     uv = gr.allblobs.detector(det_ind).uvw(ondet, 1:2);
-
-    num_omegas = numel(omegas);
-    proj_bls = struct( ...
-        'mask', cell(1, num_omegas), ...
-        'bbuim', cell(1, num_omegas), ...
-        'bbvim', cell(1, num_omegas), ...
-        'bbwim', cell(1, num_omegas), ...
-        'bbsize', cell(1, num_omegas) ...
-        );
-
-    % All functions that use the same omegas, need the same rotation tensors:
-    % Notice the sign of W !! This is going to be used for a Sam -> Lab
-    rot_comp_w = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
-    rot_omega = gtMathsRotationTensor(-omegas, rot_comp_w);
-
-    ones_verts = ones(size(verts_rec, 1), 1);
-    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);
-        verts_lab = gtGeoSam2Lab(verts_sam, rotw, labgeo, samgeo, false);
-        dvec_vers = dvecs_lab(ii .* ones_verts, :);
-
-        verts_det_pos = gtFedPredictUVMultiple([], dvec_vers', verts_lab', ...
-            detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
-            [detgeo.detrefu, detgeo.detrefv]')';
-
-%         f = figure();
-%         ax = axes('parent', f);
-%         hold(ax, 'on')
-%         scatter(ax, verts_det_pos(:, 1), verts_det_pos(:, 2));
-%         scatter(ax, uv(ii, 1), uv(ii, 2), 30, 'r');
-%         hold(ax, 'off')
-
-        verts_idxs = convhull(verts_det_pos);
-        verts_det_pos = verts_det_pos(verts_idxs, :);
-
-        bb_2d = [floor(min(verts_det_pos, [], 1)), ceil(max(verts_det_pos, [], 1))];
-        bb_2d = [bb_2d(1:2), (bb_2d(3:4) - bb_2d(1:2) + 1)];
-
-        mask = zeros(bb_2d(3:4));
-        verts_mask_pos = round(verts_det_pos) - bb_2d(ones(numel(verts_idxs), 1), 1:2) + 1;
-        verts_mask_indx = verts_mask_pos(:, 1) + bb_2d(3) * (verts_mask_pos(:, 2) - 1);
-        mask(verts_mask_indx) = 1;
-        mask = bwconvhull(mask);
-
-%         f = figure();
-%         ax = axes('parent', f);
-%         imshow(mask, 'parent', ax)
-
-        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)] / omstep);
-        proj_bls(ii).bbsize = [bb_2d(3:4), 1];
-        proj_bls(ii).mask = mask;
-
-%         pause
-    end
+    proj_bls = gtFwdSimProjectVertices2Det(gr, verts_rec, ondet, parameters, det_ind);
 end
\ No newline at end of file

zUtil_ForwardSim/gtFwdSimProjDefinition.m

+function proj = gtFwdSimProjDefinition(type, num_projs)
+    proj = struct( ...
+        'stack', [], ...
+        'num_rows', 0, ...
+        'num_cols', 0, ...
+        'bl', [], ...
+        'vol_size_x', 0, ...
+        'vol_size_y', 0, ...
+        'vol_size_z', 0, ...
+        'centerpix', zeros(1, 3), ...
+        'shift', zeros(1, 3), ...
+        'ondet', zeros(0, 1), ...
+        'included', zeros(0, 1), ...
+        'selected', zeros(0, 1) );
+
+    if (exist('type', 'var') && ~isempty(type) && strcmpi(type, 'fwd_sim'))
+        proj.geom = zeros(0, 12);
+        proj.geom_full = zeros(0, 12);
+        proj.geom_abs = zeros(0, 12);
+    end
+
+    if (exist('num_projs', 'var') && ~isempty(num_projs))
+        proj(2:num_projs) = proj;
+    end
+end

zUtil_ForwardSim/gtFwdSimProjectVertices2Det.m

+ function proj_bls = gtFwdSimProjectVertices2Det(gr, verts_rec, ondet, parameters, det_ind, verbose)
+% FUNCTION proj_bls = gtFwdSimProjectVertices2Det(gr, verts_rec, ondet, parameters, det_ind, verbose)
+%
+
+    if (~exist('verbose', 'var') || isempty(verbose))
+       verbose = false;
+    end
+
+    detgeo = parameters.detgeo(det_ind);
+    labgeo = parameters.labgeo;
+    samgeo = parameters.samgeo;
+    recgeo = parameters.recgeo(det_ind);
+
+    omstep = gtGetOmegaStepDeg(parameters, det_ind);
+
+    % Now we project the vertices to the spots and take the convex hull
+    dvecs_lab = gr.allblobs.dvec(ondet, :);
+    omegas = gr.allblobs.omega(ondet);
+    if (verbose)
+        uv = gr.allblobs.detector(det_ind).uvw(ondet, 1:2);
+    end
+
+    num_omegas = numel(omegas);
+    proj_bls = gtFwdSimBlobDefinition('blob', num_omegas);
+
+    % All functions that use the same omegas, need the same rotation tensors:
+    % Notice the sign of W !! This is going to be used for a Sam -> Lab
+    rot_l2s = gr.allblobs.srot(:, :, ondet);
+    rot_s2l = permute(rot_l2s, [2 1 3]);
+
+    ones_verts = ones(size(verts_rec, 1), 1);
+    gc_sam = gr.center(ones_verts, :);
+
+    verts_sam = gtGeoSam2Sam(verts_rec, recgeo, samgeo, false) + gc_sam;
+
+    for ii = 1:num_omegas
+        rot_s2l_w = rot_s2l(:, :, ii);
+        verts_lab = gtGeoSam2Lab(verts_sam, rot_s2l_w, labgeo, samgeo, false);
+        dvec_vers = dvecs_lab(ii .* ones_verts, :);
+
+        verts_det_pos = gtFedPredictUVMultiple([], dvec_vers', verts_lab', ...
+            detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
+            [detgeo.detrefu, detgeo.detrefv]')';
+
+        if (verbose)
+            f = figure();
+            ax = axes('parent', f);
+            hold(ax, 'on')
+            scatter(ax, verts_det_pos(:, 1), verts_det_pos(:, 2));
+            scatter(ax, uv(ii, 1), uv(ii, 2), 30, 'r');
+            hold(ax, 'off')
+        end
+
+        verts_idxs = convhull(verts_det_pos);
+        verts_det_pos = verts_det_pos(verts_idxs, :);
+
+        bb_2d = [floor(min(verts_det_pos, [], 1)), ceil(max(verts_det_pos, [], 1))];
+        bb_2d = [bb_2d(1:2), (bb_2d(3:4) - bb_2d(1:2) + 1)];
+
+        mask = zeros(bb_2d(3:4));
+        verts_mask_pos = round(verts_det_pos) - bb_2d(ones(numel(verts_idxs), 1), 1:2) + 1;
+        verts_mask_indx = verts_mask_pos(:, 1) + bb_2d(3) * (verts_mask_pos(:, 2) - 1);
+        mask(verts_mask_indx) = 1;
+        mask = bwconvhull(mask);
+
+        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)] / omstep);
+        proj_bls(ii).bbsize = [bb_2d(3:4), 1];
+        proj_bls(ii).mask = mask;
+
+        if (verbose)
+            f = figure();
+            ax = axes('parent', f);
+            imshow(mask, 'parent', ax)
+
+            pause
+        end
+    end
+end
\ No newline at end of file