gtForwardSimulate_v2.m 29.71 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 number {1}
%
% OPTIONAL INPUT (varargin):
% bb_size = <int> size of search area ([] will use default projection size)
% assemble_figure = <logical> collects all difspots and assembles them into a fullimage
% save_grain = <logical> flag to save grain_###.mat {true}
%
% OUTPUT:
% out = <struct> : updated grain structure. Fields:
% - difstack <double XxYxN> stack of normalized difspots
% - proj_geom <double Nx12> input for grain reconstruction
% / forward projection using roi (hsize, vsize)
% - proj_geom_full <double Nx12> input for grain reconstruction
% / forward projection (full images)
% - proj_geom_abs <double Nx12> input for grain reconstruction
% / forward projection of extinction spots
% - ondet <double N> indices of those spots (from the list in allblobs)
% which are expected on the detector
% - uvw <double Nx3> spot positions on the detector
% - hklsp <double Nx3(4)> signed hk(i)l of reflection
% - flag <int Nx3> status flag indicating result of forward simulation
% - allblobs <struct> forward simulation output produced by gtFedGenerateRandomGrain
%
% We will classify spots with a flag [1...5]:
% - 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
%
% Version 003 14-06-2012 by LNervo
% Use gtFsimDefaultParameters
%
% Version 004 21/06/2012 by AKing
% Trivial change to make segbb zeros(1,6)
%
% Version 005 05/03/2013 by Nicola Vigano (nicola.vigano@esrf.fr)
% Re-structured completely, to be more clear and easy to debug.
%
% SUB-FUNCTIONS:
%[sub]- sfGetRoi
if isdeployed
global GT_DB
global GT_MATLAB_HOME
load('workspaceGlobal.mat');
end
cd(workingdirectory);
gtDBConnect();
parameters = [];
load('parameters.mat');
acq = parameters.acq;
samgeo = parameters.samgeo;
labgeo = parameters.labgeo;
if isfield(parameters,'fsim')
fsim = parameters.fsim;
else
fsim = gtFsimDefaultParameters();
end
if ~isfield(fsim, 'save_grain')
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);
fsim.display_figure = false;
else
% in interactive mode we can modify parameters
fsim = parse_pv_pairs(fsim, varargin);
end
phaseID_str = sprintf('phase_%02d', phaseID);
phases_num = length(parameters.cryst); % number of phases in the material
grainfilename = fullfile(parameters.acq.dir, '4_grains', phaseID_str, 'index.mat');
fprintf('\nLoading indexing results...\n');
grains = load(grainfilename);
totGrains = length(grains.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', '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...
detdiru0 = labgeo.detdiru; % will use detdiru0 to distinguish from rotated version of this vector
detdirv0 = labgeo.detdirv;
% !!! using a mean pixel size is uncertain; ideally fwd simulation should
% use mm-s
pixelsize = mean([labgeo.pixelsizeu, labgeo.pixelsizev]);
% detector reference position in pixels
detrefpos = labgeo.detrefpos ./ pixelsize;
% the (0,0) corner of the detector in Lab system (edge of detector pixel)
detposcorner = gtGeoDetOrig(labgeo) ./ pixelsize;
rotComp = gtMathsRotationMatrixComp(labgeo.rotdir, 'row');
getRotationTensor = @(x)gtMathsRotationTensor(-x, rotComp);
fprintf('\nLoading reflections from database...\n');
% Initialization of spot table and info about pairs
segmentedSpots = initializeSpotsRecords(parameters.acq.name);
pairsTable = getPairsTable(parameters.acq.pair_tablename);
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;
for n = first : last
fprintf('\nRunning forward simulation for grain %d...\n',n);
gr = grains.grain{n}; % output from INDEXTER
% graincenter contains COM coordinates in pixels
graincenter_pix = gr.center ./ pixelsize;
%%% Predict spot positions
gr = gtCalculateGrain(gr, parameters);
u = gr.allblobs.uv(:, 1);
v = gr.allblobs.uv(:, 2);
%%% 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 = getDiffractionSpotsCommonProperties(gr, fsim, detrefpos);
% radial distances of spot center from detector center
% !!! to be corrected: detrefu and detrefv can be anywhere, not
% necessarily in the center
radiuses = sqrt((u - labgeo.detrefu).^2 + (v - labgeo.detrefv).^2);
%%% Get the indices of all spots which fully intersect the active
% region of the detector
onDet = findSpotsOnDetector(parameters, spotsCommProps, u, v, radiuses);
% Maximum offset (in number of images) from the predicted omega position
% allow for a 1/sin(eta) dependency of this value - limited to a
% maximum angular offset of fsim.MaxOmegaOffset degrees: this could be
% improved...
maxOffsetThrs = getThresholdsMaxOffset(fsim, gr, getOmegaStepSize(parameters));
image_nums = round(gr.allblobs.omega / getOmegaStepSize(parameters));
% Converts from SoA to AoS
spotsPositions = buildSpotsPositions(onDet, u, v, image_nums, maxOffsetThrs);
numspots = numel(onDet); % number of exspected diffraction spots
flag = ones(numspots, 1); % see above
% stack of zero padded diffraction spots (projections used for 3D reconstruction)
difstack = zeros(spotsCommProps.stackHSize, numspots, spotsCommProps.stackVSize);
bb = zeros(numspots, 4); % difspot BoundingBox from database
spotid = zeros(numspots, 1); % Mysql difspottable difspotID
check = zeros(numspots, 2);
omega = zeros(numspots, 1); % rotation angle in degrees
% Vector used for ROI grain reconstruction in ASTRA (grain shifted to center of roi volume)
proj_geom = zeros(numspots, 12);
% Vector describing full projection geometry (full images, grain at nominal position in sample volume)
proj_geom_full = zeros(numspots, 12);
% Vector describung absorption image projection geometry
proj_geom_abs = zeros(numspots, 12);
U = zeros(numspots, 1); % COM u positions of difspots from difspottable
V = zeros(numspots, 1); % COM v positions of difspots from difspottable
full = zeros(acq.ydet, acq.xdet);
tempIsSpotA = false(numspots, 1);
%%% 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..
for ii = 1 : numspots
reflectionPos = spotsPositions(ii);
% 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 = findCandidates(spotsCommProps, reflectionPos, segmentedSpots);
if isempty(candidateIDs)
% NO CANDIDATES FOUND:
% segmentation / overlap / intensity / ... problem ?
% If requested, we will check into the raw images later.
continue;
end
% in case there are several candidates, take the one which has
% the Area close to the expected Area. This selection "rule"
% could be improved... err appears to be not respected!
[~, area_index] = sort(abs( ...
segmentedSpots.BoundingBoxXsize(candidateIDs) .* segmentedSpots.BoundingBoxYsize(candidateIDs) ...
- spotsCommProps.Xsize * spotsCommProps.Ysize));
[~, image_index] = sort(abs(segmentedSpots.CentroidImage(candidateIDs) - reflectionPos.img_num));
matchedID = candidateIDs(image_index(1));
U(ii) = segmentedSpots.CentroidX(matchedID);
V(ii) = segmentedSpots.CentroidY(matchedID);
bb(ii, :) = [...
segmentedSpots.BoundingBoxXorigin(matchedID), segmentedSpots.BoundingBoxYorigin(matchedID), ...
segmentedSpots.BoundingBoxXsize(matchedID), segmentedSpots.BoundingBoxYsize(matchedID) ];
% Found a segmented difspot at expected position - we will check if it does match strict criteria
cent_ok = sqrt( (U(ii) - reflectionPos.u).^2 + (V(ii) - reflectionPos.v).^2 ) <= spotsCommProps.Offset;
size_ok = all( [...
between(bb(ii, 3), spotsCommProps.XsizeMin, spotsCommProps.XsizeMax), ...
between(bb(ii, 4), spotsCommProps.YsizeMin, spotsCommProps.YsizeMax) ] );
check(ii, :) = [cent_ok, size_ok];
isSpotA = ~isempty(find(gr.difspotidA == matchedID, 1));
isSpotB = ~isempty(find(gr.difspotidB == matchedID, 1));
if (isSpotA || isSpotB)
% This spot is part of a pair and has been used by INDEXTER
flag(ii) = 5;
tempIsSpotA(ii) = isSpotA;
elseif all(check(ii,:))
flag(ii) = 4;
else
flag(ii) = 3;
end
spotid(ii) = matchedID;
end
%%% Re-check flag == 1, to see if there's intensity at the predicted
% position, if requesred.
if (fsim.check_spot)
toBeChecked = find(flag == 1);
for ii = reshape(toBeChecked, [1 numel(toBeChecked)]);
reflectionPos = spotsPositions(onDet(ii));
[full, found_intensity] = checkSpotInRawImages(full, ...
reflectionPos, fsim, spotsCommProps, parameters);
if (found_intensity)
% Found intensity in FULL images at predicted position
flag(ii) = 2;
end
end
end
% Reflections that have been matched to segmented spots
matched = find(flag >= 3);
%%% We now build the diffraction stack
for ii = reshape(matched, [1 numel(matched)])
spot = buildSpot(spotid(ii), parameters, info, spotsCommProps, bb(ii, :));
if (fsim.assemble_figure)
full = gtPlaceSubImage2(spot, full, bb(ii, 1), bb(ii, 2), 'summed');
end
% Essentially zero pads the spot to make it fit into ASTRA's
% diffraction stack
shifts = getDiffStackShifts(spotsCommProps, bb(ii, :));
spot = gtPlaceSubImage2(spot, ...
zeros(spotsCommProps.stackVSize, spotsCommProps.stackHSize), ...
shifts.h + 1, shifts.v + 1, 'summed');
difstack(:, ii, :) = spot';
end
%%% We now Build ASTRA's geometry
for ii = reshape(matched, [1 numel(matched)])
% If available, we want to use the angles as determined from
% the Friedel pairs - interrogate spotpair table
%
% diffBeamDir: direction of diffracted beam used in ASTRA
if (flag(ii) == 5)
[diffBeamDir, omega(ii)] = getFriedelPair(pairsTable, spotid(ii), tempIsSpotA(ii));
omegaRotTens = getRotationTensor(omega(ii));
else
% we found a segmentedd spot - can use the theoretically
% predicted direction or apparent direction (based on COM of
% grain and spot) for backprojection
% theoretically prediceted omega angle
if (fsim.use_th)
currentReflection = onDet(ii);
% use theoretically predited directions (default)
omega(ii) = gr.allblobs.omega(currentReflection);
omegaRotTens = getRotationTensor(omega(ii));
pl = gr.allblobs.pl(currentReflection, :);
theta = gr.allblobs.theta(currentReflection);
% theoretically predicted beam direction
diffBeamDir = labgeo.beamdir * omegaRotTens + 2 * sind(theta) * pl;
else
% use experimentally determined omega angles
omega(ii) = segmentedSpots.CentroidImage(spotid(ii)) * getOmegaStepSize(parameters);
omegaRotTens = getRotationTensor(omega(ii));
% projection direction as inferred from grain / spot
% center positions
diffBeamDir = gtGeoDiffVecInSample([U(ii), V(ii)], omega(ii), gr.center, labgeo, samgeo);
end
end
%%% Now we actually build the geometry
shifts = getDiffStackShifts(spotsCommProps, bb(ii, :));
% position of Detector
detPos = detposcorner ...
+ (bb(ii, 1) - shifts.h - 1 + spotsCommProps.stackHSize/2) * detdiru0 ...
+ (bb(ii, 2) - shifts.v - 1 + spotsCommProps.stackVSize/2) * detdirv0;
detAbs = detposcorner ...
+ (acq.bb(1) - 1 + acq.bb(3)/2) * detdiru0 ...
+ (acq.bb(2) - 1 + acq.bb(4)/2) * detdirv0;
detdiru = detdiru0 * omegaRotTens;
detdirv = detdirv0 * omegaRotTens;
% take into account the coordinate shift corresponding to the offset
% of the reconstructed volume
detPos = detPos * omegaRotTens - graincenter_pix;
detAbs = detAbs * omegaRotTens;
absBeamDir = labgeo.beamdir * omegaRotTens;
% Compact version, since they're all <1x3> vectors, so we can read
% them easily
proj_geom(ii, :) = [diffBeamDir, detPos, detdiru, detdirv];
proj_geom_abs(ii, :) = [absBeamDir, detAbs, detdiru, detdirv];
% Calculate the projection geometry for the full images as well
% labgeo.detrefu: detector u (line direction) reference point (center) in lab
% labgeo.detrefv: detector v (column direction )reference point (center) in lab
detposFull = detposcorner + labgeo.detrefu * detdiru0 + labgeo.detrefv * detdirv0;
proj_geom_full(ii, :) = proj_geom(ii,:);
proj_geom_full(ii, 4:6) = (detposFull * omegaRotTens);
end % loop over predicted spot positions
included = find(flag >= fsim.included);
% indices of the selected spots in list of included spots
selected = find(flag(included) >= fsim.selected);
%%% Now fill in the information for the output
% remove reflections for which no segmented spot has been found
out.selected = selected; % vector of logicals with pre-selected spots active
out.id = gr.id;
out.phaseid = phaseID;
out.center = graincenter_pix;
out.ondet = onDet;
out.allblobs = gr.allblobs;
out.full = full;
out.flag = flag;
out.check = check;
out.pairid = gr.pairid;
out.difspotID = spotid(included);
out.omega = omega(included);
% signed hk(i)l indices of reflection
hklsp = gr.allblobs.hklsp(out.ondet, :);
out.hklsp = hklsp(included, :);
out.uvw = gr.allblobs.uvw(out.ondet(included), :);
out.uv_exp = [U(included), V(included)];
out.bb = bb(included, :);
found_reflections = length(find(flag > 1));
out.completeness = found_reflections / numspots;
out.proj.stack = difstack(:, included, :);
out.proj.geom = proj_geom(included, :);
out.proj.geom_full = proj_geom_full(included, :);
out.proj.geom_abs = proj_geom_abs(included, :);
out.proj.num_rows = spotsCommProps.stackVSize;
out.proj.num_cols = spotsCommProps.stackHSize;
% We should in the future handle properly vertical detector (general geometry)
oversizeVol = fsim.oversizeVol / fsim.oversize;
out.proj.vol_size_x = round(spotsCommProps.stackHSize * oversizeVol);
out.proj.vol_size_y = round(spotsCommProps.stackHSize * oversizeVol);
out.proj.vol_size_z = round(spotsCommProps.stackVSize * oversizeVol);
out.proj.num_iter = parameters.rec.num_iter; out.vol = [];
out.segbb = zeros(1,6);
out.threshold = 0;
out.shift(1) = round((acq.bb(3)-out.proj.vol_size_x)/2 + out.center(1)) + 1;
out.shift(2) = round((acq.bb(3)-out.proj.vol_size_y)/2 + out.center(2)) + 1;
out.shift(3) = round((acq.bb(4)-out.proj.vol_size_z)/2 + out.center(3)) + 1;
if (fsim.display_figure)
out.full = full;
gtShowFsim(out)
end
if (fsim.save_grain)
filename = fullfile(acq.dir, '4_grains', phaseID_str, ...
sprintf('grain_%04d.mat', gr.id));
save(filename, '-struct', 'out');
end
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);
sample.phases{phaseID}.setSelectedDiffspots(n, out.selected);
if (isempty(out.selected))
warning('gtForwardSimulate_v2:bad_value', ...
'No spot was selected! This will make reconstruction impossible')
elseif (numel(out.selected) < 10)
warning('gtForwardSimulate_v2:bad_value', ...
'Too few selected spots! This will give a very bad reconstruction')
end
if (~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(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('%d spots will be used for reconstruction\n', length(out.selected));
fprintf(['%d new spots have been detected in addition to the %d ' ...
'spotpairs identified by INDEXTER\n'], newspots, gr.nof_pairs);
end
spotTable.fillTable(spot2grain, 'spot2grain');
spotTable.fillTable(spot2phase, 'spot2phase');
sample.mergeToFile(filetable, samplefile, 'sample');
end % end of function
%%% sub-functions
function [spot, m] = getRawRoi(start_image, end_image, acq, bb)
% GTFORWARDSIMULATE/SFGETROI Sums a roi in raw images
nimages = end_image - start_image + 1;
spot = zeros(bb(4), bb(3), nimages);
m = zeros(nimages, 1);
centerx = round((bb(3)/2 -5) : (bb(3)/2 +5));
centery = round((bb(4)/2 -5) : (bb(4)/2 +5));
if all([bb(1) > 0, bb(2) > 0, bb(1)+bb(3) < acq.xdet, bb(2)+bb(4) < acq.ydet])
fullImgsDir = fullfile(acq.dir, '1_preprocessing', 'full'); indexes = start_image : end_image;
filename = fullfile(fullImgsDir, sprintf('full%04d.edf', indexes(1)));
info = edf_info(filename);
for img_n = 1:nimages
fullImgName = fullfile(fullImgsDir, sprintf('full%04d.edf', indexes(img_n)));
spot(:, :, img_n) = edf_read(fullImgName, bb, [], info);
m(img_n) = gtImgMeanValue(spot(centery, centerx, img_n));
end
end
end
function spotsRecords = initializeSpotsRecords(dataset_name)
difspottable = [dataset_name 'difspot'];
% Load difspot data in 'allspots' / spotpair data in 'allpairs'
query = [ ...
'SELECT difspotID, CentroidX, CentroidY, CentroidImage, '...
'BoundingBoxXsize, BoundingBoxYsize, BoundingBoxXorigin, ' ...
'BoundingBoxYorigin'...
' FROM ' difspottable ];
[difspotIDs, CentroidXs, CentroidYs, CentroidImages, BoundingBoxXsizes, ...
BoundingBoxYsizes, BoundingBoxXorigins, BoundingBoxYorigins] = mym(query);
% create vectors of size (numspots,1) (these will have some empty elements since
% not all difspotIDs exist (some are filtered out because fo size
% restrictions, for example)
numspots = max(difspotIDs);
spotsRecords = [];
spotsRecords.difspotID = NaN(numspots, 1);
spotsRecords.CentroidX = NaN(numspots, 1);
spotsRecords.CentroidY = NaN(numspots, 1);
spotsRecords.CentroidImage = NaN(numspots, 1);
spotsRecords.BoundingBoxXsize = NaN(numspots, 1);
spotsRecords.BoundingBoxYsize = NaN(numspots, 1);
spotsRecords.BoundingBoxXorigin = NaN(numspots, 1);
spotsRecords.BoundingBoxYorigin = NaN(numspots, 1);
orderedPositions = 1:numspots;
spotsRecords.difspotID(difspotIDs) = difspotIDs;
spotsRecords.CentroidX(difspotIDs) = CentroidXs(orderedPositions);
spotsRecords.CentroidY(difspotIDs) = CentroidYs(orderedPositions);
spotsRecords.CentroidImage(difspotIDs) = CentroidImages(orderedPositions);
spotsRecords.BoundingBoxXsize(difspotIDs) = BoundingBoxXsizes(orderedPositions);
spotsRecords.BoundingBoxYsize(difspotIDs) = BoundingBoxYsizes(orderedPositions);
spotsRecords.BoundingBoxXorigin(difspotIDs) = BoundingBoxXorigins(orderedPositions);
spotsRecords.BoundingBoxYorigin(difspotIDs) = BoundingBoxYorigins(orderedPositions);
end
function pairsTable = getPairsTable(pairtable_name)
query = [ ...
'SELECT ldirX, ldirY, ldirZ, omega, omegaB, difAID, difBID' ...
' FROM ' pairtable_name ];
[pairsTable.ldirX, pairsTable.ldirY, pairsTable.ldirZ, ...
pairsTable.omegaA, pairsTable.omegaB, ...
pairsTable.difAID, pairsTable.difBID ] = mym(query);
end
function thrs_max_offset = getThresholdsMaxOffset(fsim, gr, omstep)
maxOff = fsim.MaxOmegaOffset;
omegarange = fsim.omegarange;
etas = sind(gr.allblobs.eta);
spotsOff = omegarange ./ abs( etas );
thrs_max_offset = round(min(maxOff, spotsOff) / omstep);
end
function spotsProps = getDiffractionSpotsCommonProperties(gr, fsim, detrefpos_pix)
%%% Get the average difspot X and Y sizes for the current spots in the grain, spotsProps = [];
spotsProps.Xsize = gr.stat.bbxsmean; % average X Size of difspots
spotsProps.Ysize = gr.stat.bbysmean; % average Y Size of difspots
% define theshold values which will be used during the database search later
% - limit to 5 pixels minimum width
spotsProps.XsizeMin = max(5, spotsProps.Xsize - fsim.bbsize_factor * gr.stat.bbxsstd);
spotsProps.XsizeMax = spotsProps.Xsize + fsim.bbsize_factor * gr.stat.bbxsstd;
spotsProps.YsizeMin = max(5, spotsProps.Ysize - fsim.bbsize_factor * gr.stat.bbysstd);
spotsProps.YsizeMax = spotsProps.Ysize + fsim.bbsize_factor * gr.stat.bbysstd;
% allow for a lateral deviation of spot postions based on grain orientation statics
% dif_props.Offset = fsim.Rdist_factor * gr.stat.Rdiststd * norm(detrefpos_pix)
% !!! detrefpos is not necessarily the distance; depends on definition by user
spotsProps.Offset = fsim.Rdist_factor * deg2rad(gr.stat.dangstd) * norm(detrefpos_pix);
grainBBox = [max(gr.bbxs), max(gr.bbys)];
% take the biggest spot and add a relative extra margin (multiplicative factor oversize)
spotsProps.stackHSize = round(grainBBox(1) * fsim.oversize);
% ...same for vertical direction
spotsProps.stackVSize = round(grainBBox(2) * fsim.oversize);
end
function ondet = findSpotsOnDetector(parameters, dif_props, u, v, radius)
%%% Get the indices of all spots which fully intersect the active region of the detector
ondet = find(all([u > dif_props.Xsize/2, u < parameters.acq.xdet - dif_props.Xsize/2 , ...
v > dif_props.Ysize/2, v < parameters.acq.ydet - dif_props.Ysize/2, ...
radius < parameters.acq.maxradius], 2));
% Exclude spots which overlap with segmentation boundingbox
seg = find(all([( (u + dif_props.Xsize/2 > parameters.seg.bbox(1)) ...
& (u - dif_props.Xsize/2 < parameters.seg.bbox(1) + parameters.seg.bbox(3))), ...
( (v + dif_props.Ysize/2 > parameters.seg.bbox(2)) ...
& (v - dif_props.Ysize/2 < parameters.seg.bbox(2) + parameters.seg.bbox(4)))], 2));
ondet = setdiff(ondet, seg); % indices of spots intersection active area of detector
end
function [full, found_intensity] = checkSpotInRawImages(full, spotPos, ...
fsim, spotsCommProps, parameters)
% number of images summed up when searching raw data for a missing spot
imagerange = round(fsim.omegarange / getOmegaStepSize(parameters));
% in case bb_size for search BoundingBox (used for check of
% completeness) is empty, use the same size as the difspot stack
if isempty(fsim.bb_size)
bbhsize = spotsCommProps.stackHSize;
bbvsize = spotsCommProps.stackVSize;
else
bbhsize = fsim.bb_size(1);
bbvsize = fsim.bb_size(2);
end
% Adjust search bbox in order to avoid reading information
% stemming from regions outside the image border
bbox_ustart = round(spotPos.u - bbhsize/2);
bbox_vstart = round(spotPos.v - bbvsize/2);
bbox_uend = min(parameters.acq.xdet, bbox_ustart + bbhsize);
bbox_vend = min(parameters.acq.ydet, bbox_vstart + bbvsize);
bbox(1) = max(1, bbox_ustart);
bbox(2) = max(1, bbox_vstart);
bbox(3) = (bbox_uend - bbox_ustart + 1);
bbox(4) = (bbox_vend - bbox_vstart + 1);
start_image = max(0, spotPos.img_num - imagerange);
end_image = min(spotPos.img_num + imagerange, 2*parameters.acq.nproj - 1);
% read in small image stack at predicted position read ROI from
% raw images - meanVals are the mean value for each image in the stack
[raw_img, meanVals] = getRawRoi(start_image, end_image, parameters.acq, bbox); [intensity, ~] = max(meanVals);
meanIntensity = sum(raw_img(:)) / (spotsCommProps.Xsize * spotsCommProps.Xsize * spotsCommProps.Ysize);
raw_img = sum(raw_img, 3) / meanIntensity;
if (fsim.assemble_figure)
full = gtPlaceSubImage2(raw_img, full, bbox(1), bbox(2), 'summed');
end % app.assemble_figure
% Checks if found intensity in FULL images at predicted position
found_intensity = (intensity > parameters.seg.thr_grow_low);
end
function [diffBeamDir, omega] = getFriedelPair(pairsTab, spotID, isSpotA)
if (isSpotA)
idx = find(pairsTab.difAID == spotID);
omega = pairsTab.omegaA(idx);
diffBeamDir = [pairsTab.ldirX(idx), pairsTab.ldirY(idx), pairsTab.ldirZ(idx)];
else
idx = find(pairsTab.difBID == spotID);
omega = pairsTab.omegaB(idx);
% in this case we have to reverse the line direction
diffBeamDir = -[pairsTab.ldirX(idx), pairsTab.ldirY(idx), pairsTab.ldirZ(idx)];
end
end
function spotsPositions = buildSpotsPositions(ondet, u, v, image_nums, thrs_max_offset)
spotsPositions = arrayfun(@(x)assignEntry(x, u, v, image_nums, thrs_max_offset), ondet);
function spotPos = assignEntry(currentSpotID, u, v, image_nums, thrs_max_offset)
spotPos = [];
spotPos.u = u(currentSpotID);
spotPos.v = v(currentSpotID);
spotPos.img_num = image_nums(currentSpotID);
spotPos.thr_max_offset = thrs_max_offset(currentSpotID);
end
end
function spot = buildSpot(candidateID, parameters, info, spotsCommProps, bb)
info.dim_1 = bb(3);
info.dim_2 = bb(4);
spot = gtGetSummedDifSpot(candidateID, parameters, 1, info);
% here we force the intensity to be the same in all spots -
volume = spotsCommProps.Xsize * spotsCommProps.Xsize * spotsCommProps.Ysize;
spot = spot / sum(spot(:)) * volume;
end
function candidateIDs = findCandidates(spotsProps, reflPos, segmentedSpots)
image_ok = find(...
(segmentedSpots.CentroidImage >= (reflPos.img_num - reflPos.thr_max_offset) ) ...
& (segmentedSpots.CentroidImage <= (reflPos.img_num + reflPos.thr_max_offset) ));
spotsXsize = abs(segmentedSpots.CentroidX(image_ok) - reflPos.u);
spotsYsize = abs(segmentedSpots.CentroidY(image_ok) - reflPos.v);
maxXsize = max(spotsProps.Xsize, segmentedSpots.BoundingBoxXsize(image_ok) - spotsProps.Xsize) / 2;
maxYsize = max(spotsProps.Ysize, segmentedSpots.BoundingBoxYsize(image_ok) - spotsProps.Ysize) / 2;
intensity_ok = (spotsXsize < maxXsize) & (spotsYsize < maxYsize);
candidateIDs = image_ok(intensity_ok);
end
function omstep = getOmegaStepSize(parameters)
omstep = 180 / parameters.acq.nproj;
end
function shifts = getDiffStackShifts(spotsCommProps, bb)
shifts.h = round((spotsCommProps.stackHSize - bb(3)) / 2);
shifts.v = round((spotsCommProps.stackVSize - bb(4)) / 2);end