Skip to content
Snippets Groups Projects
Commit 7e722e19 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

FwdSim: adapt to some recent changes 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 666f2be5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 44 additions and 418 deletions
zUtil_Deformation/gtDefSyntheticGrainCreate.m
+ 16
− 12
View file @ 7e722e19
......@@ -245,18 +245,22 @@ for n = 1:num_dets
vol_size(3) = max(vol_size(3), 8);
end
grain.proj(n).stack = permute(stack, [1 3 2]);
grain.proj(n).num_cols = size(grain.proj(n).stack, 3);
grain.proj(n).num_rows = size(grain.proj(n).stack, 1);
grain.proj(n).bl = bl;
grain.proj(n).vol_size_x = vol_size(2);
grain.proj(n).vol_size_y = vol_size(1);
grain.proj(n).vol_size_z = vol_size(3);
grain.proj(n).centerpix = gtGeoSam2Sam(grain.center, samgeo, recgeo(n), 0);
grain.proj(n).shift = gtFwdSimComputeVolumeShifts(grain.proj(n), parameters, vol_size);
grain.proj(n).ondet = ondet;
grain.proj(n).included = included;
grain.proj(n).selected = selected;
proj = gtFwdSimProjDefinition();
proj.stack = permute(stack, [1 3 2]);
proj.num_cols = size(proj.stack, 3);
proj.num_rows = size(proj.stack, 1);
proj.bl = bl;
proj.vol_size_x = vol_size(2);
proj.vol_size_y = vol_size(1);
proj.vol_size_z = vol_size(3);
proj.centerpix = gtGeoSam2Sam(grain.center, samgeo, recgeo(n), 0);
proj.shift = gtFwdSimComputeVolumeShifts(grain.proj(n), parameters, vol_size);
proj.ondet = ondet;
proj.included = included;
proj.selected = selected;
grain.proj(n) = proj;
fprintf('\b\b(%f seconds)\n - True blob dimensions:\n', toc(c))
fprintf(' %3d) %4d %4d %4d\n', [(1:nbl)', vertcat(bl(:).mbbsize)]')
......
zUtil_Deformation/test_generation/Gt6DTwinTestCreator.m
+ 14
− 14
View file @ 7e722e19
......@@ -17,31 +17,31 @@ classdef Gt6DTwinTestCreator < Gt6DTestCreator
twin_angle = 60;
spacegroup = 225;
g_twins = gtTwinOrientations(grain.R_vector, twin_angle, twin_axis, spacegroup, lp, false);
twin_vars = gtTwinOrientations(grain.R_vector, twin_angle, twin_axis, spacegroup, lp, false);
for jj = length(g_twins):-1:1
for jj = length(twin_vars):-1:1
disp('CHANGE THE SIGN OF THE R-VECTOR')
my_g_twins(jj).R_vector = -g_twins(jj).R_vector;
my_g_twins(jj).g = g_twins(jj).g;
my_g_twins(jj).center = grain.center;
my_twin_vars(jj).R_vector = -twin_vars(jj).R_vector;
my_twin_vars(jj).g = twin_vars(jj).g;
my_twin_vars(jj).center = grain.center;
end
g_twins = my_g_twins;
twin_vars = my_twin_vars;
for jj = 1:length(g_twins)
grain.R_vector = -g_twins(jj).R_vector;
for jj = 1:length(twin_vars)
grain.R_vector = -twin_vars(jj).R_vector;
grain = gtCalculateGrain(grain, parameters);
g_twins(jj).allblobs = grain.allblobs;
twin_vars(jj).allblobs = grain.allblobs;
end
grain.g_twins = g_twins; % this will be the parent, then all the twins
grain.twin_vars = twin_vars; % this will be the parent, then all the twins
% parent grain - reset the normal stuff in grain
grain.R_vector = R_vector;
grain = gtCalculateGrain(grain, parameters);
grain.proj.geom = obj.buildGrainGeom(parameters, grain, pixelsize);
for ii = 1:numel(grain.g_twins)
grain.g_twins(ii).proj.geom = obj.buildGrainGeom(parameters, grain.g_twins(ii), pixelsize);
for ii = 1:numel(grain.twin_vars)
grain.twin_vars(ii).proj.geom = obj.buildGrainGeom(parameters, grain.twin_vars(ii), pixelsize);
end
% If we want to show the result of forward simulation
......@@ -70,7 +70,7 @@ classdef Gt6DTwinTestCreator < Gt6DTestCreator
geoms = zeros(numCols, numVols-1);
for ii = 1:numVols-1
twinGeom = grain.g_twins(ii).proj.geom(:, 1:3);
twinGeom = grain.twin_vars(ii).proj.geom(:, 1:3);
indexesTwin = ones(size(twinGeom, 1), 1);
for jj = 1:numCols
......@@ -150,7 +150,7 @@ classdef Gt6DTwinTestCreator < Gt6DTestCreator
jj = 1;
parentGeom = grain.proj.geom;
twinGeom = grain.g_twins(jj).proj.geom;
twinGeom = grain.twin_vars(jj).proj.geom;
commonIndx = find(geoms(:, jj) ~= 0);
onlyParentIndx = setdiff(1:size(parentGeom, jj), commonIndx);
......
zUtil_Drawing/gtPlotMultiScatter3DMatrix.m
+ 0
− 1
View file @ 7e722e19
......@@ -223,7 +223,6 @@ end
end % end of function
%%
% sub-functions
function sfSetInfo(hObj, xyz)
......
zUtil_Drawing/gtPlotMultiUVW.m
+ 1
− 1
View file @ 7e722e19
......@@ -187,7 +187,7 @@ for jj=1:length(grainIDs)
%gtPlotMultiScatter3D({grains}, 'ids', {grainIDs(jj)}, 'var', {'allblobs'}, 'var2', {'uvw'}, ...
% 'labels',{['fwd2 grain #' num2str(grainIDs(jj))]}, varargin{:})
case 'twin-parent'
gtPlotMultiScatter3DMatrix({grains{grainIDs(jj)}.fwd2.allblobs.detector(1).uvw, grains{grainIDs(jj)}.g_twins(fsimID).fwd2.allblobs.detector(1).uvw}, ...
gtPlotMultiScatter3DMatrix({grains{grainIDs(jj)}.fwd2.allblobs.detector(1).uvw, grains{grainIDs(jj)}.twin_vars(fsimID).fwd2.allblobs.detector(1).uvw}, ...
'labels',{['parent grain #' num2str(grainIDs(jj))],['twin variant ' num2str(fsimID)]}, varargin{:})
set(f(jj),'Name',['grain #' num2str(grainIDs(jj)) ' twin variant ' num2str(fsimID)])
case 'indexfsim-fwd2'
......
zUtil_ForwardSim/gtFwdSimPredictGrainVolumeSize.m deleted 100644 → 0
+ 0
− 24
View file @ 666f2be5
function vol_size = gtFwdSimPredictGrainVolumeSize(gr, parameters, use_polyhedron)
if (~exist('use_polyhedron', 'var') || isempty(use_polyhedron))
use_polyhedron = false;
end
if (isfield(gr, 'om_exp'))
omegas = gr.om_exp;
else
omegas = gr.omega;
end
if (isfield(gr, 'bb_exp'))
bb = gr.bb_exp;
else
bb = gr.bb;
end
stackUSize = gr.proj.num_cols;
stackVSize = gr.proj.num_rows;
selected = gr.selected;
vol_size = gtFwdSimPredictVolumeSize(gr.center, omegas(selected), ...
bb(selected, :), parameters, stackUSize, stackVSize, use_polyhedron);
end
zUtil_GUI/gtGuiGrainMontage.m
+ 13
− 8
View file @ 7e722e19
......@@ -58,15 +58,20 @@ function selected = gtGuiGrainMontage(gr, selected, use_blobs, det_ind)
if (use_blobs && ~isfield(gr.proj(det_ind), 'bl'))
gr.proj(det_ind).bl = gr.bl;
end
if (isfield(gr, 'spotid'))
spot_ids = gr.spotid(included);
if (isfield(gr, 'fwd_sim'))
spot_ids = gr.fwd_sim(det_ind).spotid(included);
intensities = gr.fwd_sim(det_ind).intensity(included);
else
spot_ids = zeros(numel(included), 1);
end
if (isfield(gr, 'intensity'))
intensities = gr.intensity;
else
intensities = zeros(numel(included), 1);
if (isfield(gr, 'spotid'))
spot_ids = gr.spotid(included);
else
spot_ids = zeros(numel(included), 1);
end
if (isfield(gr, 'intensity'))
intensities = gr.intensity;
else
intensities = zeros(numel(included), 1);
end
end
% Condition to exit from the gui
......
zUtil_Twins/gtTwinForwardSimulate.m deleted 100644 → 0
+ 0
− 358
View file @ 666f2be5
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment