FwdSim: tweaked spot selection, after the recent changes in candidate inclusion

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

zUtil_ForwardSim/gtForwardSimulate_v2.m

@@ -156,6 +156,7 @@ end
 if (isempty(conf.det_index))
     conf.det_index = 1:numel(detgeo);
 end
+num_dets = numel(conf.det_index);
 
 indexgrainModel = fullfile(parameters.acq.dir, '4_grains', 'phase_%02d', 'index.mat');
 fprintf('\nLoading indexing results...\n');
@@ -182,8 +183,8 @@ end
 % !!! using a mean pixel size is uncertain; ideally fwd simulation should
 % use mm-s
 
-segmentedSpots = cell(numel(conf.det_index), 1);
-pairsTable = cell(numel(conf.det_index), 1);
+segmentedSpots = cell(num_dets, 1);
+pairsTable = cell(num_dets, 1);
 
 fprintf('\nLoading reflections from database..');
 for ii_d = conf.det_index
@@ -240,8 +241,8 @@ for n = first : last
     %%% Predict spot positions
     gr = gtCalculateGrain(gr, parameters);
 
-    proj = gtFwdSimProjDefinition('fwd_sim', numel(conf.det_index));
-    fwd_sim = repmat(gtFwdSimDataStructureDefinition(0), [numel(conf.det_index), 1]);
+    proj = gtFwdSimProjDefinition('fwd_sim', num_dets);
+    fwd_sim = repmat(gtFwdSimDataStructureDefinition(0), [num_dets, 1]);
 
     for ii_d = conf.det_index
         uvw = gr.allblobs.detector(ii_d).uvw;
@@ -409,7 +410,7 @@ for n = first : last
             spotsCommProps, segmentedSpots, parameters);
 
         fprintf(' (Done).\n - Num spots per twin variant per detector:\n')
