gtForwardSimulate_v2.m

function out = gtForwardSimulate_v2(first, last, workingdirectory, phaseID, varargin)
% GTFORWARDSIMULATE_V2  Finds missing diffraction spots and updates grain structure
%     out = gtForwardSimulate_v2(first, last, workingdirectory, [phaseID], varargin)
%     ------------------------------------------------------------------------------
%
%     INPUT:
%       first            = <int>     first grainid
%       last             = <int>     last grainid
%       workingdirectory = <string>  working directory
%       phaseID          = <int>     phase ID
%       varargin         =           parameters.fsim plus some additional
%                                    options as in gtAppFsimExtraPars():
%
%     OPTIONAL INPUT (p/v pairs):
%       updateSample     = <logical> true/false to update sample.mat and
%                                    spot2grain.mat {true}
%       updateDBtables   = <logical> true/false to update DB table with
%                                    'grainID', 'phaseID', 'flag' {false}
%       area_factor      = <double>  Weight for area criterion {1}
%       position_factor  = <double>  Weight for UV position criterion {1}
%       image_factor     = <double>  Weight for image criterion {1}
%       int_factor       = <double>  Weight for intensity criterion {0}
%
%     OPTIONAL INPUT for twin variants calculation (p/v pairs):
%       variants         = <logical> Use fsim relative to twin variants
%                                    {false}
%       plot_variants    = <logical> Plot the variants orientation wrt to
%                                    the parent orientation as XYZ crystal
%                                    axes {false}
%       num_variant      = <double>  Variant number to forward simulate
%                                    {[]}
%       twin_angle       = <double>  twin angle {[]}
%       twin_axis        = <double>  twin axis {[]}
%       twin_center      = <double>  alternative center for the twin (mm)
%                                    or the twinID if available from the
%                                    same indexed dataset {[]}
%       twin_grain       = <struct>  twin grain structure for common spot
%                                    properties and full image {[]}
%
%     OUTPUT:
%       out  = <struct>   updated grain structure. Fields:
%          .bl             = <struct N> stack of normalized difblobs, along with metadata
%          .proj.geom      = <double>   input for grain reconstruction
%                                        / forward projection using ROI
%                                        (hsize, vsize) (Nx12)
%          .proj.geom_full = <double>   input for grain reconstruction
%                                        / forward projection (full images) (Nx12)
%          .proj.geom_abs  = <double>   input for grain reconstruction
%                                        / forward projection of extinction spots (Nx12)
%          .proj.difstack  = <double>   stack of normalized difspots (XxYxN)
%          .ondet          = <double>   indices of those spots (from the list in allblobs)
%                                       which are expected on the detector
%          .uvw            = <double>   spot positions on the detector (Nx3)
%          .hklsp          = <double>   signed hk(i)l of reflection (Nx3(4))
%          .hkl            = <double>   hk(i)l family for each reflection (Nx3(4))
%          .flag           = <int>      status flag indicating result of forward simulation
%          .allblobs       = <struct>   forward simulation output produced by gtIndexFwdSimGrains
%          .used_fam       = <double>   List of unique HKL families used
%                                       (theoretical)
%          .used_famsp     = <double>   List of all the HKL families used
%                                       (theoretical)
%
%       We will classify spots with a flag  [1...6]:
%        - 1 : No intensity at expected spot position
%        - 2 : Found intensity in FULL images at predicted position
%        - 3 : Found segmented difspot at expected position - but position
%              / size do not match strict criteria
%        - 4 : Found segmented difspot in database matching the strict
%              search criteria - we will use theoretically predicted
%              angles for reconstruction
%        - 5 : Difspot is part of a pair identified by INDEXTER - use
%              Friedelpair information for reconstruction
%        - 6 : Conflict difspot (if existing)
%
%     Version 010 21-07-2014 by LNervo
%       Extracted the ASTRA geometry and the gtFwdSimFindCandidates routines

%     Version 009 14-05-2014 by LNervo
%       Added check of used hkl families
%
%     Version 008 25-04-2014 by LNervo
%       Adapted to use twin forward simulation
%
%     Version 007 Beginning-2014 by LNervo
%       Modulized, simplified, added application parameters to update
%       sample.mat and DB difspot table
%
%     Version 006 18-12-2013 by LNervo
%       Moduling a bit more...
%
%     Version 005 05/03/2013 by Nicola Vigano (nicola.vigano@esrf.fr)
%       Re-structured completely, to be more clear and easy to debug.
%
%     Version 004 21/06/2012 by AKing
%       Trivial change to make segbb zeros(1, 6)
%
%     Version 003 14-06-2012 by LNervo
%       Use gtFsimDefaultParameters


