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

6D-reconstruction/Twins: updated twin reconstruction creation to handle raw-images 


And this also obsoletes the older function

Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent efe95b89
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrain.m deleted 100644 → 0
+ 0
− 593
View file @ efe95b89
function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFromTwinnedGrain(grs_list, phase_id, varargin)
% FUNCTION [refor] = gt6DCreateProjDataFromGrainCluster(grs_list, phase_id, varargin)
% 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
conf = struct( ...
'verbose', false, ...
'min_eta', 15, ...
'ospace_oversize', 1.1, ...
'rspace_oversize', 1.1, ...
'flat_normalization', true, ...
'check_twin_spots', false, ...
'mask_spots', true, ...
'oversize', 1.4, ...
'save', false );
conf = parse_pv_pairs(conf, varargin);
p = gtLoadParameters();
symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
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));
grs(ii_g) = gtCalculateGrain(grs(ii_g), p);
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;
grs_list = cat(2, grs(:).id);
end
if (conf.verbose)
fprintf('R_vectors:\n')
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
end
refgr = grs(1);
grs(2:end) = reorder_reflections(grs(2:end), refgr);
twingr = grs(2);
if (isfield(refgr.proj, 'ondet'))
ref_ondet = refgr.proj.ondet;
ref_included = refgr.proj.included;
ref_selected = refgr.proj.selected;
twin_ondet = twingr.proj.ondet;
twin_inc = twingr.proj.included;
twin_sel = twingr.proj.selected;
else
ref_ondet = refgr.ondet;
ref_included = refgr.included;
ref_selected = refgr.selected;
twin_ondet = twingr.ondet;
twin_inc = twingr.included;
twin_sel = twingr.selected;
end
if (conf.verbose)
try
refgr = gtComputeIncomingBeamIntensity(refgr, p);
beam_ints = refgr.beam_intensity(ref_selected);
beam_ints = beam_ints / mean(beam_ints);
catch mexc
gtPrintException(mexc, 'Skipping beam intensity renormalization')
beam_ints = 1;
end
blob_ints = refgr.intensity(ref_selected);
vis_spots = ref_ondet(ref_included(ref_selected));
thetatypes = refgr.allblobs.thetatype(vis_spots);
if (isfield(p.cryst(phase_id), 'int_exp'))
predicted_ints = p.cryst(phase_id).int_exp(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
end
predicted_ints = p.cryst(phase_id).int(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
fprintf('%d) %d, %g\n', ...
[(1:numel(blob_ints))', refgr.allblobs.thetatype(vis_spots), ...
blob_ints - predicted_ints / mean(predicted_ints) * mean(blob_ints)]')
hold(ax, 'off')
end
% let's now filter the reflections that are not available for the twin
bool_ref_ondet = false(size(refgr.allblobs.omega));
bool_ref_ondet(ref_ondet) = true;
bool_ref_inc = false(size(refgr.allblobs.omega));
bool_ref_inc(ref_ondet(ref_included)) = true;
bool_ref_sel = false(size(refgr.allblobs.omega));
bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;
bool_twin_ondet = false(size(refgr.allblobs.omega));
bool_twin_ondet(twin_ondet) = true;
bool_twin_inc = false(size(refgr.allblobs.omega));
bool_twin_inc(twin_ondet(twin_inc)) = true;
bool_twin_sel = false(size(refgr.allblobs.omega));
bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;
% bool_deleted_ref_ondet = bool_ref_ondet & ~bool_twin_ondet;
% bool_deleted_ref_inc = bool_ref_inc & ~bool_twin_ondet;
%
% ref_deleted_ondet = find(bool_deleted_ref_ondet);
% ref_deleted_included = find(bool_deleted_ref_inc(ref_deleted_ondet));
bool_ref_ondet = bool_ref_ondet & bool_twin_ondet;
bool_ref_inc = bool_ref_inc & bool_twin_ondet;
bool_ref_sel = bool_ref_sel & bool_twin_ondet;
% bool_twin_ondet = bool_ref_ondet & bool_twin_ondet;
% bool_twin_inc = bool_ref_inc & bool_twin_inc;
% bool_twin_sel = bool_ref_sel & bool_twin_sel;
ref_ondet = find(bool_ref_ondet);
ref_included = find(bool_ref_inc(ref_ondet));
ref_selected = bool_ref_sel(ref_ondet(ref_included));
if (isfield(refgr.proj, 'ondet'))
refgr.proj.ondet = ref_ondet;
refgr.proj.included = ref_included;
refgr.proj.selected = ref_selected;
else
refgr.ondet = ref_ondet;
refgr.included = ref_included;
refgr.selected = ref_selected;
end
% ref_twin_ondet_bool = bool_twin_ondet(ref_ondet);
% ref_twin_included_bool = bool_twin_inc(ref_ondet(ref_included));
% ref_twin_selected_bool = bool_twin_sel(ref_ondet(ref_included(ref_selected)));
% And now..
fwdsim_inc_refl = ref_ondet(ref_included);
if (conf.verbose)
vis_sel = fwdsim_inc_refl;
% vis_sel = ref_ondet;
fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_inds (%d<->%d), hkls [%2d %2d %2d]<->[%2d %2d %2d] <- diff: %7.2f (!! %d, L %g) [twin: i %d, s %d]\n', ...
[(1:numel(vis_sel))', ...
grs(1).allblobs.omega(vis_sel), ...
grs(2).allblobs.omega(vis_sel), ...
grs(1).allblobs.eta(vis_sel), ...
grs(2).allblobs.eta(vis_sel), ...
grs(1).allblobs.omind(vis_sel), ...
grs(2).allblobs.omind(vis_sel), ...
grs(1).allblobs.hklsp(vis_sel, [1 2 4]), ...
grs(2).allblobs.hklsp(vis_sel, [1 2 4]), ...
abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180), ...
abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180) < grs(1).allblobs.lorentzfac(vis_sel), ...
grs(1).allblobs.lorentzfac(vis_sel), ...
bool_twin_inc(vis_sel), bool_twin_sel(vis_sel)]');
end
vol_half_size = [refgr.proj.vol_size_y, refgr.proj.vol_size_x, refgr.proj.vol_size_z] / 2;
space_bbox = [ ...
refgr.proj.centerpix - vol_half_size * conf.rspace_oversize, ...
refgr.proj.centerpix + vol_half_size * conf.rspace_oversize ];
estim_space_bbox_pix = [floor(space_bbox(1:3)), ceil(space_bbox(4:6))];
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);
img_bboxes = cell(num_grains, 1);
img_sizes = cell(num_grains, 1);
for ii_g = num_grains:-1:1
sampler = GtOrientationSampling(p, grs(ii_g), 'verbose', conf.verbose);
r_vecs = sampler.guess_ODF_BB()';
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);
if (conf.verbose)
fprintf('\n');
fprintf('Estimated spatial voxel BBox: [%3d, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix);
fprintf(' BBox size: %3d, %3d, %3d (%f, %f, %f mm)\n', bbox_size_pix, bbox_size_mm);
fprintf(' Estimated orientation BBox: [%3.3f, %3.3f, %3.3f] -> [%3.3f, %3.3f, %3.3f]\n', estim_orient_bbox);
fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
fprintf('\n');
end
% 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 = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p);
refor = struct( ...
'id', refgr.id, 'phaseid', refgr.phaseid, 'gids', grs_list, ...
'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 = reorder_reflections(refor, refgr);
uvw_tab = get_uvw_tab(or, fwdsim_inc_refl);
refor_ws = refor.allblobs.omega(fwdsim_inc_refl) / p.labgeo.omstep;
uvw_tab = recenter_w_tab(uvw_tab, refor_ws, p);
img_bboxes{ii_g} = [
squeeze(floor(min(uvw_tab, [], 2))), ...
squeeze( ceil(max(uvw_tab, [], 2))) ];
img_sizes{ii_g} = img_bboxes{ii_g}(:, 4:6) - img_bboxes{ii_g}(:, 1:3) + 1;
refor.bb_ors = or;
samp_ors(ii_g) = refor;
end
% We avoid the vertical spots for convenience
bad_etas = false(size(fwdsim_inc_refl));
for ii_g = 1:num_grains
refor_ns = samp_ors(ii_g).allblobs.eta(fwdsim_inc_refl);
bad_etas = bad_etas | acosd(abs(cosd(refor_ns))) < conf.min_eta;
if (conf.verbose)
fprintf('Grain: %d\n', ii_g)
fprintf('%02d) du %8d, dv %8d, dw %8d, eta: %7.3f <- used: %d, u: [%4d %4d], v: [%4d %4d], w: [%4d %4d]\n', ...
[(1:numel(bad_etas))', img_sizes{1}, refor_ns, ~bad_etas, img_bboxes{1}(:, [1 4 2 5 3 6]) ]');
fprintf('\n')
end
end
for ii_g = 1:num_grains
img_bboxes{ii_g} = img_bboxes{ii_g}(~bad_etas, :);
img_sizes{ii_g} = img_sizes{ii_g}(~bad_etas, :);
end
[stackUSize, stackVSize] = get_stack_UV_size(cat(1, img_sizes{:}), p, conf);
inc_reflections = find(1:numel(fwdsim_inc_refl));
inc_reflections = inc_reflections(~bad_etas);
ref_BBus = cat(1, refgr.proj.bl(inc_reflections).mbbu);
ref_BBvs = cat(1, refgr.proj.bl(inc_reflections).mbbv);
ref_BBsizes = cat(1, refgr.proj.bl(inc_reflections).mbbsize);
ref_norm_BBs = [ref_BBus(:, 1), ref_BBvs(:, 1), ref_BBsizes(:, 1:2)];
ref_BBus = (ref_BBus + cat(1, refgr.proj.bl(inc_reflections).bbuim)) / 2;
ref_BBvs = (ref_BBvs + cat(1, refgr.proj.bl(inc_reflections).bbvim)) / 2;
ref_BBsizes = (ref_BBsizes + cat(1, refgr.proj.bl(inc_reflections).bbsize)) / 2;
ref_cut_BBs = [ref_BBus(:, 1), ref_BBvs(:, 1), ref_BBsizes(:, 1:2)];
inc_ref_ondet = ref_ondet(ref_included(inc_reflections));
inc_ref_omegas = samp_ors(1).allblobs.omega(inc_ref_ondet);
inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(inc_ref_ondet);
for ii_g = 1:num_grains
fprintf('Loading raw images for grain %d: ', ii_g)
if (ii_g > 1)
w_diff = inc_ref_omegas - samp_ors(ii_g).allblobs.omega(inc_ref_ondet);
% Could be improved
shared_parent = abs(mod(w_diff + 180, 360) - 180) < inc_ref_Lfac;
end
temp_gr = grs(ii_g);
temp_gr.center = refgr.center;
proj_bls = gtFwdSimPredictProjectedGrainBBox(temp_gr, ref_norm_BBs, inc_ref_omegas, inc_ref_ondet, p);
proj_cut_bls = gtFwdSimPredictProjectedGrainBBox(temp_gr, ref_cut_BBs, inc_ref_omegas, inc_ref_ondet, p);
num_blobs = numel(inc_reflections);
blobs(1:num_blobs) = gtFwdSimBlobDefinition();
for ii_b = 1:num_blobs
num_chars = fprintf('%02d/%02d', ii_b, num_blobs);
if (ii_g == 1 || ~shared_parent(ii_b))
gbl_norm = proj_bls(ii_b);
gbl_cut = proj_cut_bls(ii_b);
blobs(ii_b) = load_blob(img_bboxes{ii_g}(ii_b, :), ...
img_sizes{ii_g}(ii_b, :), stackUSize, stackVSize, gbl_norm, gbl_cut, p); %#ok<AGROW>
else
% Using the same because it gives the same size/uvw-limits
blobs(ii_b) = samp_ors(1).proj.bl(ii_b); %#ok<AGROW>
end
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;
proj.bl = blobs;
proj.stack = spots;
vol_size = bbox_size_pix + mean(bbox_size_pix) * 0.3;
proj.vol_size_x = vol_size(2);
proj.vol_size_y = vol_size(1);
proj.vol_size_z = vol_size(3);
proj.ondet = ref_ondet;
proj.included = ref_included(inc_reflections);
proj.selected = true(size(proj.included));
if (ii_g > 1)
proj.selected = samp_ors(1).proj.selected;
proj.shared_parent = shared_parent; % <- size: proj.included
else
proj.selected = true(size(proj.included));
end
samp_ors(ii_g).proj = proj;
if (conf.check_twin_spots)
samp_ors(ii_g).proj.selected = samp_ors(1).proj.selected & gtGuiGrainMontage(samp_ors(ii_g), [], true);
for ii_g_o = 1:ii_g-1
samp_ors(ii_g_o).proj.selected = samp_ors(ii_g_o).proj.selected & samp_ors(ii_g).proj.selected;
end
end
end
% partial_ints = zeros(numel(inc_reflections), num_grains);
% phys_renorm = ones(numel(inc_reflections), num_grains);
%
% partial_ints(:, 1) = cat(1, samp_ors(1).proj.bl.intensity);
% for ii_g = 2:num_grains
% to_be_summed = ~samp_ors(ii_g).proj.shared_parent;
% partial_ints(to_be_summed, ii_g) = cat(1, samp_ors(ii_g).proj.bl(to_be_summed).intensity);
% end
%
% if (~conf.flat_normalization)
% for ii_g = 1:num_grains
% try
% temp_or = gtComputeIncomingBeamIntensity(samp_ors(ii_g), p);
% beam_ints = temp_or.beam_intensity;
% catch mexc
% gtPrintException(mexc, 'Turning beam inensity correction off, because interlaced acquisition is not supported')
% beam_ints = 1;
% end
% phys_renorm(:, ii_g) = 1 ...
% ./ samp_ors(ii_g).allblobs.lorentzfac(inc_ref_ondet) ...
% ./ (beam_ints / mean(beam_ints));
% end
% end
% total_ints = sum(partial_ints .* phys_renorm, 2);
%
% % Let's renormalize
% for ii_g = 1:num_grains
% for ii_b = 1:num_blobs
% ints_renorm = phys_renorm(:, ii_g) ./ total_ints(:);
%
% if (ii_g > 1 && samp_ors(ii_g).proj.shared_parent(ii_b))
% fprintf('Shared!\n')
% else
% samp_ors(ii_g).proj.bl(ii_b).intm = samp_ors(ii_g).proj.bl(ii_b).intm * ints_renorm(ii_b);
% end
% end
% end
if (conf.verbose)
f = figure();
ax = axes('parent', f);
hold(ax, 'on');
for ii_g = 1:num_grains
gt6DPlotOrientationBBox(ax, cat(1, samp_ors(ii_g).bb_ors(:).R_vector));
scatter3(ax, grs(ii_g).R_vector(1), grs(ii_g).R_vector(2), grs(ii_g).R_vector(3), 30);
end
hold(ax, 'off');
drawnow();
end
if (conf.save)
str_ids = sprintf('_%04d', grs_list);
grain_filename = fullfile(p.acq.dir, '4_grains', ...
sprintf('phase_%02d', phase_id), ...
sprintf('grain_cluster%s.mat', str_ids));
fprintf('Saving the cluster file "%s"..', grain_filename)
save(grain_filename, 'samp_ors', '-v7.3');
fprintf('\b\b: Done.\n')
fprintf('Saving to sample.mat..')
sample = GtSample.loadFromFile();
sample.phases{phase_id}.clusters(end+1) = ...
GtPhase.makeCluster(grs_list, refor.R_vector, estim_space_bbox_pix, estim_orient_bbox, 1);
sample.saveToFile();
fprintf('\b\b: Done.\n')
end
end
function orients = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p)
orients = {};
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)
orients{end+1} = struct( ...
'id', refgr.id, 'phaseid', refgr.phaseid, ...
'center', base_x3, 'R_vector', base_r3); %#ok<AGROW>
end
end
end
end
end
end
orients = gtCalculateGrain_p(orients, p);
orients = [orients{:}];
orients = reorder_reflections(orients, refgr);
end
function orients = reorder_reflections(orients, refgr)
refgr_omind = refgr.allblobs.omind;
for ii_g = numel(orients):-1:1
orients(ii_g) = gtGrainAllblobsFilterOrder(orients(ii_g), refgr_omind);
end
end
function uvw_tab = recenter_w_tab(uvw_tab, refor_ws, p)
num_images = gtGetTotNumberOfImages(p);
% 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, :));
end
function uvw_tab = get_uvw_tab(orientations, refl_sel)
num_ors = numel(orientations);
or_abs = cat(1, orientations(:).allblobs);
uvw_tab = zeros(numel(refl_sel), num_ors, 3);
for ii_o = 1:num_ors
if (isfield(or_abs(1), 'detector'))
uvw_tab(:, ii_o, :) = or_abs(ii_o).detector(1).uvw(refl_sel, :);
else
uvw_tab(:, ii_o, :) = or_abs(ii_o).uvw(refl_sel, :);
end
end
end
function [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf)
max_img_sizes = [max(img_sizes(:, 1)), max(img_sizes(:, 2))];
stack_imgs_oversize = min(p.fsim.oversize, conf.oversize);
u_size = round(max_img_sizes(1) * stack_imgs_oversize);
v_size = round(max_img_sizes(2) * stack_imgs_oversize);
if (conf.verbose)
fprintf('\n');
fprintf(' Maximum images size: [%3d, %3d]\n', max_img_sizes);
fprintf('Stack images size (oversize: %1.2f): [%3d, %3d]\n', stack_imgs_oversize, u_size, v_size);
fprintf('\n');
end
end
function blob = load_blob(img_bboxes, img_sizes, stackUSize, stackVSize, gbl_norm, gbl_cut, p)
blob = gtFwdSimBlobDefinition();
bb = [img_bboxes(1, 1:2), img_sizes(1, 1:2)];
blob_vol = gtGetRawRoi(img_bboxes(3), img_bboxes(6), p.acq, bb);
blob_vol(blob_vol < 0) = 0;
blob_vol(isnan(blob_vol)) = 0;
blob_bb = [img_bboxes(1, 1:3), img_sizes(1, :)];
% Transposing to keep the same convention as spots
blob_vol = permute(blob_vol, [2 1 3]);
blob.mbbsize = blob_bb(4:6);
blob.mbbu = [blob_bb(1), blob_bb(1) + blob_bb(4) + 1];
blob.mbbv = [blob_bb(2), blob_bb(2) + blob_bb(5) + 1];
blob.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];
blob.intm = gtPlaceSubVolume( ...
zeros(blob_size_im, 'single'), single(blob_vol), shifts_blob);
blob.bbsize = blob_size_im;
blob_bb_im = [blob_bb(1, 1:3) - shifts_blob, blob_size_im];
blob.bbuim = [blob_bb_im(1), blob_bb_im(1) + blob_bb_im(4) - 1];
blob.bbvim = [blob_bb_im(2), blob_bb_im(2) + blob_bb_im(5) - 1];
blob.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
blob.mask = false(blob_size_im);
% spreading the flat mask over the full expected BBox in UVW
grain_mask = gbl_cut.mask(:, :, ones(size(blob.mask, 3)));
gbl_shifts_im = [gbl_cut.bbuim(1), gbl_cut.bbvim(1), blob.bbwim(1)]; % <- note w component: blob.bbwim(1)
blob_shifts_im = [blob.bbuim(1), blob.bbvim(1), blob.bbwim(1)];
mask_shifts = gbl_shifts_im - blob_shifts_im;
blob.mask = gtPlaceSubVolume(blob.mask, grain_mask, mask_shifts);
blob.intm(~blob.mask) = 0;
% 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
blob.mask = false(blob_size_im);
% spreading the flat mask over the full expected BBox in UVW
grain_mask = gbl_norm.mask(:, :, ones(size(blob.mask, 3)));
gbl_shifts_im = [gbl_norm.bbuim(1), gbl_norm.bbvim(1), blob.bbwim(1)]; % <- note w component: blob.bbwim(1)
blob_shifts_im = [blob.bbuim(1), blob.bbvim(1), blob.bbwim(1)];
mask_shifts = gbl_shifts_im - blob_shifts_im;
blob.mask = gtPlaceSubVolume(blob.mask, grain_mask, mask_shifts);
blob_int = sum(blob.intm(blob.mask));
% blob.intm = blob.intm / blob_int;
blob.intensity = blob_int;
end
zUtil_Deformation/gt6DCreateProjDataFromTwinnedGrainFwdProj.m
+ 344
− 170
View file @ 84d51cca
...@@ -21,6 +21,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat ...@@ -21,6 +21,7 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
'ospace_oversize', 1.1, ... 'ospace_oversize', 1.1, ...
'check_spots', false, ... 'check_spots', false, ...
'stack_oversize', 1.4, ... 'stack_oversize', 1.4, ...
'use_raw_images', false, ...
'save', false ); 'save', false );
conf = parse_pv_pairs(conf, varargin); conf = parse_pv_pairs(conf, varargin);
...@@ -128,72 +129,18 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat ...@@ -128,72 +129,18 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
ref_included = refgr_proj.included; ref_included = refgr_proj.included;
ref_selected = refgr_proj.selected; ref_selected = refgr_proj.selected;
if (conf.verbose) if (conf.verbose > 1 || true)
vis_sel = ref_ondet(ref_included); produce_matching_reflections_table(grs, conf, refl_matches, ref_ondet, ref_included)
w_diffs = abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180);
fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_ind %d, hkl [%2d %2d %2d] <- matches %d, diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
[(1:numel(vis_sel))', ...
grs(1).allblobs.omega(vis_sel), ...
grs(2).allblobs.omega(vis_sel), ...
grs(1).allblobs.eta(vis_sel), ...
grs(2).allblobs.eta(vis_sel), ...
grs(1).allblobs.omind(vis_sel), ...
grs(1).allblobs.hklsp(vis_sel, [1 2 4]), ...
refl_matches(vis_sel, 1), ...
w_diffs, ...
w_diffs < grs(1).allblobs.lorentzfac(vis_sel), ...
grs(1).allblobs.lorentzfac(vis_sel), ...
acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
]');
end end
if (conf.verbose && false) if (conf.verbose > 1 && false)
try produce_beam_intensity_renorm_figure(refgr, p, ref_ondet, ref_included, ref_selected)
refgr = gtComputeIncomingBeamIntensity(refgr, p);
beam_ints = refgr.beam_intensity(ref_selected);
beam_ints = beam_ints / mean(beam_ints);
catch mexc
gtPrintException(mexc, 'Skipping beam intensity renormalization')
beam_ints = 1;
end
blob_ints = refgr.intensity(ref_selected);
vis_spots = ref_ondet(ref_included(ref_selected));
thetatypes = refgr.allblobs.thetatype(vis_spots);
if (isfield(p.cryst(phase_id), 'int_exp'))
predicted_ints = p.cryst(phase_id).int_exp(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
end
predicted_ints = p.cryst(phase_id).int(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
fprintf('%d) %d, %g\n', ...
[(1:numel(blob_ints))', refgr.allblobs.thetatype(vis_spots), ...
blob_ints - predicted_ints / mean(predicted_ints) * mean(blob_ints)]')
hold(ax, 'off')
end end
fprintf('Determining real-space bounding box..')
c = tic();
% Determining real-space bounding box. It is assumed to be contiguos.
[gr_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(grs, conf.det_index, conf.rspace_oversize); [gr_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(grs, conf.det_index, conf.rspace_oversize);
gr_center_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo, p.samgeo, false, false); gr_center_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo, p.samgeo, false, false);
...@@ -202,15 +149,22 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat ...@@ -202,15 +149,22 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
gtGeoSam2Sam(estim_space_bbox_pix(4:6), 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); bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3);
if (conf.verbose) fprintf('\b\b: Done in %g seconds.\n', toc(c))
fprintf('\n');
fprintf('Estimated spatial voxel BBox: [%3d, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix); fprintf(' - Estimated extent: [%3g, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix);
fprintf(' BBox size: %3d, %3d, %3d (%f, %f, %f mm)\n', bbox_size_pix, bbox_size_mm); fprintf(' BBox size: %3d, %3d, %3d (%g, %g, %g mm)\n', bbox_size_pix, bbox_size_mm);
end
max_img_sizes = zeros(num_grains, 2);
img_sizes = zeros(num_grains, 2); extended_projs = get6DExtendedProjDefinition(num_grains);
samp_ors = cell(num_grains, 1);
% Determining orientation-space bounding boxes. They are assumed to be
% detached from each others.
for ii_g = 1:num_grains
fprintf('Determining orientation-space region: %d/%d\n', ii_g, num_grains)
for ii_g = num_grains:-1:1
sampler = GtOrientationSampling(p, grs(ii_g), ... sampler = GtOrientationSampling(p, grs(ii_g), ...
'detector_index', conf.det_index, 'verbose', conf.verbose); 'detector_index', conf.det_index, 'verbose', conf.verbose);
r_vecs = sampler.guess_ODF_BB()'; r_vecs = sampler.guess_ODF_BB()';
...@@ -225,11 +179,8 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat ...@@ -225,11 +179,8 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
bbox_size_deg = 2 * atand(bbox_size_rod); bbox_size_deg = 2 * atand(bbox_size_rod);
if (conf.verbose) fprintf(' - Estimated extent: [%g, %g, %g] -> [%g, %g, %g]\n', estim_orient_bbox);
fprintf('%d) Estimated orientation BBox: [%3.3f, %3.3f, %3.3f] -> [%3.3f, %3.3f, %3.3f]\n', ii_g, estim_orient_bbox); fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
fprintf('\n');
end
refor = struct( ... refor = struct( ...
'id', grs(ii_g).id, 'phaseid', grs(ii_g).phaseid, ... 'id', grs(ii_g).id, 'phaseid', grs(ii_g).phaseid, ...
...@@ -239,126 +190,122 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat ...@@ -239,126 +190,122 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
refor.bb_ors = gt6DCalculateFullSpaceCorners(estim_space_bbox_mm, ... refor.bb_ors = gt6DCalculateFullSpaceCorners(estim_space_bbox_mm, ...
bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p); bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p);
img_sizes(ii_g, :) = [ ... samp_ors{ii_g} = refor;
size(grs(ii_g).proj(conf.det_index).stack, 1), ...
size(grs(ii_g).proj(conf.det_index).stack, 3), ];
samp_ors(ii_g) = refor; % Populating the extended projection data-structures
if (ii_g > 1)
extended_projs(ii_g) = find_matches_with_refor(refgr, ...
grs(ii_g), refl_matches(:, ii_g-1), conf);
else
extended_projs(ii_g).ondet = ref_ondet;
extended_projs(ii_g).included = ref_included;
extended_projs(ii_g).selected = ref_selected;
end
end end
[stackUSize, stackVSize] = get_stack_UV_size(cat(1, img_sizes), p, conf); samp_ors = [samp_ors{:}];
bool_ref_inc = false(size(refgr.allblobs.omega)); fprintf('Determining UVW bounding boxes..')
bool_ref_inc(ref_ondet(ref_included)) = true; c = tic();
bool_ref_sel = false(size(refgr.allblobs.omega));
bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;
% At the moment we only support the configuration: parent + twins if (conf.use_raw_images)
refl_matches = gt6DGetMatchingReflections(refgr, grs(2:end), 2); uvw_tab = cell(num_grains, 1);
img_bboxes = cell(num_grains, 1);
img_sizes = cell(num_grains, 1);
for ii_g = 1:num_grains for ii_g = 1:num_grains
fprintf('Loading raw images for grain %d: ', ii_g) local_sel_refl = extended_projs(ii_g).ondet(extended_projs(ii_g).included(extended_projs(ii_g).selected));
if (ii_g > 1)
twingr = grs(ii_g);
twingr_proj = twingr.proj(conf.det_index);
twin_ondet = twingr_proj.ondet;
twin_inc = twingr_proj.included;
twin_sel = twingr_proj.selected;
bool_twin_ondet = false(size(refgr.allblobs.omega));
bool_twin_ondet(twin_ondet) = true;
bool_twin_inc = false(size(refgr.allblobs.omega));
bool_twin_inc(twin_ondet(twin_inc)) = true;
bool_twin_sel = false(size(refgr.allblobs.omega));
bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;
% Indices of reflections that exist on the detector for the uvw_tab{ii_g} = zeros(numel(local_sel_refl), numel(samp_ors(ii_g).bb_ors), 3);
% twin, and were included for the parent for ii_o = 1:numel(samp_ors(ii_g).bb_ors)
bool_test_inc = bool_ref_inc & bool_twin_ondet; uvw_tab{ii_g}(:, ii_o, :) = samp_ors(ii_g).bb_ors(ii_o).detector(conf.det_index).uvw(local_sel_refl, :);
end
if (false) % Let's treat those blobs at the w edge 360->0
inc_ref_omegas = samp_ors(1).allblobs.omega(bool_test_inc); % (from the sampled orientations perspective)
inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(bool_test_inc); refor_ws = refor.allblobs.omega(local_sel_refl) / gtGetOmegaStepDeg(p, conf.det_index);
inc_ref_etas = samp_ors(1).allblobs.eta(bool_test_inc); uvw_tab{ii_g}(:, :, 3) = gtGrainAnglesTabularFix360deg(uvw_tab{ii_g}(:, :, 3), refor_ws, p);
inc_twin_omegas = samp_ors(ii_g).allblobs.omega(bool_test_inc); img_bboxes{ii_g} = [ ...
inc_twin_etas = samp_ors(ii_g).allblobs.eta(bool_test_inc); min(uvw_tab{ii_g}(:, :, 1:2), [], 2), ...
max(uvw_tab{ii_g}(:, :, 1:2), [], 2) ];
w_diff = abs(mod(inc_ref_omegas - inc_twin_omegas + 180, 360) - 180); img_sizes{ii_g} = img_bboxes{ii_g}(:, 4:5) - img_bboxes{ii_g}(:, 1:2) + 1;
n_diff = abs(mod(inc_ref_etas - inc_twin_etas + 180, 360) - 180); img_sizes{ii_g} = reshape(img_sizes{ii_g}, [], 2);
test_shared = ... max_img_sizes(ii_g, :) = max(img_sizes{ii_g}, [], 1);
(w_diff < conf.angular_toll * inc_ref_Lfac) ... end
& (n_diff < conf.angular_toll); else
else for ii_g = 1:num_grains
test_shared = refl_matches(bool_test_inc, ii_g-1); max_img_sizes(ii_g, :) = [ ...
end size(grs(ii_g).proj(conf.det_index).stack, 1), ...
size(grs(ii_g).proj(conf.det_index).stack, 3), ];
end
end
% let's add the shared as included [stackUSize, stackVSize] = get_stack_UV_size(max_img_sizes, p, conf);
shared_parent = false(size(refgr.allblobs.omega));
bool_test_inc_indx = find(bool_test_inc);
shared_parent(bool_test_inc_indx(test_shared)) = true;
bool_twin_inc_redundant = bool_twin_inc & shared_parent; fprintf('\b\b: Done in %g seconds.\n', toc(c))
% So if it is shared, we use the parent's blob
bool_twin_inc = bool_twin_inc | shared_parent;
% If the parent enables it, we enable it as well
bool_twin_sel = bool_twin_sel | (shared_parent & bool_ref_sel);
% We now find the blobs in the old twin bl structure that are for ii_g = 1:num_grains
% shared with the parent fprintf('Processing orientation-space region: %d/%d\n', ii_g, num_grains)
shared_parent_inc_redundant = bool_twin_inc_redundant(twin_ondet(twin_inc));
twin_inc = find(bool_twin_inc(twin_ondet)); num_blobs = numel(extended_projs(ii_g).included);
twin_sel = bool_twin_sel(twin_ondet(twin_inc));
shared_parent_pos = bool_ref_inc & shared_parent; c = tic();
shared_parent_pos = shared_parent_pos(ref_ondet(ref_included)); if (conf.use_raw_images)
fprintf(' - Loading raw images: ')
% We cut it down to the size of the included blobs = gtFwdSimBlobDefinition('blob', num_blobs);
shared_parent = shared_parent(twin_ondet(twin_inc)); for ii_b = 1:num_blobs
num_chars = fprintf('%03d/%03d', ii_b, num_blobs);
if (ii_g > 1 && extended_projs(ii_g).shared_refl(ii_b))
blobs(ii_b) = samp_ors(1).proj.bl( ...
extended_projs(ii_g).shared_refl_pos_in_parent(ii_b));
else
blobs(ii_b) = load_blob( ...
img_bboxes{ii_g}(ii_b, :), img_sizes{ii_g}(ii_b, :), ...
stackUSize, stackVSize, p, conf);
end
local_ondet = twin_ondet; fprintf(repmat('\b', [1, num_chars]));
local_included = twin_inc; end
local_selected = twin_sel;
else else
local_ondet = ref_ondet; fprintf(' - Loading segmented blobs..')
local_included = ref_included;
local_selected = ref_selected;
end
num_blobs = numel(local_included);
if (ii_g == 1) if (ii_g == 1)
blobs = grs(ii_g).proj.bl; blobs = grs(ii_g).proj.bl;
else else
blobs = gtFwdSimBlobDefinition('blob', num_blobs); blobs = gtFwdSimBlobDefinition('blob', num_blobs);
blobs(~shared_parent) = grs(ii_g).proj.bl(~shared_parent_inc_redundant); blobs(~extended_projs(ii_g).shared_refl) = grs(ii_g).proj.bl(~extended_projs(ii_g).shared_refl_included_pos);
% Using the same because it gives the same size/uvw-limits % Using the same because it gives the same size/uvw-limits
blobs(shared_parent) = samp_ors(1).proj.bl(shared_parent_pos); blobs(extended_projs(ii_g).shared_refl) = samp_ors(1).proj.bl(extended_projs(ii_g).shared_refl_pos_in_parent);
if (conf.verbose) if (conf.verbose)
parent_ids = grs(1).fwd_sim(conf.det_index).spotid(ref_included(shared_parent_pos))'; parent_ids = grs(1).fwd_sim(conf.det_index).spotid(ref_included(extended_projs(ii_g).shared_refl_pos_in_parent))';
twin_ids = grs(ii_g).fwd_sim(conf.det_index).spotid(twin_inc(shared_parent))'; twin_ids = grs(ii_g).fwd_sim(conf.det_index).spotid(twin_inc(extended_projs(ii_g).shared_refl))';
fprintf('Shared included blobs:\n') fprintf('Shared included blobs:\n')
fprintf(' - %d: parent blob_id %d, twin blob_id %d\n', ... fprintf(' - %d: parent blob_id %d, twin blob_id %d\n', ...
[1:numel(find(shared_parent)); parent_ids; twin_ids]) [1:numel(find(extended_projs(ii_g).shared_refl)); parent_ids; twin_ids])
end
end end
end end
fprintf('\b\b: Done in %g seconds.\n - Equalizing blob sizes..', toc(c))
for ii_b = 1:num_blobs for ii_b = 1:num_blobs
num_chars = fprintf('%02d/%02d', ii_b, num_blobs); num_chars = fprintf('%02d/%02d', ii_b, num_blobs);
blobs(ii_b) = include_blob(blobs(ii_b), stackUSize, stackVSize); blobs(ii_b) = equalize_blob_size(blobs(ii_b), stackUSize, stackVSize);
fprintf(repmat('\b', [1 num_chars])); fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars])); fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars])); fprintf(repmat('\b', [1 num_chars]));
end end
fprintf('Done.\n')
fprintf('\b\b: Done in %g seconds.\n - Producing proj data-structure..', toc(c))
c = tic();
spots = arrayfun(@(x){sum(x.intm, 3)}, blobs); spots = arrayfun(@(x){sum(x.intm, 3)}, blobs);
spots = permute(cat(3, spots{:}), [1 3 2]); spots = permute(cat(3, spots{:}), [1 3 2]);
...@@ -373,21 +320,26 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat ...@@ -373,21 +320,26 @@ function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDat
proj.vol_size_y = vol_size(1); proj.vol_size_y = vol_size(1);
proj.vol_size_z = vol_size(3); proj.vol_size_z = vol_size(3);
proj.ondet = local_ondet; proj.ondet = extended_projs(ii_g).ondet;
proj.included = local_included; proj.included = extended_projs(ii_g).included;
proj.selected = local_selected; proj.selected = extended_projs(ii_g).selected;
% Temporary solution that allows the parent + twin configuration to
% work, but it should be replaced by a more advanced global
% solution for handling shared blobs by different orientation boxes
if (ii_g > 1) if (ii_g > 1)
proj.shared_bls_parent = zeros(size(shared_parent)); proj.shared_bls_parent = zeros(size(extended_projs(ii_g).shared_refl));
proj.shared_bls_parent(shared_parent) = find(shared_parent_pos); proj.shared_bls_parent(extended_projs(ii_g).shared_refl) = find(extended_projs(ii_g).shared_refl_pos_in_parent);
samp_ors(1).proj(conf.det_index).shared_bls_twins(shared_parent_pos) = true; samp_ors(1).proj(conf.det_index).shared_bls_twins(extended_projs(ii_g).shared_refl_pos_in_parent) = true;
else else
proj.shared_bls_twins = false(num_blobs, 1); proj.shared_bls_twins = false(num_blobs, 1);
end end
samp_ors(ii_g).proj(conf.det_index) = proj; samp_ors(ii_g).proj(conf.det_index) = proj;
fprintf('\b\b: Done in %g seconds.\n', toc(c))
if (conf.check_spots && ii_g > 1) if (conf.check_spots && ii_g > 1)
samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true, conf.det_index); samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true, conf.det_index);
end end
...@@ -439,7 +391,7 @@ function [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf) ...@@ -439,7 +391,7 @@ function [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf)
end end
end end
function blob = include_blob(blob, stackUSize, stackVSize) function blob = equalize_blob_size(blob, stackUSize, stackVSize)
new_bbsize = [stackUSize, stackVSize, blob.bbsize(3)]; new_bbsize = [stackUSize, stackVSize, blob.bbsize(3)];
shift = floor((new_bbsize - blob.bbsize) / 2); shift = floor((new_bbsize - blob.bbsize) / 2);
...@@ -479,5 +431,227 @@ function gr = twin_fwd_sim_to_gr(ref_gr, twin) ...@@ -479,5 +431,227 @@ function gr = twin_fwd_sim_to_gr(ref_gr, twin)
gr.intensity = twin.fwd_sim(1).intensity; gr.intensity = twin.fwd_sim(1).intensity;
end end
function produce_matching_reflections_table(grs, conf, refl_matches, ref_ondet, ref_included)
vis_sel = ref_ondet(ref_included);
w_diffs = abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180);
if (size(grs(1).allblobs.hklsp, 2) == 4)
hklsp = grs(1).allblobs.hklsp(vis_sel, [1 2 4]);
else
hklsp = grs(1).allblobs.hklsp(vis_sel, :);
end
fprintf('Reflections included by the parent:\n')
fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_ind %d, hkl [%2d %2d %2d] <- matches %d, diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
[(1:numel(vis_sel))', ...
grs(1).allblobs.omega(vis_sel), grs(2).allblobs.omega(vis_sel), ...
grs(1).allblobs.eta(vis_sel), grs(2).allblobs.eta(vis_sel), ...
grs(1).allblobs.omind(vis_sel), hklsp, refl_matches(vis_sel, 1), ...
w_diffs, w_diffs < grs(1).allblobs.lorentzfac(vis_sel), ...
grs(1).allblobs.lorentzfac(vis_sel), ...
acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
]');
vis_sel = find(refl_matches(:, 1));
w_diffs = abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180);
if (size(grs(1).allblobs.hklsp, 2) == 4)
hklsp = grs(1).allblobs.hklsp(vis_sel, [1 2 4]);
else
hklsp = grs(1).allblobs.hklsp(vis_sel, :);
end
fprintf('Reflections shared by parent and first twin:\n')
fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_ind %d, hkl [%2d %2d %2d] <- matches %d, diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
[(1:numel(vis_sel))', ...
grs(1).allblobs.omega(vis_sel), grs(2).allblobs.omega(vis_sel), ...
grs(1).allblobs.eta(vis_sel), grs(2).allblobs.eta(vis_sel), ...
grs(1).allblobs.omind(vis_sel), hklsp, refl_matches(vis_sel, 1), ...
w_diffs, w_diffs < grs(1).allblobs.lorentzfac(vis_sel), ...
grs(1).allblobs.lorentzfac(vis_sel), ...
acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
]');
end
function produce_beam_intensity_renorm_figure(refgr, p, ref_ondet, ref_included, ref_selected)
try
refgr = gtComputeIncomingBeamIntensity(refgr, p);
beam_ints = refgr.beam_intensity(ref_selected);
beam_ints = beam_ints / mean(beam_ints);
catch mexc
gtPrintException(mexc, 'Skipping beam intensity renormalization')
beam_ints = 1;
end
blob_ints = refgr.intensity(ref_selected);
vis_spots = ref_ondet(ref_included(ref_selected));
thetatypes = refgr.allblobs.thetatype(vis_spots);
if (isfield(p.cryst(refgr.phaseid), 'int_exp'))
predicted_ints = p.cryst(refgr.phaseid).int_exp(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
end
predicted_ints = p.cryst(refgr.phaseid).int(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
fprintf('%d) %d, %g\n', ...
[(1:numel(blob_ints))', refgr.allblobs.thetatype(vis_spots), ...
blob_ints - predicted_ints / mean(predicted_ints) * mean(blob_ints)]')
hold(ax, 'off')
end
function s = get6DExtendedProjDefinition(num_structs)
if (~exist('num_structs', 'var') || isempty(num_structs))
num_structs = 1;
end
repl_cell = cell(num_structs, 1);
s = struct(...
'ondet', repl_cell, ...
'included', repl_cell, ...
'selected', repl_cell, ...
'allreflbool_ondet', repl_cell, ...
'allreflbool_included', repl_cell, ...
'allreflbool_selected', repl_cell, ...
'shared_refl', repl_cell, ...
'shared_refl_pos_in_parent', repl_cell, ...
'shared_refl_included_pos', repl_cell );
end
function s = find_matches_with_refor(refor, twingr, refl_matches, conf)
size_refl_list = size(refor.allblobs.omega);
refor_proj = refor.proj(conf.det_index);
ref_ondet = refor_proj.ondet;
ref_included = refor_proj.included;
ref_selected = refor_proj.selected;
bool_ref_inc = false(size_refl_list);
bool_ref_inc(ref_ondet(ref_included)) = true;
bool_ref_sel = false(size_refl_list);
bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;
twingr_proj = twingr.proj(conf.det_index);
twin_ondet = twingr_proj.ondet;
twin_inc = twingr_proj.included;
twin_sel = twingr_proj.selected;
bool_twin_ondet = false(size_refl_list);
bool_twin_ondet(twin_ondet) = true;
bool_twin_inc = false(size_refl_list);
bool_twin_inc(twin_ondet(twin_inc)) = true;
bool_twin_sel = false(size_refl_list);
bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;
% Indices of reflections that exist on the detector for the
% twin, and were included for the parent
bool_test_inc = bool_ref_inc & bool_twin_ondet;
% Alternative code to the use of refl_matches
% inc_ref_omegas = samp_ors(1).allblobs.omega(bool_test_inc);
% inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(bool_test_inc);
% inc_ref_etas = samp_ors(1).allblobs.eta(bool_test_inc);
%
% inc_twin_omegas = samp_ors(ii_g).allblobs.omega(bool_test_inc);
% inc_twin_etas = samp_ors(ii_g).allblobs.eta(bool_test_inc);
%
% w_diff = abs(mod(inc_ref_omegas - inc_twin_omegas + 180, 360) - 180);
% n_diff = abs(mod(inc_ref_etas - inc_twin_etas + 180, 360) - 180);
%
% test_shared = ...
% (w_diff < conf.angular_toll * inc_ref_Lfac) ...
% & (n_diff < conf.angular_toll);
test_shared = refl_matches(bool_test_inc);
% let's add the shared as included
shared_refl = false(size_refl_list);
bool_test_inc_indx = find(bool_test_inc);
shared_refl(bool_test_inc_indx(test_shared)) = true;
bool_twin_inc_redundant = bool_twin_inc & shared_refl;
% So if it is shared, we use the parent's blob
bool_twin_inc = bool_twin_inc | shared_refl;
% If the parent enables it, we enable it as well
bool_twin_sel = bool_twin_sel | (shared_refl & bool_ref_sel);
s = get6DExtendedProjDefinition();
% We now find the blobs in the old twin bl structure that are
% shared with the parent
s.shared_refl_included_pos = bool_twin_inc_redundant(twin_ondet(twin_inc));
twin_inc = find(bool_twin_inc(twin_ondet));
twin_sel = bool_twin_sel(twin_ondet(twin_inc));
shared_refl_pos_in_parent = bool_ref_inc & shared_refl;
s.shared_refl_pos_in_parent = shared_refl_pos_in_parent(ref_ondet(ref_included));
% We cut it down to the size of the included
s.shared_refl = shared_refl(twin_ondet(twin_inc));
s.ondet = twin_ondet;
s.included = twin_inc;
s.selected = twin_sel;
end
function blob = load_blob(img_bboxes, img_sizes, stackUSize, stackVSize, p, conf)
blob = gtFwdSimBlobDefinition();
bb = [img_bboxes(1, 1:2), img_sizes(1, 1:2)];
blob_vol = gtGetRawRoi(img_bboxes(3), img_bboxes(6), p.acq, bb, conf.det_index);
blob_vol(blob_vol < 0) = 0;
blob_vol(isnan(blob_vol)) = 0;
blob_bb = [img_bboxes(1, 1:3), img_sizes(1, :)];
% Transposing to keep the same convention as spots
blob_vol = permute(blob_vol, [2 1 3]);
blob.mbbsize = blob_bb(4:6);
blob.mbbu = [blob_bb(1), blob_bb(1) + blob_bb(4) + 1];
blob.mbbv = [blob_bb(2), blob_bb(2) + blob_bb(5) + 1];
blob.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];
blob.intm = gtPlaceSubVolume( ...
zeros(blob_size_im, 'single'), single(blob_vol), shifts_blob);
blob.bbsize = blob_size_im;
blob_bb_im = [blob_bb(1, 1:3) - shifts_blob, blob_size_im];
blob.bbuim = [blob_bb_im(1), blob_bb_im(1) + blob_bb_im(4) - 1];
blob.bbvim = [blob_bb_im(2), blob_bb_im(2) + blob_bb_im(5) - 1];
blob.bbwim = [blob_bb_im(3), blob_bb_im(3) + blob_bb_im(6) - 1];
blob_int = sum(blob.intm(blob.mask));
% blob.intm = blob.intm / blob_int;
blob.intensity = blob_int;
end
zUtil_Deformation/gt6DMergeRealSpaceVolumes.m
+ 1
− 1
View file @ 84d51cca
...@@ -35,5 +35,5 @@ function [gr_center_pix, vol_size_pix, bbox_pix] = gt6DMergeRealSpaceVolumes(grs ...@@ -35,5 +35,5 @@ function [gr_center_pix, vol_size_pix, bbox_pix] = gt6DMergeRealSpaceVolumes(grs
vol_size_pix = ceil(vol_size_pix); vol_size_pix = ceil(vol_size_pix);
vol_half_size_pix = vol_size_pix / 2; vol_half_size_pix = vol_size_pix / 2;
bbox_pix = [gr_center_pix - vol_half_size_pix, gr_center_pix + vol_half_size_pix]; bbox_pix = [floor(gr_center_pix - vol_half_size_pix), ceil(gr_center_pix + vol_half_size_pix)];
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