-
Nicola Vigano authored
Signed-off-by:Nicola Vigano <nicola.vigano@esrf.fr>
Nicola Vigano authoredSigned-off-by:
Gt6DTwinTestCreator.m 15.63 KiB
classdef Gt6DTwinTestCreator < Gt6DTestCreator
methods (Access = public)
function obj = Gt6DTwinTestCreator(volsNum)
obj = obj@Gt6DTestCreator(volsNum);
end
function grain = generateTwinsGeometry(obj, parameters, R_vector, center)
% To be cleaned!
grain = [];
grain.phaseid = 1;
grain.R_vector = R_vector;
pixelsize = mean([parameters.labgeo.pixelsizeu, parameters.labgeo.pixelsizev]);
grain.center = center .* pixelsize;
lp = parameters.cryst(grain.phaseid).latticepar;
twin_axis = [1 1 1];
twin_angle = 60;
spacegroup = 225;
twin_vars = gtTwinOrientations(grain.R_vector, twin_angle, twin_axis, spacegroup, lp, false);
for jj = length(twin_vars):-1:1
disp('CHANGE THE SIGN OF THE R-VECTOR')
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
twin_vars = my_twin_vars;
for jj = 1:length(twin_vars)
grain.R_vector = -twin_vars(jj).R_vector;
grain = gtCalculateGrain(grain, parameters);
twin_vars(jj).allblobs = grain.allblobs;
end
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.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
% fig = figure();
% ax = axes('Parent', fig);
% hold(ax, 'on');
% for ii = 1:size(grain.proj.geom, 1)
% vect = grain.proj.geom(ii, 1:3);
% quiver3(ax, 0, 0, 0, vect(1), vect(2), vect(3));
% end
% hold(ax, 'off');
% grid(ax, 'on');
% zlim(ax, [-1, 1]);
end
function geoms = matchParentTwinGeometries(obj, grain)
% Now we need to find out which geometries are in common with
% the parent
numVols = numel(obj.singleVolumes);
detectorDistance = grain.proj.geom(1, 4:6);
detectorDistance = norm(detectorDistance);
parentGeom = grain.proj.geom(:, 1:3);
numCols = size(parentGeom, 1);
geoms = zeros(numCols, numVols-1);
for ii = 1:numVols-1
twinGeom = grain.twin_vars(ii).proj.geom(:, 1:3);
indexesTwin = ones(size(twinGeom, 1), 1);
for jj = 1:numCols
diffGeoms = parentGeom(jj .* indexesTwin, :) - twinGeom;
normDiffs = sqrt(sum(diffGeoms .* diffGeoms, 2));
[val, indx] = min(normDiffs);
angle = 2 * asind( val / 2);
if (angle < 5e-1)
geoms(jj, ii) = indx;
fprintf('Found: norm %f, angle %f (deg), offset %f\n', ...
val, angle, val * detectorDistance)
end
end
indexes = find(geoms(:, ii) ~= 0);
[parentGeom(indexes, :), twinGeom(geoms(indexes, ii), :)]
end
end
function [indexes, finalGeoms] = createTwins(obj, metaInfoGrains, grains)
if (length(grains) < 2)
error('CS:wrong_parameter', ...
'At least two twin geometries are needed');
end
if (length(grains) > 2)
error('CS:wrong_parameter', ...
'More than two is not supported, yet.');
end
indexes = [];
[indexes.common_id, indexes.common_indx1, indexes.common_indx2] ...
= intersect(grains{1}.difspotID(metaInfoGrains{1}.selectedDiffspots), ...
grains{2}.difspotID(metaInfoGrains{2}.selectedDiffspots));
[indexes.only_1_id, indexes.only_1_indx] ...
= setdiff(grains{1}.difspotID, grains{2}.difspotID);
[indexes.only_2_id, indexes.only_2_indx] ...
= setdiff(grains{2}.difspotID, grains{1}.difspotID);
finalGeoms = cell(2, 1);
finalGeoms{1} = grains{1}.proj.geom(indexes.common_indx1, :);
finalGeoms{2} = grains{1}.proj.geom(indexes.common_indx1, :);
modifiedOnly1 = grains{1}.proj.geom(indexes.only_1_indx, :);
finalGeoms{1} = [finalGeoms{1}; modifiedOnly1];
modifiedOnly1(:, 4:6) = obj.findPerpendicularDirection(modifiedOnly1(:, 4:6));
finalGeoms{2} = [finalGeoms{2}; modifiedOnly1];
% diffCenter = grains{2}.center - grains{1}.center;
modifiedOnly2 = grains{2}.proj.geom(indexes.only_2_indx, :);
% modifiedOnly2(:, 4:6) = modifiedOnly2(:, 4:6) ...
% + diffCenter(ones(size(modifiedOnly2, 1), 1), :);
finalGeoms{2} = [finalGeoms{2}; modifiedOnly2];
modifiedOnly2(:, 4:6) = obj.findPerpendicularDirection(modifiedOnly2(:, 4:6));
finalGeoms{1} = [finalGeoms{1}; modifiedOnly2];
fprintf('Common spots (num %04d):\n', length(indexes.common_id))
fprintf(' %d\n', indexes.common_id);
fprintf('First volume''s spots (num %04d):\n', length(indexes.only_1_id))
fprintf(' %d\n', indexes.only_1_id);
fprintf('Second volume''s spots (num %04d):\n', length(indexes.only_2_id))
fprintf(' %d\n', indexes.only_2_id);
end
function [detectorImages, trueDetectorImages] ...
= getSyntheticTwinDetectorImages(obj, grain, geoms, maxIndivSpots)
% GETSYNTHETICTWINDETECTORIMAGES
% [~, finalGeoms] = obj.createTwins(metaInfoGrains, grains);
numVols = numel(obj.singleVolumes);
finalGeoms = cell(1, numVols);
if (numVols > 2)
error('Gt6DTwinTestCreator:wrong_parameter', ...
'More than 2 volumes, not supported yet')
end
jj = 1;
parentGeom = grain.proj.geom;
twinGeom = grain.twin_vars(jj).proj.geom;
commonIndx = find(geoms(:, jj) ~= 0);
onlyParentIndx = setdiff(1:size(parentGeom, jj), commonIndx);
onlyTwinIndx = setdiff(1:size(twinGeom, jj), geoms(commonIndx, jj));
if (exist('maxIndivSpots', 'var'))
if (numel(onlyParentIndx) > maxIndivSpots)
onlyParentIndx = onlyParentIndx(randperm(numel(onlyParentIndx), maxIndivSpots));
end
if (numel(onlyTwinIndx) > maxIndivSpots)
onlyTwinIndx = onlyTwinIndx(randperm(numel(onlyTwinIndx), maxIndivSpots));
end
end
perpendParent = [ ...
obj.findPerpendicularDirection(twinGeom(onlyTwinIndx, 1:3)), ...
twinGeom(onlyTwinIndx, 4:end) ...
];
finalGeoms{1} = [ ...
parentGeom(onlyParentIndx, :); ...
parentGeom(commonIndx, :); ...
perpendParent ...
];
perpendTwin = [ ...
obj.findPerpendicularDirection(parentGeom(onlyParentIndx, 1:3)), ...
parentGeom(onlyParentIndx, 4:end) ...
];
finalGeoms{jj+1} = [ ...
perpendTwin; ...
parentGeom(commonIndx, :); ...
twinGeom(onlyTwinIndx, :) ...
];
for ii = 1:numel(finalGeoms)
obj.proj.addProjectionGeometry( finalGeoms{ii} );
end
% Visualization of diffraction directions
f = figure();
ax = axes('parent', f);
hold(ax, 'on');
zeroCols = zeros(numel(onlyParentIndx), 1);
quiver3(ax, zeroCols, zeroCols, zeroCols, ...
parentGeom(onlyParentIndx, 1), parentGeom(onlyParentIndx, 2), parentGeom(onlyParentIndx, 3))
zeroCols = zeros(numel(commonIndx), 1);
quiver3(ax, zeroCols, zeroCols, zeroCols, ...
parentGeom(commonIndx, 1), parentGeom(commonIndx, 2), parentGeom(commonIndx, 3))
zeroCols = zeros(numel(onlyTwinIndx), 1);
quiver3(ax, zeroCols, zeroCols, zeroCols, ...
twinGeom(onlyTwinIndx, 1), twinGeom(onlyTwinIndx, 2), twinGeom(onlyTwinIndx, 3))
hold(ax, 'off');
grid(ax, 'on'); zlim(ax, [-1 1]);
[detectorImages, trueDetectorImages] = obj.getDetectorImages();
end
function [detectorImages, trueDetectorImages] ...
= getTrueTwinDetectorImages(obj, metaInfoGrains, grains)
%
numGrains = length(grains);
volumeOversize = 1;
spotsOversize = 1.2;
if (numGrains < 2)
error('CS:wrong_parameter', ...
'At least two twin geometries are needed');
end
if (numGrains > 2)
error('CS:wrong_parameter', ...
'More than two is not supported, yet.');
end
% Let's create initial geometry, and decide which spots go in
[indexes, finalGeoms] = obj.createTwins(metaInfoGrains, grains);
% Shift geometries (relative to centers), create main voloume
% and give it a decent size
obj.fixCenterToGeometriesAndAdd(finalGeoms, grains, volumeOversize);
% First crop/pad the bigger images
for ii = 1:numGrains
sizes(ii, :) = size(grains{ii}.proj.stack);
end
maxSpotBBox = [ max([sizes(:, 1); round(obj.proj.proj_size(1) * spotsOversize)]), ...
max([sizes(:, 3); round(obj.proj.proj_size(2) * spotsOversize)])];
offsets(:, 1) = floor((maxSpotBBox(1) - sizes(:, 1)) / 2);
offsets(:, 2) = zeros(numGrains, 1);
offsets(:, 3) = floor((maxSpotBBox(2) - sizes(:, 3)) / 2);
modifiedStacks = cell(numGrains, 1);
for ii = 1:numGrains
modifiedStacks{ii} = zeros( maxSpotBBox(1), ...
size(grains{ii}.proj.stack, 2), ...
maxSpotBBox(2));
modifiedStacks{ii} = gtPlaceSubVolume(modifiedStacks{ii}, ...
grains{ii}.proj.stack, ...
offsets(ii, :), ...
0, 'assign');
end
detectorImages = [modifiedStacks{1}(:, indexes.common_indx1, :), ...
modifiedStacks{1}(:, indexes.only_1_indx, :), ...
modifiedStacks{2}(:, indexes.only_2_indx, :)];
trueDetectorImages = permute(detectorImages, [1 3 2]);
end
end
methods (Access = protected)
function proj_geom = buildGrainGeom(obj, parameters, grain, pixelsize)
u = grain.allblobs.uv(:, 1);
v = grain.allblobs.uv(:, 2);
spotsCommProps = getDiffractionSpotsCommonProperties(obj);
radiuses = sqrt((u - parameters.labgeo.detrefu).^2 + (v - parameters.labgeo.detrefv).^2);
%%% Get the indices of all spots which fully intersect the active
% region of the detector
onDet = Gt6DTwinTestCreator.findSpotsOnDetector(parameters, spotsCommProps, u, v, radiuses);
rotComp = gtMathsRotationMatrixComp(parameters.labgeo.rotdir, 'row');
getRotationTensor = @(x)gtMathsRotationTensor(-x, rotComp);
detposcorner = gtGeoDetOrig(parameters.labgeo) ./ pixelsize;
proj_geom = [];
for ii = numel(onDet):-1:1
currentReflection = onDet(ii);
% use theoretically predited directions (default)
omega = grain.allblobs.omega(currentReflection);
omegaRotTens = getRotationTensor(omega);
pl = grain.allblobs.pl(currentReflection, :);
theta = grain.allblobs.theta(currentReflection);
% theoretically predicted beam direction
diffBeamDir = parameters.labgeo.beamdir * omegaRotTens + 2 * sind(theta) * pl;
% position of Detector
% detPos = detposcorner ...
% + spotsCommProps.stackUSize/2 * parameters.labgeo.detdiru * omegaRotTens ...
% + spotsCommProps.stackVSize/2 * parameters.labgeo.detdirv * omegaRotTens;
% detPos = detPos * omegaRotTens - grain.center ./ pixelsize;
detPos = diffBeamDir * norm(detposcorner);
detdiru = parameters.labgeo.detdiru * omegaRotTens;
detdirv = parameters.labgeo.detdirv * omegaRotTens;
proj_geom(ii, :) = [diffBeamDir, detPos, detdiru, detdirv];
end
end
function spotsProps = getDiffractionSpotsCommonProperties(obj)
%%% Get the average difspot X and Y sizes for the current spots in the grain,
spotsProps = [];
spotsProps.Xsize = mean(size(obj.fullVol)); % average X Size of difspots
spotsProps.Ysize = mean(size(obj.fullVol)); % 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 - 1);
spotsProps.XsizeMax = spotsProps.Xsize + 1;
spotsProps.YsizeMin = max(5, spotsProps.Ysize - 1);
spotsProps.YsizeMax = spotsProps.Ysize + 1;
% 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 = 1;
diagonal = round(norm(size(obj.fullVol)));
% take the biggest spot and add a relative extra margin (multiplicative factor oversize)
spotsProps.stackUSize = diagonal;
% ...same for vertical direction
spotsProps.stackVSize = diagonal;
end
end
methods (Access = protected, Static)
function ondet = findSpotsOnDetector(parameters, dif_props, u, v, radius)
%%% Get the indices of all spots which fully intersect the active region of the detector
props = [ 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];
ondet = find(all(props, 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
end
end