if (isdeployed)
    global GT_DB %#ok<NUSED,TLEV>
    global GT_MATLAB_HOME %#ok<NUSED,TLEV>
    load('workspaceGlobal.mat');
end

cd(workingdirectory);
gtDBConnect();

parameters = gtLoadParameters();
if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
    parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(parameters.labgeo);
end
detgeo = parameters.detgeo;
if (~isfield(parameters.labgeo, 'omstep') || isempty(parameters.labgeo.omstep))
    parameters.labgeo.omstep = gtGetOmegaStepDeg(parameters);
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 ' ...
          'ForwardSimulation to work properly. You can now create them ' ...
          'using some of the conversion functions like "gtConvertRawImages2Blobs"' ])
end

if (~isfield(parameters.fsim, 'save_grain'))
    parameters.fsim.save_grain = true;
end
if (~exist('phaseID', 'var') || isempty(phaseID))
    phaseID = 1;
end

if (isdeployed)
    first   = str2double(first);
    last    = str2double(last);
    phaseID = str2double(phaseID);
    parameters.fsim.display_figure = false;
elseif (~isempty(varargin))
    % in interactive mode we can modify parameters
    [parameters.fsim, rej_pars] = parse_pv_pairs(parameters.fsim, varargin);
end
if (~exist('rej_pars', 'var'))
    rej_pars = {};
end

indexgrainModel = fullfile(parameters.acq.dir, '4_grains', 'phase_%02d', 'index.mat');
fprintf('\nLoading indexing results...\n');
grain = [];
load(sprintf(indexgrainModel, phaseID), 'grain');
totGrains = length(grain);

if (last > totGrains)
    warning('gtForwardSimulate_v2:wrong_parameter', ...
        'Last was bigger than the total number of grains, setting it to max')
    last = totGrains;
end

difspotfile = fullfile(parameters.acq.dir, '2_difspot', '00000', 'difspot00001.edf');
if (exist(difspotfile, 'file'))
    info = edf_info(difspotfile);
else
    info = [];
end

%%% this section should be replaced with a function dealing with different
% versions of parameters names:  fed, v1, v2...

% !!! using a mean pixel size is uncertain; ideally fwd simulation should
% use mm-s

fprintf('\nLoading reflections from database..');
% Initialization of spot table and info about pairs
segmentedSpots = gtDBLoadTable([parameters.acq.name 'difspot'], 'difspotID');
pairsTable     = gtDBLoadTable(parameters.acq.pair_tablename, 'pairID');
fprintf('(Done).\n')

difspot_num = numel(segmentedSpots.difspotID); % number of difspots in database
spot2grain  = cell(difspot_num, 1); % cell array linking difspots with grainids
spot2phase  = cell(difspot_num, 1);

filetable   = [parameters.acq.name '_filetable'];
filename    = fullfile(parameters.acq.dir, '4_grains', 'spot2grain.mat');
spotTable   = GtMergeTable(filetable, filename, {'spot2grain', 'spot2phase'});

samplefile  = fullfile(parameters.acq.dir, '4_grains', 'sample.mat');
sample      = GtSample.loadFromLockedFile(filetable, samplefile);

% initialise sample structure
phase = GtPhase(parameters.cryst(phaseID).name, totGrains);
sample.phases{phaseID} = phase;

% initialize output
out = [];

