FwdSim/Candidates-discovery: saving all the criteria for acceptance, to be able to debug it

And added multi-detector enabling as well

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

zUtil_ForwardSim/gtForwardSimulate_v2.m

@@ -143,6 +143,7 @@ elseif (~isempty(varargin))
 end
 conf = struct( ...
     'det_index', [], ...
+    'debug', false, ...
     'twin_angles', zeros(0, 1), ...
     'twin_axes', zeros(0, 3), ...
     'twin_angular_toll_deg', 0.5 ... Could be raised up to 2.5 or even 5
@@ -281,7 +282,7 @@ for n = first : last
         % 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);
+            spotsCommProps, segmentedSpots, proj_bls, parameters, conf, ii_d);
 
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         %%% Re-check flag == 1, to see if there's intensity at the
@@ -559,7 +560,7 @@ end
 %%% Forward Simulation
 
 function [fwd_sim, bl] = forward_simulate_grain(gr, ondet, ...
-        spotsCommProps, segmentedSpots, proj_bls, parameters, det_index)
+        spotsCommProps, segmentedSpots, proj_bls, parameters, conf, det_index)
 
     uvw = gr.allblobs.detector(det_index).uvw;
 
@@ -581,13 +582,17 @@ function [fwd_sim, bl] = forward_simulate_grain(gr, ondet, ...
 
         % Look for difspots appying loose search criteria - is there any intensity at the predicted spot
         % position ?  Choose the one which is closest to the expected size
-        [candidateIDs, spotsCands] = gtFwdSimFindCandidates(spot_props, gr, ondet(ii), segmentedSpots, parameters, det_index);
+        [candidateIDs, cands_info] = gtFwdSimFindCandidates(spot_props, gr, ondet(ii), segmentedSpots, parameters, det_index);
+
+        if (conf.debug)
+            fwd_sim.cands_info(ii).info = cands_info;
+        end
 
         if (isempty(candidateIDs))
             % NO CANDIDATES FOUND:
             %   segmentation / overlap / intensity / ... problem ?
 
-            if (spotsCands.found_intensity)
+            if (cands_info.found_intensity)
                 % A bigger spot was found in the region.. it may be an
                 % overlap
                 fwd_sim.flag(ii) = 2;
@@ -609,12 +614,12 @@ function [fwd_sim, bl] = forward_simulate_grain(gr, ondet, ...
 
         % Computing UV deviation/likelihood
         th_uv = uvw(refl_indx, 1:2);
-        blobs_uv = [ spotsCands.CentroidX, spotsCands.CentroidY ];
+        blobs_uv = [ cands_info.CentroidX, cands_info.CentroidY ];
         ls_uv = compute_uv_partial_likelihoods(blobs_uv, th_uv, proj_bl);
 
         % Computing W deviation/likelihood
         th_w = gr.allblobs.omega(refl_indx);
-        blob_w = spotsCands.CentroidImage * parameters.labgeo.omstep;
+        blob_w = cands_info.CentroidImage * parameters.labgeo.omstep;
         ls_w = compute_w_partial_likelihoods(blob_w, th_w);
 
         % in case there are several candidates, take the one with best
@@ -633,14 +638,14 @@ function [fwd_sim, bl] = forward_simulate_grain(gr, ondet, ...
         bl(ii) = bls(cand_inc_index);
 
         fwd_sim.uvw(ii, :) = [ ...
-            spotsCands.CentroidX(cand_inc_index), ...
-            spotsCands.CentroidY(cand_inc_index), ...
-            spotsCands.CentroidImage(cand_inc_index) * parameters.labgeo.omstep ];
+            cands_info.CentroidX(cand_inc_index), ...
+            cands_info.CentroidY(cand_inc_index), ...
+            cands_info.CentroidImage(cand_inc_index) * parameters.labgeo.omstep ];
         fwd_sim.bb(ii, :) = [...
-            spotsCands.BoundingBoxXorigin(cand_inc_index), ...
-            spotsCands.BoundingBoxYorigin(cand_inc_index), ...
-            spotsCands.BoundingBoxXsize(cand_inc_index), ...
-            spotsCands.BoundingBoxYsize(cand_inc_index) ];
+            cands_info.BoundingBoxXorigin(cand_inc_index), ...
+            cands_info.BoundingBoxYorigin(cand_inc_index), ...
+            cands_info.BoundingBoxXsize(cand_inc_index), ...
+            cands_info.BoundingBoxYsize(cand_inc_index) ];
         % The following should be beter than the previous, but it seems to
         % geneterate problems
 %         fwd_sim.bb(ii, :) = bls(cand_inc_index).mbbsize([1 2 4 5]);
@@ -718,7 +723,7 @@ function twin_vars = forward_simulate_twin(gr, gr_fwd_sim, gr_proj, conf, spotsC
                 gr_proj(ii_d), gr_fwd_sim(ii_d), twin_ondet, parameters, ii_d);
 
             [twin_fwd_sim(ii_d), bl] = forward_simulate_grain(...
-                twin, twin_ondet, spotsCommProps, segmentedSpots, proj_bls, parameters, ii_d);
+                twin, twin_ondet, spotsCommProps, segmentedSpots, proj_bls, parameters, conf, ii_d);
 
             twin_incl = find(twin_fwd_sim(ii_d).flag >= 3);
             twin_fwd_sim(ii_d) = clean_fwd_sim(twin_fwd_sim(ii_d), twin_incl);
@@ -917,8 +922,14 @@ function spot_props = predict_spot_properties(spotsCommProps, proj_bl, gr_stats,
     bb_size_factor = parameters.fsim.bbsize_factor;
 
     % Minimums could be wrong by up to 2 pixels (roundings)
-    spot_props.XsizeMin = max(5, proj_bl.mbbsize(1) - bb_size_factor * gr_stats.bbxsstd - 2);
-    spot_props.YsizeMin = max(5, proj_bl.mbbsize(2) - bb_size_factor * gr_stats.bbysstd - 2);
+    % Recentl we introduced a grain volume bounding box estimation that
+    % probably considerably over-estimates the grain volume (moving from
+    % inscribed polyhedron to convex-hull), so the new minimum sizes have
+    % to be adusted
+%     spot_props.XsizeMin = max(5, proj_bl.mbbsize(1) - bb_size_factor * gr_stats.bbxsstd - 2)
+%     spot_props.YsizeMin = max(5, proj_bl.mbbsize(2) - bb_size_factor * gr_stats.bbysstd - 2)
+    spot_props.XsizeMin = max(5, proj_bl.mbbsize(1) / bb_size_factor - bb_size_factor * gr_stats.bbxsstd - 2);
+    spot_props.YsizeMin = max(5, proj_bl.mbbsize(2) / bb_size_factor - bb_size_factor * gr_stats.bbysstd - 2);
 
     % Allowing for matching spots, even if the convex shape of the grain
     % computed from indexter's selection had a limited angle view on the
@@ -941,6 +952,7 @@ function fwd_sim = clean_fwd_sim(fwd_sim, included)
     % fwd_sim.flag = fwd_sim.flag(included);
     % fwd_sim.spotid = fwd_sim.spotid(included);
     % fwd_sim.check = fwd_sim.check(included, :);
+    % fwd_sim.check_spots = fwd_sim.check_spots(included);
 
     % And the traditionally included:
     fwd_sim.bb = fwd_sim.bb(included, :);
@@ -951,7 +963,6 @@ function fwd_sim = clean_fwd_sim(fwd_sim, 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);
 
     % And a new field:
     fwd_sim.difspotID = fwd_sim.spotid(included);

zUtil_ForwardSim/gtFwdSimDataStructureDefinition.m

@@ -15,6 +15,7 @@ function fwd_sim = gtFwdSimDataStructureDefinition(numspots)
             'likelihood_proj_mask', cell(numspots, 1), ...
             'likelihood_dev_uv', cell(numspots, 1), ...
             'likelihood_dev_w', cell(numspots, 1) ), ...
+        'cands_info', struct('info', 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)), ...

zUtil_ForwardSim/gtFwdSimFindCandidates.m

-function [candidateIDs, spotsCands] = gtFwdSimFindCandidates(spot_props, gr, spot_index, segmentedSpots, parameters, det_index)
+function [candidateIDs, cands_info] = gtFwdSimFindCandidates(spot_props, gr, spot_index, segmentedSpots, parameters, det_index)
 % FINDCANDIDATES
-% [candidateIDs, spotsCands] = gtFwdSimFindCandidates(spotsCommProps, spotsPositions, segmentedSpots)
+% [candidateIDs, spotsCands] = gtFwdSimFindCandidates(spot_props, gr, spot_index, segmentedSpots, parameters, det_index)
 %
 % considers image number, then position, then area
 % uvw : spot position
