-
Nicola Vigano authored
(multi-detector not fully implemented yet!) Signed-off-by:Nicola Vigano <nicola.vigano@esrf.fr>
Nicola Vigano authored(multi-detector not fully implemented yet!) Signed-off-by:
gtForwardSimulate_v2.m 40.57 KiB
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_angles = <double> twin angle {[]}
% twin_axes = <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;
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
conf = struct( ...
'det_index', [], ...
'twin_angles', zeros(0, 1), ...
'twin_axes', zeros(0, 3) );
if (~exist('rej_pars', 'var'))
rej_pars = {};
end
[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');
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
% Storing information about the grain in sample.mat
sample.phases{phaseID}.setR_vector(n, gr.R_vector);
sample.phases{phaseID}.setCenter(n, gr.center);
%%% 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]);
for ii_d = conf.det_index
uvw = gr.allblobs.detector(ii_d).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, ii_d);
%%% 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, ii_d);
fprintf(' (Done).\n - Preallocating structs..')
if (parameters.fsim.assemble_figure)
full = zeros(detgeo(ii_d).detsizev, detgeo(ii_d).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(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
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reflections that have been matched to segmented spots
included = find(fwd_sim(ii_d).flag >= 3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% We now build the diffraction stack
fprintf(' (Done).\n - Building difstack..');
difstack = gtFwdSimBuildDifstackSpots(fwd_sim(ii_d).spotid, ...
included, parameters, ...
spotsCommProps.stackUSize, spotsCommProps.stackVSize, ...
spotsCommProps.Xsize, spotsCommProps.Ysize, ...
fwd_sim(ii_d).bb, info);
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 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
% 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
%%% 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
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 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.fwd_sim = fwd_sim;
out.check = fwd_sim(1).check;
out.flag = fwd_sim(1).flag;
out.pairid = gr.pairid;
out.spotid = fwd_sim(1).spotid;
% only included spots info
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(1).check_spots;
% diffraction geometry
out.proj = proj;
if (~isempty(conf.twin_angles))
out.twin_vars = fp_twins;
else
out.twin_vars = [];
end
found_reflections = numel(find(fwd_sim(1).flag > 1));
out.completeness = found_reflections / 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;
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}.setCompleteness(n, out.completeness);
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(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);
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)), ...
'gv_verts', [] );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, gr, onDet(ii), segmentedSpots, parameters, det_index);
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
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_twin = gtFwdSimFindSpotsOnDetector(spotsCommProps, uvw, parameters);
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_twin, 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;
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
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, det_ind)
samgeo = parameters.samgeo;
labgeo = parameters.labgeo;
detgeo = parameters.detgeo(det_ind);
omstep = gtGetOmegaStepDeg(parameters, det_ind);
% 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) .* omstep;
gr_center = gr.center(ones(numel(difspotID_INDXTR), 1), :);
projvecs = [segmentedSpots.CentroidX(difspotID_INDXTR), segmentedSpots.CentroidY(difspotID_INDXTR)];
projvecs = gtGeoDet2Lab(projvecs, detgeo, false);
projvecs = gtGeoLab2Sam(projvecs, Ws, labgeo, samgeo, 0, 0) - gr_center;
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)
BBs = [...
segmentedSpots.BoundingBoxXorigin(difspotID_INDXTR), ...
segmentedSpots.BoundingBoxYorigin(difspotID_INDXTR), ...
segmentedSpots.BoundingBoxXsize(difspotID_INDXTR), ...
segmentedSpots.BoundingBoxYsize(difspotID_INDXTR) ];
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;
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;
proj_bls(ii_b).mask = proj_bls_original(ii_b).mask;
end
end
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);
% 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), :);
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 = gtFwdSimProjectVertices2Det(new_gr, verts_rec, new_gr_ondet, parameters, det_ind);
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 = 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;
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;
proj_bls(ii_b).mask = proj_bls_original(ii_b).mask;
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.max_w_offset_im = spotsCommProps.max_w_offset_im;
if (isfield(spotsCommProps, 'max_n_offset') && ~isempty(spotsCommProps.max_n_offset))
spot_props.max_n_offset_deg = spotsCommProps.max_n_offset_deg;
spot_props.max_t_offset_deg = spotsCommProps.max_t_offset_deg;
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
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/4 makes up for the recent change to the use of a bigger estimate of
% the grain volume
ls = (1 .* mismatch_more + 3 / 4 .* mismatch_less) ./ proj_count;
ls = exp(- (ls .^ 2));
end