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

Phantom generation: converted to multi-detector and changed output format 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent ed472b83
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
7_fed2/gtFedFwdProjExact.m
+ 12
− 7
View file @ 70d12926
function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters)
function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num)
if (~exist('det_num', 'var'))
det_num = 1;
end
fprintf('Forward projection (%s):\n', mfilename)
nbl = length(bl);
nv = size(gv.ind, 2);
detgeo = parameters.detgeo(det_num);
labgeo = parameters.labgeo;
samgeo = parameters.samgeo;
......@@ -25,7 +30,7 @@ uvw = cell(1, nbl);
fprintf(' - Computing indices and bbsizes: ')
t = tic();
% !!! parfor
for ii = 1:nbl
num_chars = fprintf('%03d/%03d', ii, nbl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -60,8 +65,8 @@ for ii = 1:nbl
gvcs_lab = gtGeoSam2Lab(gvcs', rot_s2l, labgeo, samgeo, false);
ucbl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
labgeo.detrefpos', labgeo.detnorm', labgeo.Qdet, labgeo.detrefuv, ...
om, labgeo.omstep);
detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
[detgeo.detrefu, detgeo.detrefv]', om, labgeo.omstep);
% Detector coordinates U,V in blob
bbort = bbor(ii, :)';
......@@ -126,12 +131,12 @@ end
fprintf(' Maximum BBox size: [(%d, %d, %d), (%d, %d, %d)]\n', max_bbox_proj_size(1, :), max_bbox_proj_size(2, :));
fprintf(' Maximum Blob size: [(%d, %d, %d), (%d, %d, %d)]\n', max_bbsize_blob(1, :), max_bbsize_blob(2, :));
minimum_blob_size_add = fedpars.blobsizeadd + max( [ ...
minimum_blob_size_add = fedpars.blobsizeadd(det_num, :) + max( [ ...
max_bbsize_blob(1, :) - max_bbox_proj_size(1, :); ...
max_bbox_proj_size(2, :) - max_bbsize_blob(2, :) ...
], [], 1);
fprintf(' Current Blob size add: (%d, %d, %d) -> minimum suggested: (%d, %d, %d)\n', ...
fedpars.blobsizeadd, minimum_blob_size_add);
fedpars.blobsizeadd(det_num, :), minimum_blob_size_add);
fprintf(' - Projecting volumes: ')
t = tic();
......@@ -205,4 +210,4 @@ for ii = 1:nbl
end
fprintf('Done in %f s\n', toc(t));
end % of function
end
7_fed2/gtFedTestLoadData_v2.m
+ 309
− 298
View file @ 70d12926
function [gv,gvb,bl,dmvol,fedpars, grain, blobondet] = gtFedTestLoadData_v2(fedpars, grain, parameters, dmvol)
% [gv,gvb,bl,dmvol,fedpars] = gtFedTestLoadData_v2(fedpars, grain, parameters, dmvol)
function [grain, gv, dmvol, fedpars, parameters] = gtFedTestLoadData_v2(fedpars, grain, parameters, dmvol, intvol)
% [grain, gv, dmvol, fedpars, parameters] = gtFedTestLoadData_v2(fedpars, grain, parameters, dmvol, intvol)
%
% Sets up the framework, variables, random strain distribution for simulated data.
% In the grain volume everything is defined in voxels, no dimensioning with pixel size.
% In the projections and diff. blobs, everything is defined in detector pixels (here comes the pixel size and setup geometry into play).
% In the projections and diff. blobs, everything is defined in detector pixels
% (here comes the pixel size and setup geometry into play).
% Thus the dimensioning/scaling happens in the forward projection
% functions and their derivative functions, since they relate the
% grain voxels to the detector pixels.
......@@ -18,78 +19,57 @@ function [gv,gvb,bl,dmvol,fedpars, grain, blobondet] = gtFedTestLoadData_v2(fedp
% loadbl - blobs ID-s in grain.allblobs to be generated
% dcomps - for rotations only use [1 1 1 0 0 0 0 0 0]
% blobsizeadd - blob volume padding
% vox000sam - (X,Y,Z) coordinates of voxel (0,0,0) in the Sample reference frame
% vox000sam - (X, Y, Z) coordinates of voxel (0, 0, 0) in the Sample reference frame
% grain - grain data as output from Indexter
% parameters - parameters as in DCT paramaters file
% dmvol - all deformation components per voxel for the entire volume
% (nx,ny,nz,ndef); if empty, random distribution is applied
% (nx, ny, nz, ndef); if empty, random distribution is applied
% with max limits as in fedpars.dapplim and max gradients
% as in fedpars.dappgradlim
% intvol - Volume instensities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Load and update geometry parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
labgeo = parameters.labgeo;
samgeo = parameters.samgeo;
labgeo.omstep = gtGetOmegaStepDeg(parameters);
labgeo.detnorm = gtMathsCross(labgeo.detdiru, labgeo.detdirv); % unit vector
% Detector position
% originally in detector pixels !!!
% change to grain voxels !!!
% pixel sizes in mm-s !!!
labgeo.detrefuv = [labgeo.detrefu; labgeo.detrefv];
% Rotation matrix components:
labgeo.rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
% mm/pixelsize
labgeo.detscaleu = 1 / parameters.labgeo.pixelsizeu ;
labgeo.detscalev = 1 / parameters.labgeo.pixelsizev ;
% Rotation tensor components
labgeo.Qdet = gtFedDetectorProjectionTensor( ...
labgeo.detdiru, labgeo.detdirv, labgeo.detscaleu, labgeo.detscalev);
% Update parameters
samgeo = parameters.samgeo;
recgeo = parameters.recgeo;
if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
[parameters.detgeo, labgeo] = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
end
detgeo = parameters.detgeo;
% Update labgeo
parameters.labgeo = labgeo;
% Compute detector origin from the reference position
% detorig = labgeo.detrefpos - ...
% labgeo.detrefu*labgeo.detdiru*labgeo.pixelsizeu - ...
% labgeo.detrefv*labgeo.detdirv*labgeo.pixelsizev;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Load grain parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
grsize = fedpars.grainsize; % x,y,z number of voxels in grain
gvsize = fedpars.gvsize; % x,y,z number of voxels in an element
grsize = fedpars.grainsize; % x, y, z number of voxels in grain
gvsize = fedpars.gvsize; % x, y, z number of voxels in an element
% Grain volume
grvol = gtFedTestGenerateGrainVolume(grsize);
grvol = gtFedTestGenerateGrainVolume(grsize);
fedpars.grainvol = grvol;
% Determine number of volume elements along X,Y,Z, based on grain volume
% Determine number of volume elements along X, Y, Z, based on grain volume
% and element size:
ngv(1) = ceil(size(grvol,1)/gvsize(1));
ngv(2) = ceil(size(grvol,2)/gvsize(2));
ngv(3) = ceil(size(grvol,3)/gvsize(3));
[grenv(1), grenv(2), grenv(3)] = size(grvol);
ngv = ceil(grenv ./ gvsize);
fedpars.ngv = ngv;
% Size of grain envelope
grenv(1) = ngv(1)*gvsize(1);
grenv(2) = ngv(2)*gvsize(2);
grenv(3) = ngv(3)*gvsize(3);
grenv = ngv .* gvsize;
fedpars.grainenv = grenv;
% Padded grain volume
grvol_pad = false(grenv);
grvol_pad(1:size(grvol,1),1:size(grvol,2),1:size(grvol,3)) = grvol;
grvol_pad(1:size(grvol, 1), 1:size(grvol, 2), 1:size(grvol, 3)) = grvol;
grvol = grvol_pad;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -109,343 +89,374 @@ nd = sum(fedpars.dcomps);
gv.ind = gv.ind';
% Elements' positions in the sample volume in mm!
gv.cs = fedpars.vox000sam(ones(nv,1), :)' + gv.ind .* parameters.recgeo.voxsize(ones(nv,1), :)';
%gv.cs1 = gv.cs;
gv.cs = fedpars.vox000sam(ones(nv, 1), :)' + gv.ind .* recgeo(1).voxsize(ones(nv, 1), :)';
% Grain center position in sample in mm!
%gc = mean(double(gv.cs),2);
gc = fedpars.vox000sam + (grsize/2 + [0.5, 0.5, 0.5]) .* parameters.recgeo.voxsize;
gc = fedpars.vox000sam + (grsize/2 + [0.5, 0.5, 0.5]) .* recgeo(1).voxsize;
fedpars.graincentsam = gc;
grain.center = gc;
% Element diffracting power
gv.pow = ones(1, nv);
if (~exist('intvol', 'var') || isempty(intvol))
gv.pow = ones(1, nv, fedpars.gvtype);
else
if (numel(intvol) ~= nv)
error('gtFedTestLoadData_v2:wrong_argument', ...
'The volume of voxel intensities should have the same size as dmvol''s first 3 coordinates')
end
gv.pow = reshape(intvol, 1, []);
if (strcmp(fedpars.gvtype, 'single'))
gv.pow = single(gv.pow);
end
end
% Element actual strain state
gv.d = zeros(nd, nv);
gv.d = zeros(nd, nv, fedpars.gvtype);
% Element change in strain state
gv.dd = zeros(nd, nv);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Computing allblobs info
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
grain = gtCalculateGrain(grain, parameters);
gv.dd = zeros(nd, nv, fedpars.gvtype);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Iinitial blob parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Apply strain distribution to grain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Generating blob data...')
disp('Applying random deformation...')
tic
% Blobs to load and use
lbl = fedpars.loadbl;
% Number of blobs to consider
nbl = length(fedpars.loadbl(:));
% Create a blank blob structure
tmp = gtFedBlobDefinition();
bl(1:nbl) = tmp;
% Loop through blobs
blobsdiffaway = [];
blobondet = [];
for ii = 1:nbl
% ID-s
bl(ii).blobIDfed = ii;
bl(ii).blobIDgrain = fedpars.loadbl(ii);
% Load diffraction paramaters of blobs
% Plane families
bl(ii).hkl = grain.allblobs.hkl(lbl(ii), :);
bl(ii).hklsp = grain.allblobs.hklsp(lbl(ii), :);
bl(ii).thetatype = grain.allblobs.thetatype(lbl(ii));
% Initial strained signed plane normal in SAMPLE reference
bl(ii).plsam = grain.allblobs.pl(lbl(ii), :);
disp('Active deformation components:')
disp(fedpars.dcomps)
% Initial strained non-signed plane normal in SAMPLE reference
bl(ii).plorig = grain.allblobs.plorig(lbl(ii), :);
% Measured (true, real) strain distribution:
if (~exist('dmvol', 'var') || isempty(dmvol))
[gvdm, dmvol, gradok] = gtFedTestGenerateRandomDef(fedpars, nd);
else
[gvdm, dmvol_size] = gtDefDmvol2Gvdm(dmvol);
if (dmvol_size(4) == 3)
mean_gvdm = mean(gvdm, 2);
fprintf('Recentering GVDM and DMVOL to the mean orientation: (%f, %f, %f), from: (%f, %f, %f).\n', ...
grain.R_vector, mean_gvdm)
gvdm_shift = grain.R_vector' - mean_gvdm;
gvdm = gvdm + gvdm_shift(:, ones(size(gvdm, 2), 1));
dmvol = gtDefGvdm2Dmvol(gvdm, dmvol_size);
end
gradok = [];
end
% Initial strained signed plane normal in the diffraction position in
% LAB reference
bl(ii).pl0 = grain.allblobs.pllab(lbl(ii), :);
gv.dm = gvdm;
gv.d = gv.dm;
if (strcmp(fedpars.gvtype, 'single'))
disp('Converting type to single...')
fn = fieldnames(gv);
for ii = 1:length(fn)
gv.(fn{ii}) = single(gv.(fn{ii}));
end
end
% Initial strained diffraction angles and rotation matrix
bl(ii).th0 = grain.allblobs.theta(lbl(ii)); % initial theta
bl(ii).sinth0 = grain.allblobs.sintheta(lbl(ii)); % initial sin(theta)
bl(ii).eta0 = grain.allblobs.eta(lbl(ii)); % initial eta
bl(ii).om0 = grain.allblobs.omega(lbl(ii)); % initial omega
bl(ii).omind = grain.allblobs.omind(lbl(ii)); % initial omega index (1..4)
bl(ii).S0 = grain.allblobs.srot(:, :, lbl(ii)); % initial rotation tensor
% Initial strained diffraction vector in LAB reference
bl(ii).dvec0 = grain.allblobs.dvec(lbl(ii), :);
if ~all(gradok)
disp('Gradients are higher than limit for components:')
disp(find(~gradok)')
end
% Initial strained diffraction vector in SAMPLE reference
bl(ii).dvecsam = grain.allblobs.dvecsam(lbl(ii), :);
% figure('name', 'Epsilon33')
% slice(dmvol(:, :, :, 3), grsize(1)/2, grsize(2)/2, grsize(3)/2);
% colorbar
% Initial centroid in image
bl(ii).u0im = grain.allblobs.uvw(lbl(ii), :);
disp('Largest deformation component:')
disp(max(dmvol(:)))
disp('Smallest deformation component:')
disp(min(dmvol(:)))
if isempty(bl(ii).u0im)
blobsdiffaway = [blobsdiffaway, ii];
else
% Padding around blob volume
bl(ii).addbl = fedpars.blobsizeadd;
toc
% Check if blob hits the detector area
uok = (bl(ii).u0im(1) < parameters.acq.xdet) & (bl(ii).u0im(1) > 1);
vok = (bl(ii).u0im(2) < parameters.acq.ydet) & (bl(ii).u0im(2) > 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Computing allblobs info
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (uok && vok)
blobondet = [blobondet, ii];
end
end
% Cleaning grain struct
if (isfield(grain, 'id'))
grain = struct('id', grain.id, 'phaseid', grain.phaseid, ...
'center', grain.center, 'R_vector', grain.R_vector);
else
grain = struct('id', 1, 'phaseid', grain.phaseid, ...
'center', grain.center, 'R_vector', grain.R_vector);
end
disp('Blobs inside the detector area are:')
disp(blobondet)
grain = gtCalculateGrain(grain, parameters);
if ~isempty(blobsdiffaway)
disp('The following projections (blobs) are directed away from the detector plane.')
disp(blobsdiffaway)
error('Remove projections (blobs) that do not hit the detector.')
end
% Blobs to load and use
lbl = fedpars.loadbl;
% Number of blobs to consider
nbl = numel(lbl(:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Arrows of projection vectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quiver figure of projection vectors
dvecsam = vertcat(bl(:).dvecsam);
dvecsam = grain.allblobs.dvecsam(lbl, :);
figure('name','Projection vectors')
quiver3(zeros(nbl,1), zeros(nbl,1), zeros(nbl,1), ...
dvecsam(:,1), dvecsam(:,2), dvecsam(:,3));
figure('name', 'Projection vectors')
quiver3(zeros(nbl, 1), zeros(nbl, 1), zeros(nbl, 1), ...
dvecsam(:, 1), dvecsam(:, 2), dvecsam(:, 3));
axis equal
hold on
for ii = 1:length(bl)
text(dvecsam(ii,1), dvecsam(ii,2), dvecsam(ii,3), ...
num2str(ii), 'VerticalAlignment', 'Bottom', 'Color','m')
text(dvecsam(ii,1), dvecsam(ii,2), dvecsam(ii,3), ...
num2str(bl(ii).blobIDgrain), 'VerticalAlignment','Top','Color','r')
for ii = 1:nbl
text(dvecsam(ii, 1), dvecsam(ii, 2), dvecsam(ii, 3), ...
num2str(ii), 'VerticalAlignment', 'Bottom', 'Color', 'm')
text(dvecsam(ii, 1), dvecsam(ii, 2), dvecsam(ii, 3), ...
num2str(fedpars.loadbl(ii)), 'VerticalAlignment', 'Top', 'Color', 'r')
end
title('blob ID in FED (magenta), blob ID in grain (red)')
drawnow();
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Preallocate blob volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Iinitial blob parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxdiag = ceil(sqrt(grenv*grenv'));
num_dets = numel(detgeo);
grain.proj(1:num_dets) = struct( ...
'stack', [], 'num_cols', [], 'num_rows', [], 'bl', [], ...
'vol_size_x', [], 'vol_size_y', [], 'vol_size_z', [], 'shift', [], ...
'ondet', [], 'included', [], 'selected', [], 'centerpix', [], ...
'gvb', []);
for n = 1:num_dets
ondet = find(grain.allblobs.detector(n).ondet);
nbl = numel(ondet);
is_nearfield = norm(detgeo(n).detrefpos) < 30;
% If it is farfield, we should adjust the voxel size to be 1/4 of the
% pixel size
if (~is_nearfield)
voxsize = mean([detgeo(n).pixelsizeu, detgeo(n).pixelsizev]) / 4;
recgeo(n).voxsize = voxsize([1 1 1]);
parameters.recgeo(n) = recgeo(n);
end
fprintf( ['Generating blob data (for detector: %d, blobs on the' ...
' detector: %d):\n - Allocating..'], n, nbl);
c = tic();
for ii = 1:nbl
% Bounding box sizes and boundaries for processing
bl(ii).bbsize = [maxdiag + 2*bl(ii).addbl(1), ...
maxdiag + 2*bl(ii).addbl(2), ...
2 + 2*bl(ii).addbl(3) ];
tmp = round( bl(ii).u0im(1) - bl(ii).bbsize(1)/2 );
bl(ii).bbuim = [tmp, tmp + bl(ii).bbsize(1) - 1];
tmp = round( bl(ii).u0im(2) - bl(ii).bbsize(2)/2 );
bl(ii).bbvim = [tmp, tmp + bl(ii).bbsize(2) - 1];
tmp = round( bl(ii).u0im(3) - bl(ii).bbsize(3)/2 );
bl(ii).bbwim = [tmp, tmp + bl(ii).bbsize(3) - 1];
% Blob origin and center [u,v,w]:
bl(ii).bbor = [bl(ii).bbuim(1)-1, bl(ii).bbvim(1)-1, bl(ii).bbwim(1)-1];
bl(ii).bbcent = [round( (bl(ii).bbuim(1) + bl(ii).bbuim(2))/2 ),...
round( (bl(ii).bbvim(1) + bl(ii).bbvim(2))/2 ),...
round( (bl(ii).bbwim(1) + bl(ii).bbwim(2))/2 ) ];
% Projection parameters
bl(ii).proj_geom = gtGeoProjForReconstruction( ...
bl(ii).dvecsam, bl(ii).om0, gc, ...
[bl(ii).bbuim(1) bl(ii).bbvim(1) bl(ii).bbsize(1:2)], ...
[], parameters.labgeo, parameters.samgeo, parameters.recgeo, ...
'ASTRA_grain');
%%% Needed for FED
% bl(ii).int = zeros(bl(ii).bbsize, fedpars.bltype);
% bl(ii).intd = zeros(bl(ii).bbsize, fedpars.bltype);
bl(ii).intm = zeros(bl(ii).bbsize, fedpars.bltype);
bl(ii).mask = ones(bl(ii).bbsize, fedpars.bltype);
% u in blob = u in image minus u at blob origin
bl(ii).u0bl = bl(ii).u0im - bl(ii).bbor;
end
% Create a blank blob structure
bl = repmat(gtFedBlobDefinition(), nbl, 1);
gvb = [];
toc
included = (1:nbl)';
selected = true(size(included));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Initial element parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii = 1:nbl
% Load diffraction paramaters of blobs
ii_b = ondet(included(ii));
disp('Generating element data ...')
tic
% Plane families
bl(ii).hkl = grain.allblobs.hkl(ii_b, :);
bl(ii).hklsp = grain.allblobs.hklsp(ii_b, :);
bl(ii).thetatype = grain.allblobs.thetatype(ii_b);
gvcs = gv.cs;
% Initial strained signed plane normal in SAMPLE reference
bl(ii).plsam = grain.allblobs.pl(ii_b, :);
for ii = nbl:-1:1
gvb(ii).uc0bl = zeros(size(gvcs));
gvb(ii).uc = zeros(size(gvcs));
gvb(ii).ucd = zeros(size(gvcs));
gvb(ii).ucblmod = zeros(size(gvcs));
gvb(ii).ucbl8 = zeros(1, size(gvcs,2)*8);
gvb(ii).ucbl8ini= zeros(1, size(gvcs,2)*8);
gvb(ii).ucbl8ok = 0;
% Initial strained non-signed plane normal in SAMPLE reference
bl(ii).plorig = grain.allblobs.plorig(ii_b, :);
% Expand rotation tensors
subs = {(1:3)', (1:3)', ones(1,nv)};
rot_l2s = bl(ii).S0(subs{:});
% Initial strained signed plane normal in the diffraction position in
% LAB reference
bl(ii).pl0 = grain.allblobs.pllab(ii_b, :);
% Expand diffraction vector
dvec = bl(ii).dvec0';
dvec = dvec(:,ones(1,nv));
% Initial strained diffraction angles and rotation matrix
bl(ii).th0 = grain.allblobs.theta(ii_b); % initial theta
bl(ii).sinth0 = grain.allblobs.sintheta(ii_b); % initial sin(theta)
bl(ii).eta0 = grain.allblobs.eta(ii_b); % initial eta
bl(ii).om0 = grain.allblobs.omega(ii_b); % initial omega
bl(ii).omind = grain.allblobs.omind(ii_b); % initial omega index (1..4)
bl(ii).S0 = grain.allblobs.srot(:, :, ii_b); % initial rotation tensor
% Expand omega
om = bl(ii).om0(1,ones(1,nv));
% Initial strained diffraction vector in LAB reference
bl(ii).dvec0 = grain.allblobs.dvec(ii_b, :);
rot_s2l = permute(rot_l2s, [2 1 3]);
gvcs_lab = gtGeoSam2Lab(gvcs', rot_s2l, labgeo, samgeo, false);
% Initial strained diffraction vector in SAMPLE reference
bl(ii).dvecsam = grain.allblobs.dvecsam(ii_b, :);
% Initial absolut u of gv center
uc0im = gtFedPredictUVWMultiple([], dvec, gvcs_lab', ...
labgeo.detrefpos', labgeo.detnorm', labgeo.Qdet, ...
labgeo.detrefuv, om, labgeo.omstep );
% Initial centroid in image
bl(ii).u0im = grain.allblobs.detector(n).uvw(ii_b, :);
% Initial u of gv centre relative to blob origin
bbor = bl(ii).bbor';
gvb(ii).uc0bl = uc0im - bbor(:,ones(1,nv));
end
% Padding around blob volume
bl(ii).addbl = fedpars.blobsizeadd(n, :);
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Preallocate blob volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Apply strain distribution to grain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxdiag = ceil(sqrt(grenv*grenv'));
disp('Applying random deformation...')
tic
for ii = 1:nbl
% Bounding box sizes and boundaries for processing
bl(ii).bbsize = [ maxdiag, maxdiag, 2 ] + 2 * bl(ii).addbl;
disp('Active deformation components:')
disp(fedpars.dcomps)
im_low_lims = round( bl(ii).u0im - bl(ii).bbsize/2 );
bl(ii).bbuim = [im_low_lims(1), im_low_lims(1) + bl(ii).bbsize(1) - 1];
bl(ii).bbvim = [im_low_lims(2), im_low_lims(2) + bl(ii).bbsize(2) - 1];
bl(ii).bbwim = [im_low_lims(3), im_low_lims(3) + bl(ii).bbsize(3) - 1];
% Measured (true,real) strain distribution:
if ~exist('dmvol','var') || isempty(dmvol)
[gvdm, dmvol, gradok] = gtFedTestGenerateRandomDef(fedpars, nd);
gv.dm = gvdm;
else
[gvdm, dmvol_size] = gtDefDmvol2Gvdm(dmvol);
if (dmvol_size(4) == 3)
mean_gvdm = mean(gvdm, 2);
fprintf('Recentering GVDM and DMVOL to the mean orientation: (%f, %f, %f), from: (%f, %f, %f).\n', ...
grain.R_vector, mean_gvdm)
gvdm_shift = grain.R_vector' - mean_gvdm;
gvdm = gvdm + gvdm_shift(:, ones(size(gvdm, 2), 1));
dmvol = gtDefGvdm2Dmvol(gvdm, dmvol_size);
% Blob origin and center [u, v, w]:
bl(ii).bbor = [bl(ii).bbuim(1)-1, bl(ii).bbvim(1)-1, bl(ii).bbwim(1)-1];
bl(ii).bbcent = [round( (bl(ii).bbuim(1) + bl(ii).bbuim(2))/2 ), ...
round( (bl(ii).bbvim(1) + bl(ii).bbvim(2))/2 ), ...
round( (bl(ii).bbwim(1) + bl(ii).bbwim(2))/2 ) ];
% Projection parameters
bl(ii).proj_geom = gtGeoProjForReconstruction( ...
bl(ii).dvecsam, bl(ii).om0, gc, ...
[bl(ii).bbuim(1) bl(ii).bbvim(1) bl(ii).bbsize(1:2)], ...
[], detgeo(n), labgeo, samgeo, recgeo(n), ...
'ASTRA_grain');
bl(ii).intm = zeros(bl(ii).bbsize, fedpars.bltype);
bl(ii).mask = ones(bl(ii).bbsize, fedpars.bltype);
% u in blob = u in image minus u at blob origin
bl(ii).u0bl = bl(ii).u0im - bl(ii).bbor;
end
gv.dm = gvdm;
gradok = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Initial element parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~all(gradok)
disp('Gradients are higher than limit for components:')
disp(find(~gradok)')
end
fprintf('\b\b(%f seconds).\n - Generating element data..', toc(c))
c = tic();
% figure('name','Epsilon33')
% slice(dmvol(:,:,:,3),grsize(1)/2,grsize(2)/2,grsize(3)/2);
% colorbar
gvcs = gv.cs;
disp('Largest deformation component:')
disp(max(dmvol(:)))
disp('Smallest deformation component:')
disp(min(dmvol(:)))
for ii = nbl:-1:1
% gvb(ii).uc0bl = zeros(size(gvcs), fedpars.gvtype);
gvb(ii).uc = zeros(size(gvcs), fedpars.gvtype);
gvb(ii).ucd = zeros(size(gvcs), fedpars.gvtype);
gvb(ii).ucblmod = zeros(size(gvcs), fedpars.gvtype);
gvb(ii).ucbl8 = zeros(1, size(gvcs, 2)*8, fedpars.gvtype);
gvb(ii).ucbl8ini= zeros(1, size(gvcs, 2)*8, fedpars.gvtype);
gvb(ii).ucbl8ok = 0;
toc
% Expand diffraction vector
dvec = bl(ii).dvec0';
dvec = dvec(:, ones(1, nv));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Simulate "measured" diffraction patterns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Expand omega
om = bl(ii).om0(1, ones(1, nv));
disp('Simulating true ("measured") blob intensities...')
gv.d = gv.dm;
[gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters);
% Change from sample to lab coordinates
gvcs_lab = gtGeoSam2Lab(gvcs', bl(ii).S0', labgeo, samgeo, false);
% Set back original zero actual strain
for ii = 1:nbl
bl(ii).intm = bl(ii).int;
bl(ii).int(:) = 0;
bl(ii).intd(:) = bl(ii).int - bl(ii).intm;
% Initial absolute u of gv center
uc0im = gtFedPredictUVWMultiple([], dvec, gvcs_lab', ...
detgeo(n).detrefpos', detgeo(n).detnorm', detgeo(n).Qdet, ...
[detgeo(n).detrefu, detgeo(n).detrefv]', om, labgeo.omstep );
bl(ii).comim = NaN(1,3);
bl(ii).combl = NaN(1,3);
% Initial u of gv centre relative to blob origin
bbor = bl(ii).bbor';
gvb(ii).uc0bl = cast(uc0im - bbor(:, ones(1, nv)), fedpars.gvtype);
end
gvb(ii).ucbl8ini = gvb(ii).ucbl8;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Simulate "measured" diffraction patterns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Calculate boundaries of measured blobs inside blob bounding box...')
for ii = 1:nbl
uproj = sum(sum(abs(bl(ii).intm), 2), 3);
bl(ii).mbbu = [find(uproj,1,'first'), find(uproj,1,'last')];
fprintf('\b\b(%f seconds)\n - Simulating true ("measured") blob intensities..', toc(c))
[gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, n);
vproj = sum(sum(abs(bl(ii).intm), 1), 3);
bl(ii).mbbv = [find(vproj,1,'first'), find(vproj,1,'last')];
% Set back original zero actual strain
for ii = 1:nbl
bl(ii).intm = bl(ii).int;
bl(ii).int(:) = 0;
bl(ii).intd(:) = bl(ii).int - bl(ii).intm;
wproj = sum(sum(abs(bl(ii).intm), 1), 2);
bl(ii).mbbw = [find(wproj,1,'first'), find(wproj,1,'last')];
bl(ii).comim = NaN(1, 3);
bl(ii).combl = NaN(1, 3);
bl(ii).mbbsize = [bl(ii).mbbu(2) - bl(ii).mbbu(1) + 1, ...
bl(ii).mbbv(2) - bl(ii).mbbv(1) + 1, ...
bl(ii).mbbw(2) - bl(ii).mbbw(1) + 1 ];
end
gvb(ii).ucbl8ini = gvb(ii).ucbl8;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Angle between consecutive blobs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\b\b(%f seconds)\n - Calculate boundaries of measured blobs inside blob bounding box..', toc(c))
for ii = 1:nbl
uproj = sum(sum(abs(bl(ii).intm), 2), 3);
bl(ii).mbbu = [find(uproj, 1, 'first'), find(uproj, 1, 'last')];
for ii = 1:nbl
nii = mod(ii, nbl) + 1;
vproj = sum(sum(abs(bl(ii).intm), 1), 3);
bl(ii).mbbv = [find(vproj, 1, 'first'), find(vproj, 1, 'last')];
fedpars.blobnextang(ii) = abs(bl(ii).plsam*bl(nii).plsam');
fedpars.blobnextang(ii) = min(fedpars.blobnextang(ii), 1);
fedpars.blobnextang(ii) = acosd(fedpars.blobnextang(ii));
wproj = sum(sum(abs(bl(ii).intm), 1), 2);
bl(ii).mbbw = [find(wproj, 1, 'first'), find(wproj, 1, 'last')];
fedpars.blobnextdet(ii) = abs(det(bl(ii).U0*bl(nii).U0'));
end
bl(ii).mbbsize = [bl(ii).mbbu(2) - bl(ii).mbbu(1) + 1, ...
bl(ii).mbbv(2) - bl(ii).mbbv(1) + 1, ...
bl(ii).mbbw(2) - bl(ii).mbbw(1) + 1 ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transform into single
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Angle between consecutive blobs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Transform into single precision
if strcmp(fedpars.gvtype, 'single')
disp('Converting type to single...')
fn = fieldnames(gv);
for ii = 1:length(fn)
gv.(fn{ii}) = single(gv.(fn{ii}));
for ii = 1:nbl
nii = mod(ii, nbl) + 1;
fedpars.blobnextang(ii) = abs(bl(ii).plsam*bl(nii).plsam');
fedpars.blobnextang(ii) = min(fedpars.blobnextang(ii), 1);
fedpars.blobnextang(ii) = acosd(fedpars.blobnextang(ii));
fedpars.blobnextdet(ii) = abs(det(bl(ii).U0*bl(nii).U0'));
end
fn = fieldnames(gvb);
for ii = 1:length(fn)
for jj = 1:length(gvb)
gvb(jj).(fn{ii}) = single(gvb(jj).(fn{ii}));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transform into single
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Transform into single precision
if strcmp(fedpars.gvtype, 'single')
fprintf('\b\b(%f seconds)\n - Converting type to single..', toc(c))
c = tic();
fn = fieldnames(gvb);
for ii = 1:length(fn)
for jj = 1:length(gvb)
gvb(jj).(fn{ii}) = single(gvb(jj).(fn{ii}));
end
end
end
end
disp('True blob dimensions:')
disp(vertcat(bl(:).mbbsize))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Assembling the final proj structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\b\b(%f seconds)\n - Producing..', toc(c))
c = tic();
end % of main function
stack = arrayfun(@(x){sum(x.intm, 3)}, bl);
stack = cat(3, stack{:});
vol_size = gtGeoSam2Sam(grsize, recgeo(1), recgeo(n), 1);
vol_size = ceil(vol_size * 1.2);
if (is_nearfield)
vol_size(3) = max(vol_size(3), 8);
end
grain.proj(n).stack = permute(stack, [2 3 1]);
grain.proj(n).num_cols = size(grain.proj(n).stack, 1);
grain.proj(n).num_rows = size(grain.proj(n).stack, 3);
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;
grain.proj(n).gvb = gvb;
fprintf('\b\b(%f seconds)\n - True blob dimensions:\n', toc(c))
disp(vertcat(bl(:).mbbsize))
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