@@ -14,85 +14,86 @@ function [candidateIDs, spotsCands] = gtFwdSimFindCandidates(spot_props, gr, spo
         det_index = 1;
     end
 
+    spots_info = segmentedSpots(det_index);
+
     uvw = gr.allblobs.detector(det_index).uvw(spot_index, :);
 
-    image_ok = find(...
-        (segmentedSpots.CentroidImage >= (uvw(3) - spot_props.max_w_offset_im) ) ...
-        & (segmentedSpots.CentroidImage <= (uvw(3) + spot_props.max_w_offset_im) ));
+    % centroid distances
+    dist_u_pix = abs(spots_info.BoundingBoxXorigin + spots_info.BoundingBoxXsize / 2 - uvw(1));
+    dist_v_pix = abs(spots_info.BoundingBoxYorigin + spots_info.BoundingBoxYsize / 2 - uvw(2));
+    dist_w_im = abs(spots_info.CentroidImage - uvw(3));
 
-    % centroid distance
-    Xdist = abs(segmentedSpots.BoundingBoxXorigin(image_ok) ...
-        + 0.5 * segmentedSpots.BoundingBoxXsize(image_ok) - uvw(1));
-    Ydist = abs(segmentedSpots.BoundingBoxYorigin(image_ok) ...
-        + 0.5 * segmentedSpots.BoundingBoxYsize(image_ok) - uvw(2));
+%     dist_u_ok = dist_u_pix < (spots_info.BoundingBoxXsize - spot_props.XsizeMin / 2) / 2;
+%     dist_v_ok = dist_v_pix < (spots_info.BoundingBoxYsize - spot_props.YsizeMin / 2) / 2;
+    dist_u_ok = dist_u_pix < (spot_props.XsizeMax / 2);
+    dist_v_ok = dist_v_pix < (spot_props.YsizeMax / 2);
+    dist_w_ok = dist_w_im <=  spot_props.max_w_offset_im;
 
     if (isfield(spot_props, 'max_n_offset') && ~isempty(spot_props.max_n_offset))
         labgeo = parameters.labgeo;
         samgeo = parameters.samgeo;
-        detgeo = parameters.detgeo(det_index);
 
-        spots_w = segmentedSpots.CentroidImage(image_ok) * labgeo.omstep;
+        omstep = gtGetOmegaStepDeg(parameters, det_index);
+
+        spots_w = spots_info.CentroidImage * omstep;
         gr_center = gr.center(ones(numel(spots_w), 1), :);
         gr_pos = gtGeoSam2Lab(gr_center, spots_w, labgeo, samgeo, false);
 
-        spots_uv = [ ...
-            segmentedSpots.BoundingBoxXorigin(image_ok), ...
-            segmentedSpots.BoundingBoxYorigin(image_ok) ];
-        spots_diff_vec = gtGeoDet2Lab(spots_uv, detgeo, false);
-        spots_diff_vec = spots_diff_vec - gr_pos;
-        spots_diff_vec = gtMathsNormalizeVectorsList(spots_diff_vec);
-        spots_eta = gtGeoEtaFromDiffVec(spots_diff_vec, labgeo);
-
-        n = gr.allblobs.eta(spot_index);
-        eta_ok = abs(spots_eta - n) < spot_props.max_n_offset_deg;
+        spots_diff_vec = gtMathsNormalizeVectorsList(spots_info.ScatteringVector - gr_pos);
 
+        spots_eta = gtGeoEtaFromDiffVec(spots_diff_vec, labgeo);
         spots_theta = gtGeoThetaFromDiffVec(spots_diff_vec, labgeo);
 
+        n = gr.allblobs.eta(spot_index);
         t = gr.allblobs.theta(spot_index);
-        theta_ok = abs(spots_theta - t) < spot_props.max_t_offset_deg;
 
-        Xdist_ok = (spot_props.XsizeMax * 2) > Xdist;
-        Ydist_ok = (spot_props.YsizeMin * 2) > Ydist;
-%         position_ok = Xdist_ok & Ydist_ok & eta_ok & theta_ok;
-        position_ok = eta_ok & theta_ok;
+        dist_eta_ok = abs(spots_eta - n) < spot_props.max_n_offset_deg;
+        dist_theta_ok = abs(spots_theta - t) < spot_props.max_t_offset_deg;
+
+        position_ok = dist_w_ok & dist_eta_ok & dist_theta_ok;
     else
-        Xdist_ok = segmentedSpots.BoundingBoxXsize(image_ok) > 2*(Xdist + 0.5 * spot_props.XsizeMin / 2);
-        Ydist_ok = segmentedSpots.BoundingBoxYsize(image_ok) > 2*(Ydist + 0.5 * spot_props.YsizeMin / 2);
-        position_ok = Xdist_ok & Ydist_ok;
+        position_ok = dist_w_ok & dist_u_ok & dist_v_ok;
     end
 
-    image_pos_ok = image_ok(position_ok);
-
     % candidates must have a minimum and maximum size
-    size_x_min_ok = segmentedSpots.BoundingBoxXsize(image_pos_ok) >= spot_props.XsizeMin;
-    size_x_max_ok = segmentedSpots.BoundingBoxXsize(image_pos_ok) <= spot_props.XsizeMax;
-    size_y_min_ok = segmentedSpots.BoundingBoxYsize(image_pos_ok) >= spot_props.YsizeMin;
-    size_y_max_ok = segmentedSpots.BoundingBoxYsize(image_pos_ok) <= spot_props.YsizeMax;
+    size_x_min_ok = spots_info.BoundingBoxXsize >= spot_props.XsizeMin;
+    size_x_max_ok = spots_info.BoundingBoxXsize <= spot_props.XsizeMax;
+    size_y_min_ok = spots_info.BoundingBoxYsize >= spot_props.YsizeMin;
+    size_y_max_ok = spots_info.BoundingBoxYsize <= spot_props.YsizeMax;
 
     area_min_ok = size_x_min_ok & size_y_min_ok;
     area_max_ok = size_x_max_ok & size_y_max_ok;
-    image_area_ok = image_pos_ok(area_min_ok & area_max_ok);
-
-    candidateIDs = image_area_ok;
-
-    spotsCands = sfGetTableIndexes(segmentedSpots, candidateIDs);
-    spotsCands.thr.area_thr  = 0.5 * spot_props.XsizeMin .* spot_props.YsizeMin;
-    spotsCands.thr.image_thr = spot_props.max_w_offset_im;
-    spotsCands.thr.xdist_thr = (0.5 * spot_props.XsizeMin / 2)*2;
-    spotsCands.thr.ydist_thr = (0.5 * spot_props.YsizeMin / 2)*2;
-    spotsCands.diff.Xdist         = Xdist;
-    spotsCands.diff.Ydist         = Ydist;
-    spotsCands.diff.image_ok      = image_ok;
-    spotsCands.diff.Xdist_ok      = Xdist_ok;
-    spotsCands.diff.Ydist_ok      = Ydist_ok;
-    spotsCands.diff.position_ok   = position_ok;
-    spotsCands.diff.image_pos_ok  = image_pos_ok ;
-    spotsCands.diff.area_ok       = area_min_ok & area_max_ok;
-    spotsCands.diff.image_area_ok = image_area_ok;
+    image_area_ok = area_min_ok & area_max_ok;
+
+    candidateIDs = find(position_ok & image_area_ok);
+
+    cands_info = sfGetTableIndexes(segmentedSpots, candidateIDs);
+    cands_info.thr.area_min_thr = spot_props.XsizeMin .* spot_props.YsizeMin;
+    cands_info.thr.area_max_thr = spot_props.XsizeMax .* spot_props.YsizeMax;
+    cands_info.thr.w_dist_im_thr = spot_props.max_w_offset_im;
+    cands_info.thr.u_dist_pix_thr = spot_props.XsizeMax / 2;
+    cands_info.thr.v_dist_pix_thr = spot_props.YsizeMax / 2;
+    cands_info.diff.dist_u_pix    = dist_u_pix;
+    cands_info.diff.dist_v_pix    = dist_v_pix;
+    cands_info.diff.dist_w_im     = dist_w_im;
+    cands_info.diff.dist_u_ok     = dist_u_ok;
+    cands_info.diff.dist_v_ok     = dist_v_ok;
+    cands_info.diff.dist_w_ok     = dist_w_ok;
+    if (isfield(spot_props, 'max_n_offset') && ~isempty(spot_props.max_n_offset))
+        cands_info.diff.dist_n_ok = dist_eta_ok;
+        cands_info.diff.dist_t_ok = dist_theta_ok;
+    else
+        cands_info.diff.dist_n_ok = [];
+        cands_info.diff.dist_t_ok = [];
+    end
+    cands_info.diff.position_ok   = position_ok;
+    cands_info.diff.area_min_ok   = area_min_ok;
+    cands_info.diff.area_max_ok   = area_max_ok;
+    cands_info.diff.area_ok       = area_min_ok & area_max_ok;
 
     % This is for completeness estimation: even if the candidate is not
     % good, but if it has a minimum size, it might just be an overlap
-    spotsCands.found_intensity = any(area_min_ok);
+    cands_info.found_intensity = any(position_ok(area_min_ok));
 end
 
 function spotsCands = sfGetTableIndexes(segmentedSpots, candidateIDs)