FwdSim: added twin variants search

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

zUtil_ForwardSim/gtForwardSimulate_v2.m

@@ -143,6 +143,10 @@ elseif (~isempty(varargin))
     % in interactive mode we can modify parameters
     [parameters.fsim, rej_pars] = parse_pv_pairs(parameters.fsim, varargin);
 end
+twin_params = struct( ...
+    'twin_angles', zeros(0, 1), ...
+    'twin_axes', zeros(0, 3) );
+[twin_params, rej_pars] = parse_pv_pairs(twin_params, rej_pars);
 if (~exist('rej_pars', 'var'))
     rej_pars = {};
 end
@@ -295,6 +299,8 @@ for n = first : last
             parameters, spotsCommProps, fwd_sim.bb(included, :));
     end
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%% We now build the Geometry
     fprintf(' (Done).\n - Building geometry..');
 
     %%% Diffraction geometry
@@ -339,6 +345,30 @@ for n = first : last
         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, :);
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%% We now look for possible twins if asked
+    if (~isempty(twin_params.twin_angles))
+        fprintf(' (Done).\n - Searching twins..');
+
+        g_twins = forward_simulate_twin(gr, fwd_sim, proj, twin_params, ...
+            spotsCommProps, segmentedSpots, parameters, det_index);
+
+        fprintf(' (Done).\n - Num spots per twin:\n')
+        for ii_t = 1:numel(g_twins)
+            fprintf('   %d) %d\n', ii_t, numel(g_twins(ii_t).proj.included));
+        end
+        fprintf(' - Assembling final data-structure..')
+    else
+        fprintf(' (Done).\n - Assembling final data-structure..')
+    end
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %%% Now fill in the information for the output
@@ -371,14 +401,11 @@ for n = first : last
     out.checkSpots = fwd_sim.check_spots;
 
     % diffraction geometry
-    out.proj       = proj;
-    % Selections
-    out.proj.ondet    = onDet;
-    out.proj.included = included;
-    out.proj.selected = selected; % logicals with pre-selected active spots
-    % blob/spot stacks
-    out.proj.bl    = bl(included);
-    out.proj.stack = difstack(:, included, :);
+    out.proj = proj;
+
+    if (~isempty(twin_params.twin_angles))
+        out.g_twins = g_twins;
+    end
 
     found_reflections = numel(find(fwd_sim.flag > 1));
     out.completeness  = found_reflections / numspots;
@@ -616,6 +643,69 @@ function [fwd_sim, bl] = forward_simulate_grain(gr, onDet, ...
     end
 end
 
+function g_twins = forward_simulate_twin(gr, fwd_sim, proj, twin_params, ...
+        spotsCommProps, segmentedSpots, parameters, det_index)
+    spacegroup = parameters.cryst(gr.phaseid).spacegroup;
+    lp = parameters.cryst(gr.phaseid).latticepar;
+
+    g_twins = gtTwinOrientations(gr.R_vector, ...
+        twin_params.twin_angles, twin_params.twin_axes, spacegroup, lp);
+
+    spotsCommProps.XsizeMin = 5;
+    spotsCommProps.YsizeMin = 5;
+    spotsCommProps.max_w_offset_im = 2.5 / parameters.labgeo.omstep;
+    spotsCommProps.max_n_offset_deg = 2.5;
+    spotsCommProps.max_t_offset_deg = 1;
+
+    num_twins = numel(g_twins);
+    fprintf('\b\b, variant: ')
+    for ii = 1:num_twins
+        num_chars = fprintf('%d/%d', ii, num_twins);
+
+        twin = struct('id', gr.id, 'phaseid', gr.phaseid, ...
+            'R_vector', g_twins(ii).R_vector, 'center', gr.center);
+        twin = gtCalculateGrain(twin, parameters);
+
+        twin.stat = gr.stat;
+        twin.stat.bbxsstd = twin.stat.bbxsstd * 0.1;
+        twin.stat.bbysstd = twin.stat.bbysstd * 0.1;
+        twin.difspotidA = [];
+        twin.difspotidB = [];
+
+        uvw = twin.allblobs.detector(det_index).uvw;
+        ondet = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
+
+        proj_bls = predict_spot_masks_from_fwdsimgrain(twin, proj, ondet, parameters);
+
+        [twin.fwd_sim, bl] = forward_simulate_grain(...
+            twin, ondet, spotsCommProps, segmentedSpots, proj_bls, parameters, det_index);
+
+        included = find(twin.fwd_sim.flag >= 3);
+
+        twin.proj(det_index).bl = bl(included);
+        if (~isempty(twin.proj(det_index).bl))
+            spots = arrayfun(@(x){sum(x.intm, 3)}, bl(included));
+            twin.proj(det_index).stack = permute(cat(3, spots{:}), [1 3 2]);
+        end
+
+        twin.proj(det_index).ondet = ondet;
+        twin.proj(det_index).included = included;
+        twin.proj(det_index).selected = false(numel(included), 1);
+
+        g_twins(ii).allblobs = twin.allblobs;
+        g_twins(ii).fwd_sim = clean_fwd_sim(twin.fwd_sim, included);
+        g_twins(ii).proj = twin.proj;
+
+        % Just to make gtGuiGrainMontage happy
+        g_twins(ii).id = gr.id;
+        g_twins(ii).phaseid = gr.phaseid;
+
+        fprintf(repmat('\b', [1 num_chars]))
+    end
+
+    fprintf('%d/%d', num_twins, num_twins)
+end
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Support Functions
 
@@ -705,6 +795,42 @@ 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)
+
+    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)];
+
+    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);
+
+    % 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)];
+
+    proj_bls_original = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);
+
+    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;
+    end
+end
+
 function spot_props = predict_spot_properties(spotsCommProps, proj_bl, gr_stats, parameters)
     bb_size_factor = parameters.fsim.bbsize_factor;
 
@@ -728,6 +854,22 @@ function spot_props = predict_spot_properties(spotsCommProps, proj_bl, gr_stats,
     end
 end
 
+function fwd_sim = clean_fwd_sim(fwd_sim, included)
+%     fwd_sim.flag = fwd_sim.flag(included);
+    fwd_sim.bb = fwd_sim.bb(included, :);
+    fwd_sim.cent = fwd_sim.cent(included);
+%     fwd_sim.spotid = fwd_sim.spotid(included);
+    fwd_sim.difspotID = fwd_sim.spotid(included);
+%     fwd_sim.check = fwd_sim.check(included, :);
+    fwd_sim.intensity = fwd_sim.intensity(included);
+    fwd_sim.avg_pixel_int = fwd_sim.avg_pixel_int(included);
+    fwd_sim.likelihood = fwd_sim.likelihood(included);
+    fwd_sim.cands = fwd_sim.cands(included);
+    fwd_sim.uvw = fwd_sim.uvw(included, :);
+    fwd_sim.spot_type = fwd_sim.spot_type(included);
+    fwd_sim.check_spots = fwd_sim.check_spots(included);
+end
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Likelihood computation