function grain = gtTwinForwardSimulate(phaseID, grain, twin_center, varargin)
% GTTWINFORWARDSIMULATE Finds diffraction spots in the database relative to the twin variant
% grain = gtTwinForwardSimulate(phaseID, grain, twin_center, varargin)
% --------------------------------------------------------------------------------------
%
% INPUT:
% phaseID = <int> phase number {1}
% grain = <struct> grain structure of interest
%
% OPTIONAL INPUT (varargin):
% save_grain = <logical> flag to save or not the twins_####.mat
% variantIDs = <double> array of ID for variants {all}
%
% OUTPUT:
% grain = cell structure with fields for each variant:
% - 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
% - status <string 5> explanation string corresponding to flag
%
% Version 001 16-10-2012 by LNervo
%
% SUB-FUNCTIONS:
%[sub]- sfGetRoi
if isdeployed
global GT_DB
global GT_MATLAB_HOME
load('workspaceGlobal.mat');
end
gtDBConnect();
parameters = [];
load('parameters.mat');
acq = parameters.acq;
samgeo = parameters.samgeo;
labgeo = parameters.labgeo;
fsim = parameters.fsim;
app.save_grain = false;
app.variantIDs = [];
if ~exist('phaseID','var') || isempty(phaseID)
phaseID = 1;
end
if ~exist('twin_center','var')
twin_center = [];
end
if isdeployed
phaseID = str2double(phaseID);
fsim.display_figure = false;
else
[app,rej_pars] = parse_pv_pairs(app, varargin);
end
fsim.check_spot = true;
fsim.display_figure = true;
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...
detrefu = labgeo.detrefu; % detector u (line direction) reference point (center) in lab
detrefv = labgeo.detrefv; % detector v (column direction )reference point (center) in lab
detdiru0 = labgeo.detdiru; % will use detdiru0 to distinguish from rotated version of this vector
detdirv0 = labgeo.detdirv;
rotdir = labgeo.rotdir;
beamdir = labgeo.beamdir;
pixelsize = mean([labgeo.pixelsizeu, labgeo.pixelsizev]);
detrefpos = labgeo.detrefpos ./ pixelsize; % detector reference position in pixels
detposcorner = gtGeoDetOrig(labgeo) ./ pixelsize; % the (0,0) corner of the detector in Lab system (edge of detector pixel)
detrefpos_pix = parameters.labgeo.detrefpos ./ pixelsize;
rotcomp = gtMathsRotationMatrixComp(rotdir,'row');
omstep = 180 / acq.nproj;
imagerange = round(fsim.omegarange/omstep); % number of images summed up when searching raw data for a missing spot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load info from database
fprintf('\nLoading reflections from database...\n');
segmentedSpots = gtDBLoadTable([parameters.acq.name 'difspot'], 'difspotID');
fprintf('Done!\n')
graincenter = grain.center ./ pixelsize; % graincenter contains COM coordinates in pixels
spotsCommProps = getDiffractionSpotsCommonProperties(grain, parameters.fsim, detrefpos_pix, omstep);
% calc radius of u,v coordinate
getradius = @(u,v) sqrt((u - labgeo.detrefu).^2 + (v - labgeo.detrefv).^2);
if ~isempty(twin_center) && isRowVectorWithLength(twin_center,3)
grain.center_parent = grain.center;
grain.center = twin_center;
fprintf('Using twin center (%f %f %f) for simulation of grain %d\n',grain.center,grain.id)
fprintf('parent grain center (%f %f %f)\n',grain.center_parent)
end
if ~isfield(grain, 'g_twins') || ( isfield(grain, 'g_twins') )
for kk = 1:length(grain.g_twins)
grain = gtCalculateGrain(grain, parameters, kk, 'variants', true, varargin{:});
end
end
if isempty(app.variantIDs)
app.variantIDs = 1:length(grain.g_twins);
end
gr = grain;
% loop over the variants
for kk = app.variantIDs
fprintf('\nRunning forward simulation for twin variant %d for grain %d ...\n',kk,grain.id);
variants = grain.g_twins(kk);
% assume same center for parent and twin
graincenter = grain.center ./ pixelsize; % graincenter contains COM coordinates in pixels
out.g = variants.g;
out.R_vector = variants.R_vector;
%%% Get the average difspot X and Y sizes for the current spots in the grain,
% define theshold values which will be used during the database search later
dif_Xsize = mean(segmentedSpots.BoundingBoxXsize);
dif_Ysize = mean(segmentedSpots.BoundingBoxYsize);
% take the biggest spot and add a relative extra margin (multiplicative factor oversize)
hsize = round(max(segmentedSpots.BoundingBoxXsize)*fsim.oversize);
% ...same for vertical direction
vsize = round(max(segmentedSpots.BoundingBoxYsize)*fsim.oversize);
u = variants.allblobs.detector(1).uvw(:, 1);
v = variants.allblobs.detector(1).uvw(:, 2);
uvw = variants.allblobs.detector(1).uvw;
radiuses = getradius(uvw(:,1),uvw(:,2));
ondet = findSpotsOnDetector(spotsCommProps, uvw(:,1), uvw(:,2), radiuses, parameters.acq, parameters.seg.bbox);
uvw = uvw(ondet, :);
% number of exspected diffraction spots
numspots = numel(ondet);
grain.spotid = zeros(numspots,1);
% 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'
flag = zeros(numspots, 1); % see above
difstack = zeros(hsize, numspots, vsize); % stack of zero padded diffraction spots (projections used for 3D reconstruction)
bb = zeros(numspots, 4); % difspot BoundingBox from database
status = cell(numspots, 1); % string corresponding to flag (summarizing result of forward simulation)
spotid = zeros(numspots, 1); % Mysql difspottable difspotID
check = zeros(numspots, 2);
hklsp = zeros(numspots,size(variants.allblobs.hklsp,2)); % signed hk(i)l indices of reflection
omega = zeros(numspots, 1); % rotation angle in degrees
proj_geom = zeros(numspots, 12); % Vector used for ROI grain reconstruction in ASTRA (grain shifted to center of roi volume)
proj_geom_full = zeros(numspots, 12); % Vector describing full projection geometry (full images, grain at nominal position in sample volume)
proj_geom_abs = zeros(numspots, 12); % Vector describung absorption image projection geometry
r = zeros(numspots, 3); % direction of diffracted beam used in ASTRA
r_th = zeros(numspots, 3); % theoretically predicted direction of diffracted beam (based on R_vector)
r_ex = zeros(numspots, 3); % estimated direction of diffracted beam, based on grain & difspot COM positions
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.xdet, acq.ydet);
%%% Main loop: examine all hklsp and search corresponding spots / intensity on the detector..
for ii = 1 : numspots
currentSpotID = ondet(ii);
pl = variants.allblobs.pl(currentSpotID, :);
theta = variants.allblobs.theta(currentSpotID);
om_th = variants.allblobs.omega(currentSpotID);
image_n = round(om_th/omstep);
eta = variants.allblobs.eta(currentSpotID);
hklsp(ii,:) = variants.allblobs.hklsp(currentSpotID, :);
% 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...
thr_max_offset = round( fsim.omegarange/(abs(sind(eta))) /omstep );
% 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
%if isempty(candidateIDs)
% no candidates have been found: segmentation / overlap / intensity / ... problem ?
flag(ii) = 1; % No intensity at expected spot position
if fsim.check_spot
% 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 = hsize;
bbvsize = vsize;
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(u(currentSpotID) - bbhsize/2);
bbox_vstart = round(v(currentSpotID) - bbvsize/2);
bbox_uend = min(acq.xdet, bbox_ustart + bbhsize);
bbox_vend = min(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, image_n - imagerange);
end_image = min(image_n + imagerange, 2*acq.nproj - 1);
% read in small image stack at predicted position read ROI from
% raw images - m is the mean value for each image in the stack
[roi, m] = sfGetRoi(start_image, end_image, acq, bbox);
[intensity, ~] = max(m);
%test = sum(roi,3)/sum(roi(:))*dif_Xsize*dif_Xsize*dif_Ysize;
test = sum(roi,3);
if (intensity > 0.2)
flag(ii) = 2; % Found intensity in FULL images at predicted position
end
if (fsim.assemble_figure)
full = gtPlaceSubImage2(test, full, bbox(1), bbox(2), 'summed');
end % app.assemble_figure
end % app.check_spot
% use theoretically predited directions (default)
omega(ii) = om_th;
Rottens = gtMathsRotationTensor(-omega(ii), rotcomp);
% theoretically predicted beam direction
r(ii,:) = beamdir*Rottens + 2*sind(theta)*pl;
end % loop over predicted spot positions
found_reflections = length(find(flag > 1));
completeness = found_reflections / numspots;
newspots = length(find(flag == 2));
included = find(flag == 2);
selected = find(flag(included) == 3); % indices of the selected spots in list of included spots
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', grain.id);
fprintf(['Forward simulation found intensity at %d out of %d expected ' ...
'spot positions \n'], found_reflections, numspots);
fprintf('Completeness estimate: %f\n', completeness);
fprintf('%d spots will be used for reconstruction\n', length(selected));
%%% Now fill in the information for the output
% remove reflections for which no segmented spot has been found
out.included = flag(included) == 3;
out.selected = flag(included) == 4;%>= fsim.selected; % vector of logicals with pre-selected spots active
out.id = grain.id;
out.parent_R_vector = grain.R_vector;
out.phaseid = phaseID;
out.center = graincenter*pixelsize;
out.ondet = ondet;
out.allblobs = variants.allblobs;
out.full = full;
out.flag = flag;
out.status = status;
out.check = check;
out.spotid = grain.spotid;
out.difspotID = spotid(included);
out.omega = omega(included);
out.hklsp = hklsp(included,:);
out.uvw = variants.allblobs.detector(1).uvw(ondet(included),:);
out.uv_exp = [U(included), V(included)];
out.bb = bb(included,:);
out.r_th = r_th(included,:);
out.completeness = completeness;
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 = vsize;
out.proj.num_cols = hsize;
out.proj.vol_size_x = round(hsize / sqrt(2));
out.proj.vol_size_y = round(hsize / sqrt(2));
out.proj.vol_size_z = round(vsize / sqrt(2));
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)-hsize)/2 + graincenter(1)) + 1;
out.shift(2) = round((acq.bb(3)-hsize)/2 + graincenter(2)) + 1;
out.shift(3) = round((acq.bb(4)-vsize)/2 + graincenter(3)) + 1;
if fsim.display_figure
out.full = full;
gtShowFsim(out,out.phaseid,'clims',[-200 1000],'f_title',sprintf('twin variant %d for grain %d ...\n',kk,grain.id))
end
grain.g_twins(kk).out = out;
end % loop over the variants
%
% if app.save_grain
% filename = fullfile(acq.dir, '4_grains', sprintf('phase_%02d',phaseID), ...
% sprintf('twins_%04d.mat', grain.id));
% save(filename, 'grain');
% end
end % end of function
%%
% sub-functions
function [spot, m] = sfGetRoi(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