if (parameters.fsim.verbose)
    disp('FwdSim parameters:')
    print_structure(parameters.fsim, 'parameters.fsim', false, true)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = first : last
    fprintf('\nRunning forward simulation for grain %d:\n', n);
    fprintf(' - Computing grain properties..');

    gr = grain{n}; % output from INDEXTER

    %%% 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);

    %%% Get the indices of all spots which fully intersect the active
    % region of the detector
    onDet = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);

    fprintf(' (Done).\n - Producing blob masks..')

    proj_bls = predict_spot_masks_from_indexed( ...
        gr, grain, segmentedSpots, spotsCommProps, onDet, parameters);

    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);

    % 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, bl] = forward_simulate_grain(gr, onDet, spotsCommProps, ...
        segmentedSpots, proj_bls, parameters, det_index);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% 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 = gtPlaceSubImage2(spotInfo.spot, full, spotInfo.bb(1), spotInfo.bb(2), 'summed');
                end
            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

    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);

    % 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

    %%% 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);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Now fill in the information for the output
    %  remove reflections for which no segmented spot has been found

    out.id        = gr.id;
    out.phaseid   = phaseID;
    out.R_vector  = gr.R_vector;
    out.center    = gr.center;
    out.allblobs  = gr.allblobs;

    out.check     = fwd_sim.check;
    out.flag      = fwd_sim.flag;
    out.pairid    = gr.pairid;
    out.spotid    = fwd_sim.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);

    % save full image
    out.full      = full;

    % check spot info
    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, :);

    found_reflections = numel(find(fwd_sim.flag > 1));
    out.completeness  = found_reflections / numspots;
    out.completeness_incl = numel(find(fwd_sim.flag >= 3)) / numspots;
    out.goodness      = goodness;
    out.goodness_sel  = mean(out.likelihood(selected));
    out.confidence    = out.goodness * out.completeness_incl;

    if (parameters.fsim.save_grain)
        fprintf(' (Done).\n - Saving ouput..');

        gtSaveGrain(phaseID, gr.id, out);
    end

    if (parameters.fsim.display_figure)
        fprintf(' (Done).\n - Displaying figure..');

        displayShowFsim(phaseID, out, rej_pars{:})
    end

    fprintf(' (Done).\n\n');

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Now update spotTables and sample
    spot2grain(out.difspotID) = {n};
    spot2phase(out.difspotID) = {phaseID};

    sample.phases{phaseID}.setR_vector(n, gr.R_vector);
    sample.phases{phaseID}.setCompleteness(n, out.completeness);
    sample.phases{phaseID}.setCenter(n, out.center); % center in pixels

    if (isempty(out.proj.selected))
        warning('gtForwardSimulate_v2:bad_value', ...
            'No spot was selected! This will make reconstruction impossible')
    elseif (numel(out.proj.selected) < 10)
        warning('gtForwardSimulate_v2:bad_value', ...
            'Too few selected spots! This will give a very bad reconstruction')
    end

    if (out.goodness < 0.25)
        warning('gtForwardSimulate_v2:very_low_goodness', ...
            [ 'Very low goodness (included %f, selected %f) indicates a' ...
              ' probably poor quality of grain reconstruction' ], ...
            out.goodness, out.goodness_sel)
    elseif (out.goodness < 0.5)
        warning('gtForwardSimulate_v2:low_goodness', ...
            [ 'Low goodness (included %f, selected %f) indicates a' ...
              ' possibly poor quality of grain reconstruction' ], ...
            out.goodness, out.goodness_sel)
    end

    if (~parameters.fsim.check_spot)
        fprintf(['Completeness estimates may be too pessimistic - activate ' ...
                 'option ''check_spot'' for more reliable results\n'])
    end

    fprintf('\n Grain number %d:\n', n);

    newspots = length(find(fwd_sim.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);
    fprintf('Goodness estimate: included %f, selected %f\n', ...
        out.goodness, out.goodness_sel);
    fprintf('%d spots will be used for reconstruction\n', length(out.proj.selected));
    fprintf(['%d new spots have been detected in addition to the %d ' ...
             'spotpairs identified by INDEXTER\n'], newspots, gr.nof_pairs);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (parameters.fsim.save_grain)
        % update difspot table in the DB
        % keep list of difspots
        spotID = out.spotid;
        spotID(out.spotid == 0) = [];
        flags = out.flag;
        flags(out.spotid == 0) = [];

        gtDBFillUpdateColumn([parameters.acq.name 'difspot'], ...
            'grainID', 'int', 'integral', 'difspotID', {out.difspotID}, n)
        gtDBFillUpdateColumn([parameters.acq.name 'difspot'], ...
            'phaseID', 'int', 'integral', 'difspotID', {out.difspotID}, phaseID)
        gtDBFillUpdateColumn([parameters.acq.name 'difspot'], ...
            'flag', 'int', 'integral', 'difspotID', {spotID}, {flags});
    end
end

if (parameters.fsim.save_grain)
    % fill spotTables
    spotTable.fillTable(spot2grain, 'spot2grain');
    spotTable.fillTable(spot2phase, 'spot2phase');
    % save sample.mat
    sample.mergeToFile(filetable, samplefile, 'sample');
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sub-functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Forward Simulation

function fwd_sim = gtFwdSimDataStructureDefinition(numspots)
    fwd_sim = struct( ...
        'flag', ones(numspots, 1), ... % see above
    	'bb', zeros(numspots, 4), ... % difspot BoundingBox from database
    	'cent', zeros(numspots, 1), ... % centroid distance from database added 06-02-2014 LNervo
    	'spotid', zeros(numspots, 1), ... % Mysql difspottable difspotID
        'check', false(numspots, 2), ...
    	'intensity', zeros(numspots, 1), ... % Integrated blob intensity
        'avg_pixel_int', zeros(numspots, 1), ... % Average pixel intensity in blob
        'likelihood', zeros(numspots, 1), ... % Blob likelihood
        'cands', struct( ...
            'id', cell(numspots, 1), ...
            'likelihood_tot', cell(numspots, 1), ...
            'likelihood_proj_mask', cell(numspots, 1), ...
            'likelihood_dev_uv', cell(numspots, 1), ...
            '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)) );
end

function [fwd_sim, bl] = forward_simulate_grain(gr, onDet, ...
        spotsCommProps, segmentedSpots, proj_bls, parameters, det_index)

    uvw = gr.allblobs.detector(det_index).uvw;

    % number of exspected diffraction spots
    numspots = numel(onDet);

    % stack of zero padded diffraction spots (projections used for 3D reconstruction)
    bl(1:numspots) = gtFwdSimBlobDefinition();
    fwd_sim = gtFwdSimDataStructureDefinition(numspots);

    % 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
    for ii = 1:numspots

        % Producing predicted spot's properties
        spot_props = predict_spot_properties(spotsCommProps, proj_bls(ii), gr.stat, parameters);

        % 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, uvw(onDet(ii), :), segmentedSpots);

        if (isempty(candidateIDs))
            % NO CANDIDATES FOUND:
            %   segmentation / overlap / intensity / ... problem ?

            if (spotsCands.found_intensity)
                % A bigger spot was found in the region.. it may be an
                % overlap
                fwd_sim.flag(ii) = 2;
            end

            % If requested, we will check into the raw images later.
            continue;
        end

        bls = gtFwdSimBuildDifstackBlobs(candidateIDs, [], parameters, ...
            spotsCommProps.stackUSize, spotsCommProps.stackVSize);

        proj_bl = proj_bls(ii * ones(numel(bls), 1));

        % Computing spot area/shape deviation/likelihood
        ls_pb = compute_proj_bbs_partial_likelihoods(bls, proj_bl);

        refl_indx = onDet(ii) * ones(numel(bls), 1);

        % Computing UV deviation/likelihood
        th_uv = uvw(refl_indx, 1:2);
        blobs_uv = [ spotsCands.CentroidX, spotsCands.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;
        ls_w = compute_w_partial_likelihoods(blob_w, th_w);

        % 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;
        [~, cand_inc_index] = max(ls);

        fwd_sim.cands(ii).id = candidateIDs';
        fwd_sim.cands(ii).likelihood_tot = ls';
        fwd_sim.cands(ii).likelihood_proj_mask = ls_pb';
        fwd_sim.cands(ii).likelihood_dev_uv = ls_uv';
        fwd_sim.cands(ii).likelihood_dev_w = ls_w';

        fwd_sim.spotid(ii) = candidateIDs(cand_inc_index);

        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 ];
        fwd_sim.bb(ii, :) = [...
            spotsCands.BoundingBoxXorigin(cand_inc_index), ...
            spotsCands.BoundingBoxYorigin(cand_inc_index), ...
            spotsCands.BoundingBoxXsize(cand_inc_index), ...
            spotsCands.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]);
        fwd_sim.intensity(ii) = bls(cand_inc_index).intensity;
        fwd_sim.avg_pixel_int(ii) = fwd_sim.intensity(ii) / sum(bls(cand_inc_index).mask(:));
        fwd_sim.likelihood(ii) = ls(cand_inc_index);

        % Found a segmented difspot at expected position:
        % we will check if it does match strict criteria
        uv_ok = fwd_sim.cent(ii, :) <= spotsCommProps.Offset;
        bb_ok = all( [...
            between(fwd_sim.bb(ii, 3), spotsCommProps.XsizeMin, spotsCommProps.XsizeMax), ...
            between(fwd_sim.bb(ii, 4), spotsCommProps.YsizeMin, spotsCommProps.YsizeMax) ] );
        fwd_sim.check(ii, :)  = [uv_ok, bb_ok];

        isSpotA = any(gr.difspotidA == fwd_sim.spotid(ii));
        isSpotB = any(gr.difspotidB == fwd_sim.spotid(ii));

        if (isSpotA || isSpotB)
            % This spot is part of a pair and has been used by INDEXTER
            fwd_sim.flag(ii) = 5;
            fwd_sim.spot_type(ii) = char('a' * isSpotA + 'b' * isSpotB);
        elseif all(fwd_sim.check(ii, :))
            fwd_sim.flag(ii) = 4; % This flag is kinda obsolete
        else
            fwd_sim.flag(ii) = 3;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Support Functions

function full = sfBuildFull(spotid, parameters, spotsCommProps, bb)
% full = sfBuildFull(spotid, parameters, spotsCommProps, bb)
%
% spotid         : spotid from grain
% parameters     : parameters.mat
% spotsCommProps : spots common properties
% bb             : bounding boxes from DB

    full = zeros(parameters.detgeo.detsizev, parameters.detgeo.detsizeu, 'single');

    for ii = 1:numel(spotid)
        [spot, ~] = gtFwdSimBuildSpot(spotid(ii), parameters, spotsCommProps.Xsize, spotsCommProps.Ysize, bb(ii, :));

        full = gtPlaceSubImage2(spot, full, bb(ii, 1), bb(ii, 2), 'summed');
    end
end

function [diff_beam_dir, omega] = sfGetFriedelPair(pairsTable, spotID, isSpotA)
% [diff_beam_dirs, omega] = sfGetFriedelPair(pairsTable, spotID, isSpotA)

    for ii = numel(spotID):-1:1
        if (isSpotA(ii))
            idx = find(pairsTable.difAID == spotID(ii));
            omega(ii) = pairsTable.omega(idx);
            diff_beam_dir(ii, :) = [pairsTable.ldirX(idx), pairsTable.ldirY(idx), pairsTable.ldirZ(idx)];
        else
            idx = find(pairsTable.difBID == spotID(ii));
            omega(ii) = pairsTable.omegaB(idx);
            % in this case we have to reverse the line direction
            diff_beam_dir(ii, :) = -[pairsTable.ldirX(idx), pairsTable.ldirY(idx), pairsTable.ldirZ(idx)];
        end
    end
end

function displayShowFsim(phaseID, grain, varargin)
% displayShowFsim(phaseID, grain, varargin)

    if isempty(varargin)
        [vars_pars, vars_list] = gtAppFsimDefaultPars();
        ind = findValueIntoCell(fieldnames(vars_pars), {'fsimtype', 'f_title'});
        vars_list(ind(:, 1), 4) = {2};

        vars_pars = gtModifyStructure(vars_pars, vars_list, 1, 'Options for displaying the Forward Sim:');
        varargin = struct2pv(vars_pars);
    end
    gtShowFsim(grain, phaseID, varargin{:});
end

function proj_bls = predict_spot_masks_from_indexed(gr, grain_INDX, ...
    segmentedSpots, spotsCommProps, onDet, parameters)

    % First argument: included_INDXTR, might be needed in the future for
    % some reason...
    [~, difspotID_INDXTR] = gtFwdSimGetIndexterReflections(gr, grain_INDX{gr.id}, segmentedSpots);

    % Producing blobs' boundingboxes (as by dilated segmentation info)
    bls_INDXTR = gtFwdSimBuildDifstackBlobs(difspotID_INDXTR, [], ...
        parameters, spotsCommProps.stackUSize, spotsCommProps.stackVSize);
    % INDEXTER Data
    BBus = cat(1, bls_INDXTR(:).mbbu);
    BBvs = cat(1, bls_INDXTR(:).mbbv);
    BBsizes = cat(1, bls_INDXTR(:).mbbsize);
    BBs = [BBus(:, 1), BBvs(:, 1), BBsizes(:, 1:2)];

    Ws = segmentedSpots.CentroidImage(difspotID_INDXTR) .* parameters.labgeo.omstep;

    proj_bls = gtFwdSimPredictProjectedGrainBBox(gr, BBs, Ws, onDet, parameters);

    % Producing original blobs' boundingboxes (as by original segmentation
    % info)
    BBs = [...
        segmentedSpots.BoundingBoxXorigin(difspotID_INDXTR), ...
        segmentedSpots.BoundingBoxYorigin(difspotID_INDXTR), ...
        segmentedSpots.BoundingBoxXsize(difspotID_INDXTR), ...
        segmentedSpots.BoundingBoxYsize(difspotID_INDXTR) ];

    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;

    % 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);

    % 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;

    spot_props.thr_max_offset = spotsCommProps.thr_max_offset;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Likelihood computation

