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

6D-PhantomGeneration: added creation of projection data from cluster of grains 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 97cb0b75
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m 0 → 100644
+ 286
− 0
View file @ c0e55a16
function [proj, refor, or] = gt6DCreateProjDataFromGrainCluster(grs_list, phase_id)
% FUNCTION [proj, refor, or] = gt6DCreateProjDataFromGrainCluster(grs_list, phase_id)
% proj: is a grain.proj structure
% refor: is a grain structure for the average orientation in the average
% point
% or: is a set of grain structures for the extreeme corners of the
% orientation space that should be considered
if (~exist('phase_id', 'var'))
phase_id = 1;
end
num_grains = numel(grs_list);
if (~isstruct(grs_list))
fprintf('Loading grains: ')
for ii_g = num_grains:-1:1
num_chars = fprintf('%02d/%02d (%d)', num_grains-ii_g+1, num_grains, grs_list(ii_g));
grs(ii_g) = gtLoadGrain(phase_id, grs_list(ii_g));
fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars]));
end
fprintf('Done.\n')
else
grs = grs_list;
end
p = gtLoadParameters();
symm = gtCrystGetSymmetryOperators(p.cryst.crystal_system);
cat(1, grs(:).R_vector)
fprintf('Reciprocal disorientations:\n')
for ii_g_1 = 1:num_grains
grs_to_check = (ii_g_1+1):num_grains;
if (numel(grs_to_check) > 0)
fprintf(' - Disorientations from grain: %d\n', grs(ii_g_1).id)
for ii_g_2 = grs_to_check
dis_angle = gtDisorientation(grs(ii_g_1).R_vector', grs(ii_g_2).R_vector', symm);
fprintf(' + Grain %d: %f\n', grs(ii_g_2).id, dis_angle);
end
end
end
refgr = grs(1);
refgr_omind = refgr.allblobs.omind;
refgr_omind_ind = [ ...
find(refgr_omind == 1); find(refgr_omind == 2); ...
find(refgr_omind == 3); find(refgr_omind == 4) ];
for ii_g = num_grains:-1:2
grs(ii_g) = gtGrainAllblobsFilterOrder(grs(ii_g), refgr_omind_ind);
end
if (isfield(refgr.proj, 'ondet'))
ref_ondet = refgr.proj.ondet;
ref_included = ref_ondet(refgr.proj.included);
ref_selected = refgr.proj.selected;
else
ref_ondet = refgr.ondet;
ref_included = ref_ondet(refgr.included);
ref_selected = refgr.selected;
end
space_bboxes = zeros(num_grains, 6);
for ii_g = 1:num_grains
vol_half_size = [grs(ii_g).proj.vol_size_y, grs(ii_g).proj.vol_size_x, grs(ii_g).proj.vol_size_z] / 2;
space_bboxes(ii_g, :) = [ ...
grs(ii_g).proj.centerpix - vol_half_size, ...
grs(ii_g).proj.centerpix + vol_half_size ];
end
estim_space_bbox_pix = [floor(min(space_bboxes(:, 1:3), [], 1)), ceil(max(space_bboxes(:, 4:6), [], 1))]
bbox_size_pix = estim_space_bbox_pix(4:6) - estim_space_bbox_pix(1:3)
estim_space_bbox_mm = [ ...
gtGeoSam2Sam(estim_space_bbox_pix(1:3), p.recgeo, p.samgeo, false, false), ...
gtGeoSam2Sam(estim_space_bbox_pix(4:6), p.recgeo, p.samgeo, false, false) ];
bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3)
for ii_g = num_grains:-1:1
sampler = GtOrientationSampling(p, grs(ii_g));
odfw_R_vectors{ii_g} = sampler.guess_ODF_BB()';
end
r_vecs = cat(1, odfw_R_vectors{:});
estim_orient_bbox = [min(r_vecs, [], 1), max(r_vecs, [], 1)]
bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
% oversizing the orienation a bit
delta_bbox_size_rod = bbox_size_rod * 0.05 ;
estim_orient_bbox = estim_orient_bbox + [-delta_bbox_size_rod, delta_bbox_size_rod]
bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
bbox_size_deg = 2 * atand(bbox_size_rod)
% Let's now compute the bb on the images, by computing for each corner
% of the space bb, the position on the detector of each corner of the
% orientation bb.
or = {};
for ii_x1 = 0:1
base_x1 = estim_space_bbox_mm(1:3) + [ii_x1 * bbox_size_mm(1), 0, 0];
for ii_x2 = 0:1
base_x2 = base_x1 + [0, ii_x2 * bbox_size_mm(2), 0];
for ii_x3 = 0:1
base_x3 = base_x2 + [0, 0, ii_x3 * bbox_size_mm(3)];
for ii_r1 = 0:1
base_r1 = estim_orient_bbox(1:3) + [ii_r1 * bbox_size_rod(1), 0, 0];
for ii_r2 = 0:1
base_r2 = base_r1 + [0, ii_r2 * bbox_size_rod(2), 0];
for ii_r3 = 0:1
base_r3 = base_r2 + [0, 0, ii_r3 * bbox_size_rod(3)];
% fprintf(' center: %f, %f, %f, r_vec: %f, %f, %f\n', base_x3, base_r3)
or{end+1} = struct( ...
'id', refgr.id, 'phaseid', refgr.phaseid, ...
'center', base_x3, 'R_vector', base_r3); %#ok<AGROW>
end
end
end
end
end
end
or = gtCalculateGrain_p(or, p);
or = [or{:}];
num_ors = numel(or);
for ii_g = num_ors:-1:1
or(ii_g) = gtGrainAllblobsFilterOrder(or(ii_g), refgr_omind_ind);
end
or_abs = cat(1, or(:).allblobs);
uvw_tab = zeros(numel(ref_included(ref_selected)), num_ors, 3);
for ii_g = 1:num_ors
if (isfield(or_abs(1), 'detector'))
uvw_tab(:, ii_g, :) = or_abs(ii_g).detector(1).uvw(ref_included(ref_selected), :);
else
uvw_tab(:, ii_g, :) = or_abs(ii_g).uvw(ref_included(ref_selected), :);
end
end
refor = struct( ...
'id', refgr.id, 'phaseid', refgr.phaseid, ...
'center', estim_space_bbox_mm(1:3) + bbox_size_mm / 2, ...
'R_vector', estim_orient_bbox(1:3) + bbox_size_rod / 2 );
refor = gtCalculateGrain(refor, p);
refor = gtGrainAllblobsFilterOrder(refor, refgr_omind_ind);
num_images = gtGetTotNumberOfImages(p);
refor_ws = refor.allblobs.omega(ref_included(ref_selected)) / p.labgeo.omstep;
% Let's treat those blobs at the w edge 360->0
% (from the sampled orientations perspective)
opposite_ws = mod(refor_ws + num_images / 2, num_images);
gt_ref_int = opposite_ws > num_images / 2;
lt_ref_int = ~gt_ref_int;
opposite_ws_repmat = opposite_ws(:, ones(1, size(uvw_tab(:, :, 3), 2)));
uvw_tab(gt_ref_int, :, 3) = uvw_tab(gt_ref_int, :, 3) ...
- num_images .* (uvw_tab(gt_ref_int, :, 3) > opposite_ws_repmat(gt_ref_int, :));
uvw_tab(lt_ref_int, :, 3) = uvw_tab(lt_ref_int, :, 3) ...
+ num_images .* (uvw_tab(lt_ref_int, :, 3) < opposite_ws_repmat(lt_ref_int, :));
img_bboxes = [
squeeze(floor(min(uvw_tab, [], 2))), ...
squeeze( ceil(max(uvw_tab, [], 2))) ];
refor_ns = refor.allblobs.eta(ref_included(ref_selected));
img_sizes = img_bboxes(:, 4:6) - img_bboxes(:, 1:3) + 1;
fprintf('du %8d, dv %8d, dw %8d, eta: %5.3f\n', ...
[img_sizes, refor_ns]');
% We avoid the vertical spots for convenience
inconvenient_etas = acosd(abs(cosd(refor_ns))) < 15;
img_bboxes = img_bboxes(~inconvenient_etas, :);
img_sizes = img_sizes(~inconvenient_etas, :);
fprintf('Loading raw images: ')
num_blobs = size(img_bboxes, 1);
blobs(1:num_blobs) = struct( ...
'intm', [], ... % The normalized blob
'mask', [], ... % The segmentation mask
'bbuim', [], ... % Image BB on U
'bbvim', [], ... % Image BB on V
'bbwim', [], ... % Image BB on W <- Not strictly w
'bbsize', [], ... % Image BoundingBox (includes margins)
'mbbsize', [], ... % Real (Measured) Blob Bounding Box
'mbbu', [], ... % Real (Measured) Blob BB on U
'mbbv', [], ... % Real (Measured) Blob BB on V
'mbbw', [], ... % Real (Measured) Blob BB on W <- Not strictly w
'intensity', [] ); % Blob original intensity
stackUSize = max(img_sizes(:, 1));
stackVSize = max(img_sizes(:, 2));
sel_reflections = find(ref_selected);
sel_reflections = sel_reflections(~inconvenient_etas);
for ii_i = 1:num_blobs
num_chars = fprintf('%02d/%02d', ii_i, num_blobs);
bb = [img_bboxes(ii_i, 1:2), img_sizes(ii_i, 1:2)];
blob_vol = gtGetRawRoi(img_bboxes(ii_i, 3), img_bboxes(ii_i, 6), p.acq, bb);
blob_bb = [img_bboxes(ii_i, 1:3), img_sizes(ii_i, :)];
% Transposing to keep the same convention as spots
blob_vol = permute(blob_vol, [2 1 3]);
blobs(ii_i).mbbsize = blob_bb(4:6);
blobs(ii_i).mbbu = [blob_bb(1), blob_bb(1) + blob_bb(4) + 1];
blobs(ii_i).mbbv = [blob_bb(2), blob_bb(2) + blob_bb(5) + 1];
blobs(ii_i).mbbw = [blob_bb(3), blob_bb(3) + blob_bb(6) + 1];
blob_size_im = [stackUSize, stackVSize, blob_bb(6)+2];
shifts_blob = gtFwdSimGetStackShifts(stackUSize, stackVSize, blob_bb, false);
shifts_blob = [shifts_blob.u, shifts_blob.v, 1];
blobs(ii_i).intm = gtPlaceSubVolume( ...
zeros(blob_size_im, 'single'), blob_vol, shifts_blob);
blobs(ii_i).bbsize = blob_size_im;
blob_bb_im = [blob_bb(1:3) - shifts_blob, blob_size_im];
blobs(ii_i).bbuim = [blob_bb_im(1), blob_bb_im(1) + blob_bb_im(4) - 1];
blobs(ii_i).bbvim = [blob_bb_im(2), blob_bb_im(2) + blob_bb_im(5) - 1];
blobs(ii_i).bbwim = [blob_bb_im(3), blob_bb_im(3) + blob_bb_im(6) - 1];
% We should build the mask from the segmented reflections of the
% grains
% for the moment we only use one (the reference one), to
% renormalize the different reflections
blobs(ii_i).mask = false(blob_size_im);
for ii_g = 1
gbl = refgr.proj.bl(sel_reflections(ii_i));
% gbl = grs(ii_g).proj.bl(sel_reflections(ii_i));
grain_mask = gbl.mask;
gbl_shifts_im = [gbl.bbuim(1), gbl.bbvim(1), gbl.bbwim(1)];
blob_shifts_im = [ ...
blobs(ii_i).bbuim(1), blobs(ii_i).bbvim(1), blobs(ii_i).bbwim(1)];
mask_shifts = gbl_shifts_im - blob_shifts_im;
blobs(ii_i).mask = gtPlaceSubVolume(blobs(ii_i).mask, ...
grain_mask, mask_shifts);
end
blob_int = sum(blobs(ii_i).intm(blobs(ii_i).mask));
blobs(ii_i).intm = blobs(ii_i).intm / blob_int;
% Since we need to build the mask from the reflections of the
% grains, and we might not have them all, we in the end reset the
% mask to 1 but we used it for renormalization purpouses
if (true)
blobs(ii_i).mask = gtPlaceSubVolume( ...
false(blob_size_im), true(size(blob_vol)), shifts_blob);
end
blobs(ii_i).intensity = blob_int;
fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars]));
end
fprintf('Done.\n')
spots = arrayfun(@(x){sum(x.intm, 3)}, blobs);
spots = permute(cat(3, spots{:}), [1 3 2]);
proj.centerpix = estim_space_bbox_pix(1:3) + bbox_size_pix([2 1 3]) / 2;
proj.bl = blobs;
proj.stack = spots;
proj.vol_size_x = bbox_size_pix(2);
proj.vol_size_y = bbox_size_pix(1);
proj.vol_size_z = bbox_size_pix(3);
proj.shift = gtFwdSimComputeVolumeShifts(proj, p);
proj.ondet = ref_ondet;
proj.included = ref_included;
finally_selected = false(size(proj.included));
finally_selected(sel_reflections) = true;
proj.selected = finally_selected;
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