-        for ii_d = 1:numel(conf.det_index)
+        for ii_d = 1:num_dets
             fprintf('   + Detector: %d\n', conf.det_index(ii_d))
             for ii_t = 1:numel(fp_twins)
                 fprintf('     %d) %d\n', ii_t, ...
@@ -624,7 +625,7 @@ function [fwd_sim, bl] = forward_simulate_grain(gr, ondet, ...
 
         % in case there are several candidates, take the one with best
         % likelihood in terms of: shape, UV position, and W position.
-        ls = ls_pb .* ls_uv .* ls_w;
+        ls = sqrt(ls_pb) .* ls_uv .* ls_w;
         [~, cand_inc_index] = max(ls);
 
         fwd_sim.cands(ii).id = candidateIDs';
@@ -680,6 +681,8 @@ function twin_vars = forward_simulate_twin(gr, gr_fwd_sim, gr_proj, conf, spotsC
     sg = parameters.cryst(gr.phaseid).spacegroup;
     lp = parameters.cryst(gr.phaseid).latticepar;
 
+    num_dets = numel(conf.det_index);
+
     % Computing all twin variants for all the given pairs angle-axis
     num_angles = numel(conf.twin_angles);
     twin_vars = cell(num_angles, 1);
@@ -709,8 +712,8 @@ function twin_vars = forward_simulate_twin(gr, gr_fwd_sim, gr_proj, conf, spotsC
         twin.difspotidA = [];
         twin.difspotidB = [];
 
-        twin_fwd_sim = repmat(gtFwdSimDataStructureDefinition(0), [numel(conf.det_index), 1]);
-        twin_proj = gtFwdSimProjDefinition('fwd_sim', numel(conf.det_index));
+        twin_fwd_sim = repmat(gtFwdSimDataStructureDefinition(0), [num_dets, 1]);
+        twin_proj = gtFwdSimProjDefinition('fwd_sim', num_dets);
 
         for ii_d = conf.det_index
             om_step = gtGetOmegaStepDeg(parameters, ii_d);
@@ -921,24 +924,22 @@ end
 function spot_props = predict_spot_properties(spotsCommProps, proj_bl, gr_stats, parameters)
     bb_size_factor = parameters.fsim.bbsize_factor;
 
-    % Minimums could be wrong by up to 2 pixels (roundings)
     % 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
-    % grain volume
-    spot_props.XsizeMin = min(spot_props.XsizeMin, spotsCommProps.XsizeMin);
-    spot_props.YsizeMin = min(spot_props.YsizeMin, spotsCommProps.YsizeMin);
-
-    spot_props.XsizeMax = proj_bl.mbbsize(1) + bb_size_factor * gr_stats.bbxsstd;
-    spot_props.YsizeMax = proj_bl.mbbsize(2) + bb_size_factor * gr_stats.bbysstd;
+    % to be adusted. We chose to divide them by the same factor that then
+    % multiplies the standard deviation
+    spot_props.XsizeMin = min(proj_bl.mbbsize(1) / bb_size_factor, spotsCommProps.XsizeMin);
+    spot_props.YsizeMin = min(proj_bl.mbbsize(2) / bb_size_factor, spotsCommProps.YsizeMin);
+
+    % Even if the convex shape of the grain computed from indexter's
+    % selection had a limited angle view on the grain volume, we add some
+    % reduction of the minimum spot size by user definition
+    spot_props.XsizeMin = max(spot_props.XsizeMin - bb_size_factor * gr_stats.bbxsstd, 5);
+    spot_props.YsizeMin = max(spot_props.YsizeMin - bb_size_factor * gr_stats.bbysstd, 5);
+
+    spot_props.XsizeMax = max(proj_bl.mbbsize(1), spotsCommProps.XsizeMax) + bb_size_factor * gr_stats.bbxsstd;
+    spot_props.YsizeMax = max(proj_bl.mbbsize(2), spotsCommProps.YsizeMax) + bb_size_factor * gr_stats.bbysstd;
 
     spot_props.max_w_offset_im = spotsCommProps.max_w_offset_im;
     if (isfield(spotsCommProps, 'max_n_offset') && ~isempty(spotsCommProps.max_n_offset))

zUtil_ForwardSim/gtFsimDefaultParameters.m

@@ -29,10 +29,10 @@ function par_fsim = gtFsimDefaultParameters()
     par_fsim.save_grain      = true;
     % threshold for looking for new spots in full images
     par_fsim.thr_check       = 0.1;
-    % Selects spots for reconsruction that are in the top 4.5 sigma of the likelihood rankings
-    par_fsim.sel_lh_stdev    = 1.5;
+    % Selects spots for reconsruction that are in the top 4 sigma of the likelihood rankings
+    par_fsim.sel_lh_stdev    = 1;
     % Selects spots for reconsruction that are in the top 4 sigma of the intensity rankings
     par_fsim.sel_int_stdev   = 1;
-    % Selects spots for reconsruction that are in the top 3.5 sigma of the average intensity per pixel rankings
-    par_fsim.sel_avint_stdev = 0.5;
+    % Selects spots for reconsruction that are in the top 4 sigma of the average intensity per pixel rankings
+    par_fsim.sel_avint_stdev = 1;
 end

zUtil_ForwardSim/gtFwdSimSelectSpots.m

@@ -2,13 +2,13 @@ function [selected, goodness] = gtFwdSimSelectSpots(likelihood, intensity, avg_p
 % Flag has to be filtered by the included vector
     fsim = parameters.fsim;
     if (~isfield(fsim, 'sel_lh_stdev') || isempty(fsim.sel_lh_stdev))
-        fsim.sel_lh_stdev = 1.5;
+        fsim.sel_lh_stdev = 1;
     end
     if (~isfield(fsim, 'sel_int_stdev') || isempty(fsim.sel_int_stdev))
         fsim.sel_int_stdev = 1;
     end
     if (~isfield(fsim, 'sel_avint_stdev') || isempty(fsim.sel_avint_stdev))
-        fsim.sel_avint_stdev = 0.5;
+        fsim.sel_avint_stdev = 1;
     end
     % Let's determine the indices of the selected spots in list of included
     % spots, for different categories.
@@ -27,15 +27,12 @@ end
 
 function [selected, avg_lh] = select_best_likelihoods(likelihood, std_factor, lower_thr)
     if (~exist('lower_thr', 'var'))
-        lower_thr = 1e-1;
-    end
-    if (~exist('std_factor', 'var'))
-        std_factor = 1.5;
+        lower_thr = 1 / std_factor;
     end
 
     avg_lh = mean(likelihood);
 
-    std_dev_thr = mean(likelihood) - std_factor * std(likelihood);
+    std_dev_thr = avg_lh - std_factor * std(likelihood);
     avg_thr = avg_lh * lower_thr;
 
     selected = (likelihood > std_dev_thr) & (likelihood > avg_thr);