function ls = compute_uv_partial_likelihoods(blob_uv, th_uv, proj_bls)
    proj_bbox(:, [1 3]) = cat(1, proj_bls(:).bbuim);
    proj_bbox(:, [2 4]) = cat(1, proj_bls(:).bbvim);

    blob_expected_sizes = proj_bbox(:, [3 4]) - proj_bbox(:, [1 2]) + 1;

    ls = (th_uv - blob_uv) ./ max(blob_expected_sizes ./ 10, 10);
    ls = exp(- (ls .^ 2) );
    ls = ls(:, 1) .* ls(:, 2);
end

function ls = compute_w_partial_likelihoods(blob_w, th_w)
    ls = 2 .* (th_w - blob_w);
    ls = exp(- (ls .^ 2));
end

function ls = compute_proj_bbs_partial_likelihoods(bls, proj_bls)
    bls_bboxes(:, [1, 3]) = cat(1, bls(:).bbuim);
    bls_bboxes(:, [2, 4]) = cat(1, bls(:).bbvim);

    proj_bls_bboxes(:, [1, 3]) = cat(1, proj_bls(:).bbuim);
    proj_bls_bboxes(:, [2, 4]) = cat(1, proj_bls(:).bbvim);

    cont_bbxes = [ ...
        min(bls_bboxes(:, 1:2), proj_bls_bboxes(:, 1:2)), ...
        max(bls_bboxes(:, 3:4), proj_bls_bboxes(:, 3:4)) ...
        ];

    cont_bbxes = [cont_bbxes(:, 1:2), (cont_bbxes(:, 3:4) - cont_bbxes(:, 1:2) + 1)];

    bls_shifts = bls_bboxes(:, 1:2) - cont_bbxes(:, 1:2);
    bls_shifts(:, 3) = 0;
    proj_bls_shifts = proj_bls_bboxes(:, 1:2) - cont_bbxes(:, 1:2);
    proj_bls_shifts(:, 3) = 0;

    num_bls = numel(proj_bls);
    mismatch_more = zeros(num_bls, 1);
    mismatch_less = zeros(num_bls, 1);
    proj_count = zeros(num_bls, 1);

    for ii = 1:num_bls
        cont_size = cont_bbxes(ii, 3:4);

        proj_count(ii) = sum(proj_bls(ii).mask(:));

        bl_mask = bls(ii).mask;
        bl_mask = sum(bl_mask, 3);
        bl_mask = bl_mask > 0;

        proj_bl_mask = logical(proj_bls(ii).mask);

        bl_mask = gtPlaceSubVolume( ...
            zeros(cont_size, 'uint8'), bl_mask, bls_shifts(ii, :));
        proj_bl_mask = gtPlaceSubVolume( ...
            zeros(cont_size, 'uint8'), proj_bl_mask, proj_bls_shifts(ii, :));

        overlap = bl_mask & proj_bl_mask;

        bl_mask(overlap) = 0;
        proj_bl_mask(overlap) = 0;

        mismatch_more(ii) = sum(bl_mask(:));
        mismatch_less(ii) = sum(proj_bl_mask(:));
    end

    % 3 and 2 are temporary
    ls = (1 .* mismatch_more + 1 .* mismatch_less) ./ proj_count;

    ls = exp(- (ls .^ 2));
end