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

gtFedTestLoadData: fixed previous commit to use less memory (in exchange for being more verbose) 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 5c250fee
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
7_fed2/gtFedFwdProjExact.m
+ 190
− 194
View file @ e5cb15c0
......@@ -8,6 +8,7 @@ fprintf('Forward projection (%s):\n', mfilename)
nbl = length(bl);
uinds = gv.used_ind;
nv = numel(uinds);
detgeo = parameters.detgeo(det_num);
labgeo = parameters.labgeo;
......@@ -19,19 +20,11 @@ defmethod = fedpars.defmethod;
gvd = gv.d(:, uinds);
gvpow = gv.pow(uinds);
% Deformation tensor (relative to reference state) from its components
defT = gtFedDefTensorFromComps(gvd, dcomps, defmethod, 0);
num_oversampling = size(gv.cs, 3);
if (num_oversampling == 1)
gvcs = gv.cs(:, uinds);
else
gvcs = reshape(gv.cs(:, uinds, :), size(gv.cs, 1), []);
gvpow = repmat(gvpow, [num_oversampling, 1]) / num_oversampling;
gvpow = gvpow / num_oversampling;
defT = repmat(defT, [1, 1, num_oversampling]);
end
nv = numel(uinds) * (num_oversampling ^ 3);
% Deformation tensor (relative to reference state) from its components
defT = gtFedDefTensorFromComps(gvd, dcomps, defmethod, 0);
plorig = vertcat(bl(:).plorig);
sinth0 = vertcat(bl(:).sinth0);
......@@ -39,198 +32,201 @@ omind = vertcat(bl(:).omind);
bbor = vertcat(bl(:).bbor);
bbsize = vertcat(bl(:).bbsize);
uvw = cell(1, nbl);
fprintf(' - Computing indices and bbsizes: ')
t = tic();
for ii = 1:nbl
num_chars = fprintf('%03d/%03d', ii, nbl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate new detector coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% New deformed plane normals and relative elongations (relative to
% reference state)
plsam = plorig(ii, :)'; % ! use plorig and not plsam to keep omega order below!
plsam = plsam(:, ones(1, nv));
[pl_samd, drel] = gtStrainPlaneNormals(plsam, defT); % unit column vectors
% New sin(theta)
sinth = sinth0(ii) ./ drel;
% The plane normals need to be brought in the Lab reference where the
% beam direction and rotation axis are defined.
% Use the Sample -> Lab orientation transformation assuming omega=0;
% (vector length preserved for free vectors)
pl_labd = gtGeoSam2Lab(pl_samd', eye(3), labgeo, samgeo, true, false)';
% Predict omega angles: 4 for each plane normal
[om, pllab, ~, rot_l2s] = gtFedPredictOmegaMultiple(pl_labd, ...
sinth, labgeo.beamdir', labgeo.rotdir', labgeo.rotcomp, omind(ii));
% Delete those where no reflection occurs
if any(isnan(om));
error('No diffraction from some elements after deformation.')
for ii_ss = 1:num_oversampling
gvcs = gv.cs(:, uinds, ii_ss);
uvw = cell(1, nbl);
fprintf(' - Computing indices and bbsizes: ')
t = tic();
for ii = 1:nbl
num_chars = fprintf('%03d/%03d', ii, nbl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate new detector coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% New deformed plane normals and relative elongations (relative to
% reference state)
plsam = plorig(ii, :)'; % ! use plorig and not plsam to keep omega order below!
plsam = plsam(:, ones(1, nv));
[pl_samd, drel] = gtStrainPlaneNormals(plsam, defT); % unit column vectors
% New sin(theta)
sinth = sinth0(ii) ./ drel;
% The plane normals need to be brought in the Lab reference where the
% beam direction and rotation axis are defined.
% Use the Sample -> Lab orientation transformation assuming omega=0;
% (vector length preserved for free vectors)
pl_labd = gtGeoSam2Lab(pl_samd', eye(3), labgeo, samgeo, true, false)';
% Predict omega angles: 4 for each plane normal
[om, pllab, ~, rot_l2s] = gtFedPredictOmegaMultiple(pl_labd, ...
sinth, labgeo.beamdir', labgeo.rotdir', labgeo.rotcomp, omind(ii));
% Delete those where no reflection occurs
if any(isnan(om));
error('No diffraction from some elements after deformation.')
end
% Diffraction vector
dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
rot_s2l = permute(rot_l2s, [2 1 3]);
gvcs_lab = gtGeoSam2Lab(gvcs', rot_s2l, labgeo, samgeo, false);
ucbl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
[detgeo.detrefu, detgeo.detrefv]', om, labgeo.omstep);
% Detector coordinates U,V in blob
bbort = bbor(ii, :)';
ucbl = ucbl - bbort(:, ones(1, nv));
% U,V shift relative to reference state
gvb(ii).ucd = ucbl - gvb(ii).uc0bl - gvb(ii).uc;
gvb(ii).uc = ucbl - gvb(ii).uc0bl;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sum intensities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Distribute value over 8 pixels
% Modulus of uvw position
gvb(ii).ucblmod = mod(ucbl, 1);
% This way of rounding is correct also for round values when mu=mod...
% is used above:
ul = floor(ucbl); % lower pixel indices
uh = ul + 1; % higher pixel indices
% uvw indices of 8 nearest neighbours (3xnx8)
ucbl8 = [ ul, ...
[ul([1,2],:); uh(3,:)], ...
[ul(1,:); uh(2,:); ul(3,:)], ...
[ul(1,:); uh([2,3],:)], ...
[uh(1,:); ul([2,3],:)], ...
[uh(1,:); ul(2,:); uh(3,:)], ...
[uh([1,2],:); ul(3,:)], ...
uh ];
% uvw indices of 8 nearest neighbours (nx3)
uvw{ii} = ucbl8';
% Don't use sub2ind, it's slow
gvb(ii).ucbl8 = ucbl8(1, :) ...
+ bbsize(ii, 1) .* (ucbl8(2, :)-1) ...
+ bbsize(ii, 1) .* bbsize(ii, 2) .* (ucbl8(3, :)-1);
fprintf(repmat('\b', [1, num_chars]));
end
fprintf('Done in %f s\n', toc(t));
% Diffraction vector
dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
rot_s2l = permute(rot_l2s, [2 1 3]);
gvcs_lab = gtGeoSam2Lab(gvcs', rot_s2l, labgeo, samgeo, false);
ucbl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
[detgeo.detrefu, detgeo.detrefv]', om, labgeo.omstep);
% Detector coordinates U,V in blob
bbort = bbor(ii, :)';
ucbl = ucbl - bbort(:, ones(1, nv));
% U,V shift relative to reference state
gvb(ii).ucd = ucbl - gvb(ii).uc0bl - gvb(ii).uc;
gvb(ii).uc = ucbl - gvb(ii).uc0bl;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sum intensities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Distribute value over 8 pixels
% Modulus of uvw position
gvb(ii).ucblmod = mod(ucbl, 1);
% This way of rounding is correct also for round values when mu=mod...
% is used above:
ul = floor(ucbl); % lower pixel indices
uh = ul + 1; % higher pixel indices
% uvw indices of 8 nearest neighbours (3xnx8)
ucbl8 = [ ul, ...
[ul([1,2],:); uh(3,:)], ...
[ul(1,:); uh(2,:); ul(3,:)], ...
[ul(1,:); uh([2,3],:)], ...
[uh(1,:); ul([2,3],:)], ...
[uh(1,:); ul(2,:); uh(3,:)], ...
[uh([1,2],:); ul(3,:)], ...
uh ];
% uvw indices of 8 nearest neighbours (nx3)
uvw{ii} = ucbl8';
% Don't use sub2ind, it's slow
gvb(ii).ucbl8 = ucbl8(1, :) ...
+ bbsize(ii, 1) .* (ucbl8(2, :)-1) ...
+ bbsize(ii, 1) .* bbsize(ii, 2) .* (ucbl8(3, :)-1);
fprintf(repmat('\b', [1, num_chars]));
end
fprintf('Done in %f s\n', toc(t));
fprintf(' - Computing max BBox size and feasibilities:\n')
max_bbox_proj_size = [inf, inf, inf; 0, 0, 0];
max_bbsize_blob = [1, 1, 1; 0, 0, 0];
for ii = 1:nbl
bbsize_blob = bbsize(ii,:);
min_ind_blob = min(uvw{ii}, [], 1);
max_ind_blob = max(uvw{ii}, [], 1);
fprintf(' - Computing max BBox size and feasibilities:\n')
max_bbox_proj_size = [inf, inf, inf; 0, 0, 0];
max_bbsize_blob = [1, 1, 1; 0, 0, 0];
for ii = 1:nbl
bbsize_blob = bbsize(ii,:);
min_ind_blob = min(uvw{ii}, [], 1);
max_ind_blob = max(uvw{ii}, [], 1);
max_bbsize_blob(2, :) = max(max_bbsize_blob(2, :), bbsize_blob);
max_bbsize_blob(2, :) = max(max_bbsize_blob(2, :), bbsize_blob);
max_bbox_proj_size(1, :) = min(max_bbox_proj_size(1, :), min_ind_blob);
max_bbox_proj_size(2, :) = max(max_bbox_proj_size(2, :), max_ind_blob);
max_bbox_proj_size(1, :) = min(max_bbox_proj_size(1, :), min_ind_blob);
max_bbox_proj_size(2, :) = max(max_bbox_proj_size(2, :), max_ind_blob);
if (any(min_ind_blob < 1) || any(max_ind_blob > bbsize_blob))
fprintf(' + Blob: %d, Blobsize: (%d, %d, %d), bbox projection [(%d, %d, %d), (%d, %d, %d)]\n', ...
ii, bbsize_blob, min_ind_blob, max_ind_blob);
if (any(min_ind_blob < 1) || any(max_ind_blob > bbsize_blob))
fprintf(' + Blob: %d, Blobsize: (%d, %d, %d), bbox projection [(%d, %d, %d), (%d, %d, %d)]\n', ...
ii, bbsize_blob, min_ind_blob, max_ind_blob);
end
end
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, :));
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, :));
if (isfield(fedpars, 'detector'))
blob_size_add = fedpars.detector(det_num).blobsizeadd;
else
blob_size_add = fedpars.blobsizeadd(det_num, :);
end
minimum_blob_size_add = blob_size_add + 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', ...
blob_size_add(det_num, :), minimum_blob_size_add);
fprintf(' - Projecting volumes: ')
t = tic();
for ii = 1:nbl
num_chars = fprintf('%03d/%03d', ii, nbl);
mu = gvb(ii).ucblmod;
mu1 = mu(1, :);
mu2 = mu(2, :);
mu3 = mu(3, :);
mu12 = mu1 .* mu2;
mu13 = mu1 .* mu3;
mu23 = mu2 .* mu3;
mu123 = mu12 .* mu3;
ints = zeros(1, numel(gvpow) * 8, 'like', bl(ii).int);
ints( 1 : nv) = gvpow.*(-mu1 - mu2 - mu3 + mu12 + mu13 + mu23 - mu123 + 1);
ints( nv+1 : 2*nv) = gvpow.*(mu3 - mu13 - mu23 + mu123);
ints(2*nv+1 : 3*nv) = gvpow.*(mu2 - mu12 - mu23 + mu123);
ints(3*nv+1 : 4*nv) = gvpow.*(mu23 - mu123);
ints(4*nv+1 : 5*nv) = gvpow.*(mu1 - mu12 - mu13 + mu123);
ints(5*nv+1 : 6*nv) = gvpow.*(mu13 - mu123);
ints(6*nv+1 : 7*nv) = gvpow.*(mu12 - mu123);
ints(7*nv+1 : 8*nv) = gvpow.*mu123;
% Accumulate all intensities in the blob voxel-wise
% 'uvw' needs to be nx3; 'ints' is now 1xnx8
try
bl(ii).int = accumarray(uvw{ii}, ints, bbsize(ii,:));
catch err
disp(' ')
disp('ERROR')
disp(['The error is probably caused by the projected intensities' ...
' falling outside the blob volume. Try increasing the blob' ...
' volume padding: fedpars.detector(det_ind).blobsizeadd'])
disp('Blob ID:')
disp(ii)
disp('Blob size:')
disp(bbsize(ii,:))
disp('Min projected U,V,W coordinates:')
disp(min(uvw{ii},[],1))
disp('Max projected U,V,W coordinates:')
disp(max(uvw{ii},[],1))
rethrow(err)
if (isfield(fedpars, 'detector'))
blob_size_add = fedpars.detector(det_num).blobsizeadd;
else
blob_size_add = fedpars.blobsizeadd(det_num, :);
end
% Smooth blob
if ~strcmp(fedpars.intsmoothalg, 'none')
bl(ii).int = smooth3(bl(ii).int, fedpars.intsmoothalg, ...
fedpars.intsmoothsize, fedpars.intsmoothstd);
minimum_blob_size_add = blob_size_add + 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', ...
blob_size_add(det_num, :), minimum_blob_size_add);
fprintf(' - Projecting volumes: ')
t = tic();
for ii = 1:nbl
num_chars = fprintf('%03d/%03d', ii, nbl);
mu = gvb(ii).ucblmod;
mu1 = mu(1, :);
mu2 = mu(2, :);
mu3 = mu(3, :);
mu12 = mu1 .* mu2;
mu13 = mu1 .* mu3;
mu23 = mu2 .* mu3;
mu123 = mu12 .* mu3;
ints = zeros(1, numel(gvpow) * 8, 'like', bl(ii).int);
ints( 1 : nv) = gvpow.*(-mu1 - mu2 - mu3 + mu12 + mu13 + mu23 - mu123 + 1);
ints( nv+1 : 2*nv) = gvpow.*(mu3 - mu13 - mu23 + mu123);
ints(2*nv+1 : 3*nv) = gvpow.*(mu2 - mu12 - mu23 + mu123);
ints(3*nv+1 : 4*nv) = gvpow.*(mu23 - mu123);
ints(4*nv+1 : 5*nv) = gvpow.*(mu1 - mu12 - mu13 + mu123);
ints(5*nv+1 : 6*nv) = gvpow.*(mu13 - mu123);
ints(6*nv+1 : 7*nv) = gvpow.*(mu12 - mu123);
ints(7*nv+1 : 8*nv) = gvpow.*mu123;
% Accumulate all intensities in the blob voxel-wise
% 'uvw' needs to be nx3; 'ints' is now 1xnx8
try
bl(ii).int = bl(ii).int + accumarray(uvw{ii}, ints, bbsize(ii,:));
catch err
disp(' ')
disp('ERROR')
disp(['The error is probably caused by the projected intensities' ...
' falling outside the blob volume. Try increasing the blob' ...
' volume padding: fedpars.detector(det_ind).blobsizeadd'])
disp('Blob ID:')
disp(ii)
disp('Blob size:')
disp(bbsize(ii,:))
disp('Min projected U,V,W coordinates:')
disp(min(uvw{ii},[],1))
disp('Max projected U,V,W coordinates:')
disp(max(uvw{ii},[],1))
rethrow(err)
end
% % Smooth blob
% if ~strcmp(fedpars.intsmoothalg, 'none')
% bl(ii).int = smooth3(bl(ii).int, fedpars.intsmoothalg, ...
% fedpars.intsmoothsize, fedpars.intsmoothstd);
% end
% % Intensity difference: "measured" - simulated
% bl(ii).intd = bl(ii).int - bl(ii).intm;
%
% % Mean error
% bl(ii).err = sum(abs(bl(ii).intd(:)));
%
% % Correct no. of pixels
% gvb(ii).ucbl8ok = sum(gvb(ii).ucbl8ini == gvb(ii).ucbl8);
%
% % Non-zero pixels
% gvb(ii).ucbl8nonz = sum( bl(ii).intm(gvb(ii).ucbl8) > 0 );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Center of mass and extension of projected blobs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bl(ii).combl = NaN(1,3);
bl(ii).comim = NaN(1,3);
fprintf(repmat('\b', [1, num_chars]));
end
% Intensity difference: "measured" - simulated
bl(ii).intd = bl(ii).int - bl(ii).intm;
% Mean error
bl(ii).err = sum(abs(bl(ii).intd(:)));
% Correct no. of pixels
gvb(ii).ucbl8ok = sum(gvb(ii).ucbl8ini == gvb(ii).ucbl8);
% Non-zero pixels
gvb(ii).ucbl8nonz = sum( bl(ii).intm(gvb(ii).ucbl8) > 0 );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Center of mass and extension of projected blobs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bl(ii).combl = NaN(1,3);
bl(ii).comim = NaN(1,3);
fprintf(repmat('\b', [1, num_chars]));
fprintf('Done in %f s\n', toc(t));
end
fprintf('Done in %f s\n', toc(t));
end
7_fed2/gtFedTestLoadData_v2.m
+ 4
− 3
View file @ e5cb15c0
......@@ -337,7 +337,7 @@ for n = 1:num_dets
[], detgeo(n), labgeo, samgeo, recgeo(n), ...
'ASTRA_grain');
bl(ii).intm = zeros(bl(ii).bbsize, fedpars.bltype);
bl(ii).int = 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
......@@ -391,8 +391,9 @@ for n = 1:num_dets
% 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;
bl(ii).int(:) = [];
% bl(ii).int(:) = 0;
% bl(ii).intd(:) = bl(ii).int - bl(ii).intm;
bl(ii).comim = NaN(1, 3);
bl(ii).combl = NaN(1, 3);
......
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