Skip to content
Snippets Groups Projects

Create forward prjection structure for twin cluster.

Open Create forward prjection structure for twin cluster.
2-add-TRD-fwd-proj into master
Open Zheheng Liu requested to merge 2-add-TRD-fwd-proj into master
Viewing commit fc1d0109
Next
Show latest version
2 files
+ 978
− 0
Compare changes
  • Side-by-side
  • Inline
Files
2
zUtil_Twins/gt6DCreateProjDataFromTwinClusterFwdProj.m 0 → 100755
+ 940
− 0
function [cl, estim_space_bbox_pix, estim_orient_bbox_rod] = gt6DCreateProjDataFromTwinClusterFwdProj(grs_list, phase_id, p, varargin)
% FUNCTION cl = gt6DCreateProjDataFromTwinClusterFwdProj(parent_id, variants, phase_id, varargin)
% INPUT:
% parent_id: Parent grain ID (it could be the parent grain structure)
% variants: an array of structures defined like the following:
% struct('type', {'gr_id' | 'gr_str' | 'var_num'}, 'data', {...} )
% OUTPUT:
% samp_ors: is an array of grain structures containing the parent and its
% twins
if (~exist('phase_id', 'var') || isempty(phase_id))
phase_id = 1;
end
conf = struct( ...
'verbose', false, ...
'det_index', 1, ...
'min_eta', 15, ...
'angular_tol', 1, ...
'rspace_oversize', 1, ...
'ospace_oversize', 1, ...
'check_spots', false, ...
'stack_oversize', 1.4, ...
'merge_dist_deg', 1, ...
'use_raw_images', 1, ...
'include_all', true, ...
'select_all', true, ...
'save', false, ...
'desel_global_pos', [], ...
'sel_global_pos', []);
conf = parse_pv_pairs(conf, varargin);
if (~exist('p', 'var') || isempty(p))
p = gtLoadParameters();
end
symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
if isstruct(grs_list)
grs = grs_list;
grs_list = cat(2, grs(:).id);
[grs_list, tmp_ind] = unique(grs_list);
grs = grs(tmp_ind);
else
fprintf('Loading grains: ')
grs_list = unique(grs_list(:))';
num_grains = numel(grs_list);
grs(1) = gtLoadGrain(phase_id, grs_list(1));
grs(1) = gtCalculateGrain(grs(1), p);
for ii_g = num_grains:-1:2
grs(ii_g) = gtLoadGrain(phase_id, grs_list(ii_g));
grs(ii_g) = gtCalculateGrain(grs(ii_g), p);
end
end
num_grains = numel(grs_list);
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);
refgr_omind = refgr.allblobs(conf.det_index).omind;
for ii_g = num_grains:-1:2
grs(ii_g) = gtGrainAllblobsFilterOrder(grs(ii_g), refgr_omind, conf.det_index);
end
size_refl_list = size(grs(1).allblobs(conf.det_index).omega);
if (numel(conf.use_raw_images) == 1)
conf.use_raw_images = conf.use_raw_images(ones(num_grains, 1));
end
% % % deselect the spots that are shared by the grains not belonging to
% the cluster
% Remove the other phases
spots_shared_by_externel_gr_to_desel = sfDetermineTheBlobsSharedByExternalGrainsToDeselect(p.acq(1).dir, phase_id, grs_list, grs, size_refl_list(1), num_grains, conf.det_index);
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_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo(conf.det_index), p.samgeo, false, false);
estim_space_bbox_mm = [ ...
gtGeoSam2Sam(estim_space_bbox_pix(1:3), p.recgeo(conf.det_index), p.samgeo, false, false), ...
gtGeoSam2Sam(estim_space_bbox_pix(4:6), p.recgeo(conf.det_index), p.samgeo, false, false) ];
bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3);
fprintf('\b\b: Done in %g seconds.\n', toc(c))
fprintf(' - Estimated extent: [%3g, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix);
fprintf(' BBox size: %3d, %3d, %3d (%g, %g, %g mm)\n', bbox_size_pix, bbox_size_mm);
estim_orient_bbox_rod = zeros(num_grains, 6);
bbox_size_rod = zeros(num_grains, 3);
% 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)
sampler = GtOrientationSampling(p, grs(ii_g), ...
'detector_index', conf.det_index, 'verbose', conf.verbose);
r_vecs = sampler.guess_ODF_BB()';
estim_orient_bbox_rod(ii_g, :) = [min(r_vecs, [], 1), max(r_vecs, [], 1)];
bbox_size_rod(ii_g, :) = estim_orient_bbox_rod(ii_g, 4:6) - estim_orient_bbox_rod(ii_g, 1:3);
% oversizing the orientation a bit
delta_bbox_size_rod = bbox_size_rod(ii_g, :) .* (conf.ospace_oversize - 1) / 2;
estim_orient_bbox_rod(ii_g, :) = estim_orient_bbox_rod(ii_g, :) + [-delta_bbox_size_rod, delta_bbox_size_rod];
bbox_size_rod(ii_g, :) = estim_orient_bbox_rod(ii_g, 4:6) - estim_orient_bbox_rod(ii_g, 1:3);
bbox_size_deg = 2 * atand(bbox_size_rod(ii_g, :));
fprintf(' - Estimated extent: [%g, %g, %g] -> [%g, %g, %g]\n', estim_orient_bbox_rod(ii_g, :));
fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
end
% include_all and select_all
grs = sfIncludeAllAndSelectAllAndUniformBlStruct(grs, conf);
% Determining the regions to be merged!
merge_lists = determine_merge_list(estim_orient_bbox_rod, num_grains, conf);
num_grains_new = numel(merge_lists);
new_grs_list_indices = zeros(1, num_grains_new);
num_created_grains = 0;
for ii_m = 1:num_grains_new
if (numel(merge_lists{ii_m}) > 1)
num_created_grains = num_created_grains + 1;
new_grs_list_indices(ii_m) = num_grains + num_created_grains;
else
new_grs_list_indices(ii_m) = merge_lists{ii_m};
end
end
samp_ors = cell(num_grains+num_created_grains, 1);
max_img_sizes = zeros(num_grains+num_created_grains, 2);
inconvenient_etas = cell(num_grains_new, 1);
if (num_grains_new < num_grains)
conf.use_raw_images = cellfun(@(x) max(conf.use_raw_images(x)), merge_lists);
% Merging regions
new_grs = cell(num_created_grains, 1);
for ii_m = 1:num_grains_new
if (new_grs_list_indices(ii_m) > num_grains)
new_gr = grs(merge_lists{ii_m}(end));
for ii_g = merge_lists{ii_m}(end-1:-1:1)
new_gr = gtGrainMerge(new_gr, grs(ii_g), p, [], 1 + (conf.use_raw_images(ii_m) == 1));
end
new_grs{new_grs_list_indices(ii_m) - num_grains} = new_gr;
grs(merge_lists{ii_m}) = sfUniformProj(grs(merge_lists{ii_m}), new_gr, conf.det_index, size_refl_list(1));
end
end
grs(num_grains+1:num_grains+num_created_grains) = [new_grs{:}];
estim_orient_bbox_rod = [estim_orient_bbox_rod; zeros(num_created_grains, 6)];
for ii_m = 1:num_grains_new
if (new_grs_list_indices(ii_m) > num_grains)
estim_orient_bbox_rod(new_grs_list_indices(ii_m), :) = [ ...
min(estim_orient_bbox_rod(merge_lists{ii_m}, 1:3), [], 1), ...
max(estim_orient_bbox_rod(merge_lists{ii_m}, 4:6), [], 1) ...
];
end
end
bbox_size_rod = estim_orient_bbox_rod(:, 4:6) - estim_orient_bbox_rod(:, 1:3);
fprintf('New orientation-space regions:\n')
for ii_g = 1:num_grains+num_created_grains
bbox_size_deg = 2 * atand(bbox_size_rod(ii_g, :));
fprintf(' - Estimated extent: [%g, %g, %g] -> [%g, %g, %g]\n', estim_orient_bbox_rod(ii_g, :));
fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
end
end
% Computing the new grain structures
for ii_g = 1:num_grains+num_created_grains
refor = struct( ...
'id', grs(ii_g).id, 'phaseid', grs(ii_g).phaseid, ...
'center', gr_center_mm, 'R_vector', grs(ii_g).R_vector );
refor = gtCalculateGrain(refor, p, 'ref_omind', refgr_omind);
refor.bb_ors = gt6DCalculateFullSpaceCorners(estim_space_bbox_mm, ...
bbox_size_mm, estim_orient_bbox_rod(ii_g, :), ...
bbox_size_rod(ii_g, :), refgr, p, conf.det_index);
samp_ors{ii_g} = refor;
end
samp_ors = [samp_ors{:}];
matches_ref_to_twin = cell(num_grains_new-1, 1);
matches_twin_to_ref = cell(num_grains_new-1, 1);
for ii_g = 1:num_grains_new-1
% refl_matches links the positions (i.e the hklsp type) of shared reflections
% i.e. refl_matches= [2, 0; 0, 1] indicates that the reflection in
% twin(1).allblobs(2) shares the spot with parent.allblobs(1) and
% twin(2).allblobs(1) shares the spot with parent.allblobs(1)
mis_angle_diff = zeros(1, num_grains_new-ii_g);
for ii_g_2 = (ii_g+1):num_grains_new
twin_info = gtCrystTwinTest(...
samp_ors(new_grs_list_indices(ii_g)).R_vector, ...
samp_ors(new_grs_list_indices(ii_g_2)).R_vector, phase_id);
if (twin_info.sigmaAnd > 0)
mis_angle_diff(ii_g_2-ii_g) = min(abs(twin_info.mis_angle - twin_info.theo_mis_angle));
end
end
angular_tolerance = max(conf.angular_tol, mis_angle_diff * 1.05);
if (ii_g == 1)
fprintf(' - Determining matching master reflections with other regions..')
else
angular_tolerance = 2 * angular_tolerance;
fprintf(' - Determining matching twin (%d) reflections with other regions..', ii_g)
end
[matches_ref_to_twin{ii_g}, matches_twin_to_ref{ii_g}] = gt6DGetMatchingReflections(...
samp_ors(new_grs_list_indices(ii_g)), ...
samp_ors(new_grs_list_indices((ii_g+1):end)), angular_tolerance, conf.det_index);
fprintf('\b\b: Done.\n')
if (conf.verbose > 1)
produce_matching_reflections_table(samp_ors(new_grs_list_indices(ii_g:end)), conf, angular_tolerance, matches_ref_to_twin{ii_g})
end
end
extended_projs = find_matches_with_refor(...
grs(new_grs_list_indices), matches_ref_to_twin, matches_twin_to_ref, ...
size_refl_list(1), num_grains_new, conf);
% updating grs
for ii_m = num_grains_new:-1:1
for ii_g = unique([merge_lists{ii_m}, new_grs_list_indices(ii_m)])
gr_proj = grs(ii_g).proj(conf.det_index);
gr_proj.included = extended_projs(ii_m).included;
gr_proj.selected = extended_projs(ii_m).selected;
gr_proj.bl = gtFwdSimBlobDefinition('blob', numel(gr_proj.included));
gr_proj.bl(extended_projs(ii_m).included_old_inc) = grs(ii_g).proj(conf.det_index).bl;
grs(ii_g).proj(conf.det_index) = gr_proj;
end
end
for ii_g = 1:num_grains_new
local_ondet = extended_projs(ii_g).ondet;
local_included = extended_projs(ii_g).included;
refor_ns = refor.allblobs(conf.det_index).eta(local_ondet(local_included));
% We avoid the vertical spots for convenience
inconvenient_etas{ii_g} = acosd(abs(cosd(refor_ns))) < conf.min_eta;
end
fprintf('Determining UVW bounding boxes..')
c = tic();
if any(conf.use_raw_images)
uvw_tab = cell(num_grains+num_created_grains, 1);
img_bboxes = cell(num_grains+num_created_grains, 1);
img_sizes = cell(num_grains+num_created_grains, 1);
end
for ii_m = 1:num_grains_new
if conf.use_raw_images(ii_m)
local_ondet = extended_projs(ii_m).ondet;
local_included = extended_projs(ii_m).included;
local_inc_refl = local_ondet(local_included);
for ii_g = unique([merge_lists{ii_m}, new_grs_list_indices(ii_m)])
uvw_tab{ii_g} = zeros(numel(local_inc_refl), numel(samp_ors(ii_g).bb_ors), 3);
for ii_o = 1:numel(samp_ors(ii_g).bb_ors)
uvw_tab{ii_g}(:, ii_o, :) ...
= samp_ors(ii_g).bb_ors(ii_o).allblobs(conf.det_index).detector.uvw(local_inc_refl, :);
end
% Let's treat those blobs at the w edge 360->0
% (from the sampled orientations perspective)
or_ws = samp_ors(ii_g).allblobs(conf.det_index).omega(local_inc_refl) / gtAcqGetOmegaStep(p, conf.det_index);
uvw_tab{ii_g}(:, :, 3) = gtGrainAnglesTabularFix360deg(uvw_tab{ii_g}(:, :, 3), or_ws, p);
img_bboxes{ii_g} = [ ...
reshape(floor(min(uvw_tab{ii_g}, [], 2)), [], 3), ...
reshape(ceil(max(uvw_tab{ii_g}, [], 2)), [], 3) ];
img_sizes{ii_g} = img_bboxes{ii_g}(:, 4:6) - img_bboxes{ii_g}(:, 1:3) + 1;
if conf.use_raw_images(ii_m) == 2
tmp_bbsizes = img_sizes{ii_g};
else
tmp_bbsizes = cat(1, img_sizes{ii_g}, grs(ii_g).proj(conf.det_index).bl.bbsize);
end
max_img_sizes(ii_g, :) = max(tmp_bbsizes(:, 1:2), [], 1);
end
if (conf.verbose || true)
or_ns = samp_ors(new_grs_list_indices(ii_m)).allblobs(conf.det_index).eta(local_inc_refl);
tbl_out_format = ['%02d)du %8d, dv %8d, dw %8d, eta: %7.3f, omega: %7.3f ' ...
'<- used: %d, u: [%4d %4d], v: [%4d %4d], w: [%4d %4d], sel: %d'];
tbl_output = [(1:numel(or_ns))', ...
img_sizes{new_grs_list_indices(ii_m)}, or_ns, or_ws * gtAcqGetOmegaStep(p, conf.det_index), ...
~inconvenient_etas{ii_m}, ...
img_bboxes{new_grs_list_indices(ii_m)}(:, [1 4 2 5 3 6]), ...
extended_projs(ii_m).selected ];
fprintf('\nImages for orientation %d:\n', ii_m)
if (ii_m == 1)
fprintf([tbl_out_format '\n'], ...
tbl_output');
else
pos_of_shared_blobs = zeros(size(extended_projs(ii_m).shared_refl_prevgrs_inc, 1), 2);
for ii_prevg = ii_m - 1:-1:1 % this way or loop for blobs?
pos_of_shared_blobs(extended_projs(ii_m).shared_refl_prevgrs_inc(:, ii_prevg), 1) = extended_projs(ii_m).shared_refl_pos_in_prevgrs_inc{ii_prevg};
pos_of_shared_blobs(extended_projs(ii_m).shared_refl_prevgrs_inc(:, ii_prevg), 2) = ii_prevg;
end
fprintf([tbl_out_format ', shared: %d in grain %d\n'], ...
[tbl_output(), pos_of_shared_blobs]');
end
end
else
for ii_g = unique([merge_lists{ii_m}, new_grs_list_indices(ii_m)])
tmp_bbsizes = cat(1, grs(ii_g).proj(conf.det_index).bl.bbsize);
max_img_sizes(ii_g, :) = max(tmp_bbsizes(:, 1:2), [], 1);
end
end
end
[stackUSize, stackVSize] = get_stack_UV_size(max_img_sizes, p, conf);
fprintf('\b\b: Done in %g seconds.\n', toc(c))
for ii_m = 1:num_grains_new
fprintf('Processing orientation-space region: %d/%d\n', ii_m, num_grains_new)
num_blobs = numel(extended_projs(ii_m).included);
c = tic();
switch conf.use_raw_images(ii_m)
case 2
fprintf(' - Loading raw images: ')
blobs = gtFwdSimBlobDefinition('blob', num_blobs);
for ii_b = 1:num_blobs
num_chars = fprintf('%03d/%03d', ii_b, num_blobs);
tmp_blob = blobs(ii_b * ones(1, numel(merge_lists{ii_m})));
for ii_subg = 1:numel(tmp_blob)
tmp_blob(ii_subg) = load_blob( ...
img_bboxes{merge_lists{ii_m}(ii_subg)}(ii_b, :), img_sizes{merge_lists{ii_m}(ii_subg)}(ii_b, :), ...
stackUSize, stackVSize, p, conf);
end
tmp_blob = gtGrainMergeProjBl(tmp_blob);
blobs(ii_b) = tmp_blob;
fprintf(repmat('\b', [1, num_chars]));
end
case 1
fprintf(' - Loading raw images for missing blobs: ')
blobs = grs(new_grs_list_indices(ii_m)).proj(conf.det_index).bl;
for ii_b = 1:num_blobs
num_chars = fprintf('%03d/%03d', ii_b, num_blobs);
if isempty(blobs(ii_b).intm)
for ii_subg = 1:numel(merge_lists{ii_m})
tmp_blob = grs(merge_lists{ii_m}(ii_subg)).proj(conf.det_index).bl(ii_b);
if isempty(tmp_blob.intm)
tmp_blob = load_blob( ...
img_bboxes{merge_lists{ii_m}(ii_subg)}(ii_b, :), img_sizes{merge_lists{ii_m}(ii_subg)}(ii_b, :), ...
stackUSize, stackVSize, p, conf);
end
blobs(ii_b) = gtGrainMergeProjBl([blobs(ii_b), tmp_blob]);
end
end
fprintf(repmat('\b', [1, num_chars]));
end
end
grs(new_grs_list_indices(ii_m)).proj(conf.det_index).bl = blobs;
fprintf('\b\b: Done in %g seconds.\n - Producing proj data-structure..', toc(c))
proj = gtFwdSimProjDefinition();
proj.centerpix = gr_center_pix;
vol_size = bbox_size_pix;
% 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 = extended_projs(ii_m).ondet;
proj.included = extended_projs(ii_m).included;
proj.selected = extended_projs(ii_m).selected;
if (ii_m > 1)
shared_inc = extended_projs(ii_m).shared_refl_prevgrs_inc;
shared_inc = any(shared_inc, 2);
non_shared_blobs = blobs(~shared_inc);
shared_inc = find(shared_inc);
proj.global_pos = extended_projs(ii_m).refl_pos_inc_global(proj.ondet(proj.included));
for ii_b = 1:numel(shared_inc)
all_blobs_list(proj.global_pos(shared_inc(ii_b))) = gtGrainMergeProjBl(...
[all_blobs_list(proj.global_pos(shared_inc(ii_b))), blobs(shared_inc(ii_b))]);
end
all_blobs_list = cat(1, all_blobs_list(:), non_shared_blobs(:));
else
proj.global_pos = (1:numel(blobs));
all_blobs_list = blobs;
end
samp_ors(new_grs_list_indices(ii_m)).proj(conf.det_index) = proj;
fprintf('\b\b: Done in %g seconds.\n', toc(c))
end
samp_ors = samp_ors(new_grs_list_indices);
fprintf('\n - Equalizing blob sizes: ')
c = tic;
tot_num_blobs = numel(all_blobs_list);
for ii_b = 1:tot_num_blobs
if (stackUSize < size(all_blobs_list(ii_b).intm, 1))
stackUSize = size(all_blobs_list(ii_b).intm, 1);
end
if (stackVSize < size(all_blobs_list(ii_b).intm, 2))
stackVSize = size(all_blobs_list(ii_b).intm, 2);
end
end
for ii_b = 1:tot_num_blobs
num_chars = fprintf('%02d/%02d', ii_b, tot_num_blobs);
all_blobs_list(ii_b) = equalize_blob_size(all_blobs_list(ii_b), stackUSize, stackVSize);
fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars]));
end
fprintf('\b\b: Done in %g seconds.\n', toc(c))
for ii_m = 1:num_grains_new
proj = samp_ors(ii_m).proj(conf.det_index);
proj.bl = all_blobs_list(proj.global_pos(:));
spots = arrayfun(@(x){sum(x.intm, 3)}, proj.bl);
spots = permute(cat(3, spots{:}), [1 3 2]);
proj.stack = spots;
samp_ors(ii_m).proj(conf.det_index) = proj;
if (conf.check_spots && ii_m > 1)
samp_ors(ii_m).proj.selected = gtGuiGrainMontage(samp_ors(ii_m), [], true, conf.det_index);
end
end
if (conf.verbose)
f = figure();
ax = axes('parent', f);
hold(ax, 'on');
for ii_g = 1:num_grains_new
gt6DPlotOrientationBBox(ax, cat(1, samp_ors(ii_g).bb_ors(:).R_vector));
scatter3(ax, samp_ors(ii_g).R_vector(1), samp_ors(ii_g).R_vector(2), samp_ors(ii_g).R_vector(3), 30);
end
hold(ax, 'off');
drawnow();
end
cl = struct(...
'samp_ors', samp_ors, ...
'blobs', {all_blobs_list}, ...
'conf', conf);
% % % deselect and select the spots % % %
c = tic;
fprintf('Select and deselect the spots..')
for ii_g = 1:num_grains_new
tmp_ind_desel = any(bsxfun(@eq, samp_ors(ii_g).id(:), spots_shared_by_externel_gr_to_desel.grs_list(:)'), 1);
tmp_ind_desel = any(spots_shared_by_externel_gr_to_desel.bool_deselect(:, tmp_ind_desel), 2);
tmp_ind_desel = tmp_ind_desel(samp_ors(ii_g).proj(conf.det_index).ondet(samp_ors(ii_g).proj(conf.det_index).included));
tmp_ind_desel = samp_ors(ii_g).proj(conf.det_index).global_pos(tmp_ind_desel);
spots_shared_by_externel_gr_to_desel.global_pos = unique([spots_shared_by_externel_gr_to_desel.global_pos(:); tmp_ind_desel]);
end
sel_desel_global_pos = {spots_shared_by_externel_gr_to_desel.global_pos(:), conf.sel_global_pos(:), conf.desel_global_pos(:)};
cl = gtTwinClusterFwdProjUpdateSelected(cl, conf.det_index, sel_desel_global_pos{2}, [sel_desel_global_pos{1}; sel_desel_global_pos{3}], 'flag_selected_all', false, 'deselect_eta_tol', conf.min_eta);
fprintf('\b\b: Done in %g seconds.\n', toc(c))
if (conf.save)
str_ids = sprintf('_%04d', sort(grs_list));
grain_filename = fullfile(p.acq(conf.det_index).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, '-struct', 'cl', '-v7.3');
fprintf('\b\b: Done.\n')
% fprintf('Saving to sample.mat..')
% sample = GtSample.loadFromFile();
% cl_info = GtPhase.makeCluster(grs_list, ...
% cat(1, samp_ors(:).R_vector), ...
% estim_space_bbox_pix, estim_orient_bbox_rod, 2);
% sample.phases{phase_id}.setClusterInfo(grs_list, cl_info);
% sample.saveToFile();
% fprintf('\b\b: Done.\n')
end
end
function grs = sfIncludeAllAndSelectAllAndUniformBlStruct(grs, conf)
blob_struct = gtFwdSimBlobDefinition('blob', 1);
for ii_g = numel(grs):-1:1
if (conf.use_raw_images(ii_g) && conf.include_all)
ond_ind = grs(ii_g).proj(conf.det_index).ondet;
inc_ind = grs(ii_g).proj(conf.det_index).included;
sel_bool = grs(ii_g).proj(conf.det_index).selected;
gr_proj_bl = gtFwdSimBlobDefinition('blob', numel(ond_ind));
for ii_b = 1:numel(inc_ind)
gr_proj_bl(inc_ind(ii_b)) = gtAddMatFile(blob_struct, grs(ii_g).proj(conf.det_index).bl(ii_b), true, false, false, false, false);
gr_proj_bl(inc_ind(ii_b)).nvar = [];
end
gr_selected = false(size(ond_ind));
gr_selected(inc_ind(sel_bool)) = true;
grs(ii_g).proj(conf.det_index).selected = gr_selected;
grs(ii_g).proj(conf.det_index).included = 1:numel(ond_ind);
grs(ii_g).proj(conf.det_index).bl = gr_proj_bl;
else
gr_proj_bl = gtFwdSimBlobDefinition('blob', numel(grs(ii_g).proj(conf.det_index).included));
for ii_b = 1:numel(gr_proj_bl)
gr_proj_bl(ii_b) = gtAddMatFile(blob_struct, grs(ii_g).proj(conf.det_index).bl(ii_b), true, false, false, false, false);
gr_proj_bl(ii_b).nvar = [];
end
grs(ii_g).proj(conf.det_index).bl = gr_proj_bl;
end
end
if conf.select_all
for ii_g = numel(grs):-1:1
grs(ii_g).proj(conf.det_index).selected = true(size(grs(ii_g).proj(conf.det_index).included));
end
end
end
function spots_shared_by_externel_gr_to_desel = sfDetermineTheBlobsSharedByExternalGrainsToDeselect(dir, phase_id, grs_list, grs, height_refl_list, num_grains, det_ind)
% % % deselect the spots that are shared by the grains not belonging to
% the cluster
% Remove the other phases
spot2grainfilename = fullfile(dir,'4_grains','spot2grain.mat');
spots_shared_by_externel_gr_to_desel = struct('grs_list', grs_list, ...
'bool_deselect', false(height_refl_list, num_grains), ...
'global_pos', []);
if exist(spot2grainfilename, 'file')
c = tic;
fprintf('Find the spots shared by the grain not belonging to the cluster..')
spot2grain = load(spot2grainfilename);
for ii_s = 1:numel(spot2grain.spot2phase)
tmp_ind_spot2grain = (spot2grain.spot2phase{ii_s} ~= phase_id);
if ~isempty(tmp_ind_spot2grain)
spot2grain.spot2phase{ii_s}(tmp_ind_spot2grain) = [];
spot2grain.spot2grain{ii_s}(tmp_ind_spot2grain) = [];
end
end
spot2grain = spot2grain.spot2grain;
% find the spots to deselect
for ii_g = num_grains:-1:1
inc_ind = grs(ii_g).proj(det_ind).ondet(grs(ii_g).proj(det_ind).included);
spot2grain_i = spot2grain(grs(ii_g).fwd_sim(det_ind).difspotID(:));
for ii_spot = 1:numel(spot2grain_i)
check_existing_other_grain = any(~any(bsxfun(@eq, spot2grain_i{ii_spot}(:), grs_list(:)'), 2));
if check_existing_other_grain
spots_shared_by_externel_gr_to_desel.bool_deselect(inc_ind(ii_spot), ii_g) = true;
end
end
end
fprintf('\b\b: Done in %g seconds. %d spots found.\n', toc(c), sum(spots_shared_by_externel_gr_to_desel.bool_deselect(:)))
end
end
function grs = sfUniformProj(grs, refgr, det_ind, totnum)
for ii_g = numel(grs):-1:1
proj = grs(ii_g).proj(det_ind);
proj.bl = gtFwdSimBlobDefinition('blob', numel(refgr.proj(det_ind).included));
bool_inc = false(totnum, 1);
bool_inc(proj.ondet(proj.included)) = true;
bool_inc_ondet = false(totnum, 1);
bool_inc_ondet(refgr.proj(det_ind).ondet) = true;
bool_inc_ondet = bool_inc_ondet & bool_inc;
bool_included_old = bool_inc(refgr.proj(det_ind).ondet(refgr.proj(det_ind).included));
proj.bl(bool_included_old) = grs(ii_g).proj(det_ind).bl(bool_inc_ondet(proj.ondet(proj.included)));
proj.ondet = refgr.proj(det_ind).ondet;
proj.included = refgr.proj(det_ind).included;
proj.selected = refgr.proj(det_ind).selected;
grs(ii_g).proj(det_ind) = proj;
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))];
if any(conf.use_raw_images)
% Because in this case we saturate the volume field of view from
% the spots already, and we only make up for what will probably be
% used as oversize later on for the volume
stack_imgs_oversize = p.rec.grains.options.rspace_oversize;
else
stack_imgs_oversize = min(p.fsim.oversize, conf.stack_oversize);
end
u_size = round(max_img_sizes(1) * stack_imgs_oversize);
v_size = round(max_img_sizes(2) * stack_imgs_oversize);
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
function blob = equalize_blob_size(blob, stackUSize, stackVSize)
new_bbsize = [stackUSize, stackVSize, blob.bbsize(3)];
shift = floor((new_bbsize - blob.bbsize) / 2);
new_bbuim = blob.bbuim - shift(1);
new_bbuim(2) = new_bbuim(1) + new_bbsize(1) - 1;
new_bbvim = blob.bbvim - shift(2);
new_bbvim(2) = new_bbvim(1) + new_bbsize(2) - 1;
blob.bbuim = new_bbuim;
blob.bbvim = new_bbvim;
blob.bbsize = new_bbsize;
blob.intm = gtPlaceSubVolume(zeros(new_bbsize, 'like', blob.intm), blob.intm, shift);
blob.mask = logical(gtPlaceSubVolume(zeros(new_bbsize, 'uint8'), blob.mask, shift));
end
function produce_matching_reflections_table(grs, conf, tols, refl_matches)
for ii = 1 : size(refl_matches, 2)
vis_sel = find(refl_matches(:, ii));
w_diffs = abs(mod(grs(1).allblobs(conf.det_index).omega(vis_sel) - grs(ii+1).allblobs(conf.det_index).omega(refl_matches(vis_sel, ii)) + 180, 360) - 180);
n_diffs = abs(mod(grs(1).allblobs(conf.det_index).eta(vis_sel) - grs(ii+1).allblobs(conf.det_index).eta(refl_matches(vis_sel, ii)) + 180, 360) - 180);
% pre-allocate
hklsp = zeros([size(grs(1).allblobs(conf.det_index).hklsp, 1), 3, numel(grs)]);
if (size(grs(1).allblobs(conf.det_index).hklsp, 2) == 4)
for jj = 1 : numel(grs)
hklsp(:,:,jj) = grs(jj).allblobs(conf.det_index).hklsp(:, [1 2 4]);
end
else
for jj = 1 : numel(grs)
hklsp(:,:,jj) = grs(jj).allblobs(conf.det_index).hklsp(:, :);
end
end
fprintf('Reflections shared by reference grain and twin number %d :\n', ii)
fprintf(['%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_ind %d, hkl [%2d %2d %2d] ' ...
'<- matches w_ind %2d, [%2d %2d %2d], diff: %5.2f, %5.2f (!! %d, L %g) [bad n %d]\n'], ...
[(1:numel(vis_sel))', ...
grs(1).allblobs(conf.det_index).omega(vis_sel), grs(ii+1).allblobs(conf.det_index).omega(refl_matches(vis_sel, ii)), ...
grs(1).allblobs(conf.det_index).eta(vis_sel), grs(ii+1).allblobs(conf.det_index).eta(refl_matches(vis_sel, ii)), ...
grs(1).allblobs(conf.det_index).omind(vis_sel), hklsp(vis_sel, :, 1), refl_matches(vis_sel, ii), ...
hklsp(refl_matches(vis_sel, ii), :, ii+1), ...
w_diffs, n_diffs, ...
n_diffs < tols(ii) & w_diffs < tols(ii) * grs(1).allblobs(conf.det_index).lorentzfac(vis_sel), ...
grs(1).allblobs(conf.det_index).lorentzfac(vis_sel), ...
acosd(abs(cosd(grs(1).allblobs(conf.det_index).eta(vis_sel)))) < conf.min_eta, ...
]');
end
end
function eproj = get6DExtendedProjDefinition(num_structs)
if (~exist('num_structs', 'var') || isempty(num_structs))
num_structs = 1;
end
repl_cell = cell(num_structs, 1);
eproj = struct(...
'ondet', repl_cell, ...
'included', repl_cell, ...
'selected', repl_cell, ...
'allreflbool_ondet', repl_cell, ...
'allreflbool_included', repl_cell, ...
'allreflbool_selected', repl_cell, ...
'included_old_to_keep', repl_cell, ...
'included_replaced', repl_cell, ...
'included_old_inc', repl_cell, ...
'shared_refl_prevgrs_inc', repl_cell, ...
'shared_refl_nextgrs_inc', repl_cell, ...
'shared_refl_pos_in_prevgrs_inc', repl_cell, ...
'shared_refl_pos_in_nextgrs_inc', repl_cell, ...
'refl_pos_inc_global', repl_cell);
end
function eproj = find_matches_with_refor(grs, matches_ref_to_twin, matches_twin_to_ref, length_allbls_refl_list, num_grains_new, conf)
eproj = get6DExtendedProjDefinition(num_grains_new);
for ii_g = 1:num_grains_new
gr_proj = grs(ii_g).proj(conf.det_index);
gr_ondet = gr_proj.ondet;
gr_included = gr_proj.included;
gr_selected = gr_proj.selected;
gr_allreflbool_ondet = false(length_allbls_refl_list, 1);
gr_allreflbool_included = false(length_allbls_refl_list, 1);
gr_allreflbool_selected = false(length_allbls_refl_list, 1);
gr_allreflbool_ondet(gr_ondet) = true;
gr_allreflbool_included(gr_ondet(gr_included)) = true;
gr_allreflbool_selected(gr_ondet(gr_included(gr_selected))) = true;
eproj(ii_g).ondet = gr_ondet;
eproj(ii_g).included = gr_included;
eproj(ii_g).selected = gr_selected;
eproj(ii_g).allreflbool_ondet = gr_allreflbool_ondet;
eproj(ii_g).allreflbool_included = gr_allreflbool_included;
eproj(ii_g).allreflbool_selected = gr_allreflbool_selected;
end
matches_othergrs_to_twin = zeros(length_allbls_refl_list, num_grains_new - 1);
matches_twin_to_othergrs = zeros(length_allbls_refl_list, num_grains_new - 1);
for ii_g = 1:num_grains_new
for ii_otherg = 1:ii_g - 1
matches_othergrs_to_twin(:, ii_otherg) = matches_ref_to_twin{ii_otherg}(:, ii_g - ii_otherg);
matches_twin_to_othergrs(:, ii_otherg) = matches_twin_to_ref{ii_otherg}(:, ii_g - ii_otherg);
end
if ii_g < num_grains_new
matches_othergrs_to_twin(:, ii_g:num_grains_new - 1) = matches_twin_to_ref{ii_g};
matches_twin_to_othergrs(:, ii_g:num_grains_new - 1) = matches_ref_to_twin{ii_g};
end
othergrs_eprojs = [eproj(1:ii_g - 1); eproj(ii_g + 1:num_grains_new)];
othergrs_ondet = cat(2, {othergrs_eprojs.ondet});
othergrs_included = cat(2, {othergrs_eprojs.included});
% bool_* vectors are all of length numel(gr.allblobs(conf.det_index).omega)
bool_othergrs_inc = cat(2, zeros(numel(eproj(ii_g).allreflbool_included), 0), othergrs_eprojs.allreflbool_included);
bool_othergrs_sel = cat(2, zeros(numel(eproj(ii_g).allreflbool_selected), 0), othergrs_eprojs.allreflbool_selected);
% list of twin reflections for which a shared spot has been predicted
bool_twin_shared = false(length_allbls_refl_list, num_grains_new - 1);
bool_twin_shared(matches_twin_to_othergrs > 0) = true;
twin_ondet = eproj(ii_g).ondet;
twin_inc = eproj(ii_g).included;
twin_sel = eproj(ii_g).selected;
% Conversion of original ondet, included, and shared to bool list
bool_twin_ondet = false(length_allbls_refl_list, 1);
bool_twin_ondet(twin_ondet) = true;
bool_twin_inc = false(length_allbls_refl_list, 1);
bool_twin_inc(twin_ondet(twin_inc)) = true;
bool_twin_sel = false(length_allbls_refl_list, 1);
bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;
% Other grains' indices of shared reflections that are ondet for the twin
othergrs_shared_ondet = subsubfunc_bool_arrays_indices_mapping(...
matches_twin_to_othergrs, bsxfun(@and, bool_twin_shared, bool_twin_ondet), ...
[length_allbls_refl_list, num_grains_new - 1]);
% ...and which were included for the previous grains
othergrs_shared_ondet_inc = bool_othergrs_inc & othergrs_shared_ondet;
othergrs_shared_selected = othergrs_shared_ondet_inc & bool_othergrs_sel;
% translate this back into twin indices
twin_shared_othergrs_inc = subsubfunc_bool_arrays_indices_mapping(...
matches_othergrs_to_twin, othergrs_shared_ondet_inc, ...
[length_allbls_refl_list, num_grains_new - 1]);
% % %
% % % % These will tell which ones to keep from the fwd-simulation
% % % bool_twin_inc_keep = bsxfun(@and, bool_twin_inc, ~any(twin_shared_othergrs_inc(:, 1:ii_g - 1), 2));
% % % eproj(ii_g).included_old_to_keep = bool_twin_inc_keep(twin_ondet(twin_inc));
% These wll be the ones that will be included from the other grains
bool_twin_inc_old = bool_twin_inc;
bool_twin_inc_import = bsxfun(@and, bool_twin_ondet, twin_shared_othergrs_inc);
bool_twin_inc(any(bool_twin_inc_import, 2)) = true;
twin_inc = find(bool_twin_inc(twin_ondet));
% These are the ones that were included
eproj(ii_g).included_old_inc = bool_twin_inc_old(bool_twin_inc);
% If a shared & twin included spot is selected in the other grains, select
% it for the twin, as well
for ii_otherg = 1:num_grains_new - 1
othergr_sel_mapping_to_twin = matches_othergrs_to_twin(othergrs_shared_selected(:, ii_otherg), ii_otherg);
bool_twin_sel(othergr_sel_mapping_to_twin) = bool_twin_sel(othergr_sel_mapping_to_twin) | bool_twin_inc(othergr_sel_mapping_to_twin);
end
twin_sel = bool_twin_sel(twin_ondet(twin_inc));
% We now find the blobs in the old twin bl structure that are
% shared with the previous grains
eproj(ii_g).shared_refl_prevgrs_inc = twin_shared_othergrs_inc(twin_ondet(twin_inc), 1:ii_g - 1);
eproj(ii_g).shared_refl_nextgrs_inc = twin_shared_othergrs_inc(twin_ondet(twin_inc), ii_g:num_grains_new - 1);
% We have to find the included local and global IDs of the corresponding
% shared and non shared reflections
twin_refl_pos_inc_global = zeros(length_allbls_refl_list, 1);
twin_shared_refl_pos_in_prevgrs_inc = cell(ii_g - 1, 1);
twin_shared_refl_pos_in_nextgrs_inc = cell(num_grains_new - ii_g, 1);
if ii_g == 1
twin_refl_pos_inc_global(twin_ondet(twin_inc)) = (1:numel(twin_inc));
twin_index_of_last_non_shared_blob = numel(twin_inc);
else
twin_non_shared_refl_prevgrs_inc = ~any(eproj(ii_g).shared_refl_prevgrs_inc, 2);
twin_refl_pos_inc_global(twin_ondet(twin_inc(twin_non_shared_refl_prevgrs_inc))) = twin_index_of_last_non_shared_blob + (1:sum(twin_non_shared_refl_prevgrs_inc));
twin_index_of_last_non_shared_blob = twin_index_of_last_non_shared_blob + sum(twin_non_shared_refl_prevgrs_inc);
for ii_prevg = 1:ii_g - 1
twin_inc_import_mapping_to_othergr = matches_twin_to_othergrs(bool_twin_inc_import(:, ii_prevg), ii_prevg);
twin_refl_pos_inc_global(bool_twin_inc_import(:, ii_prevg)) = eproj(ii_prevg).refl_pos_inc_global(twin_inc_import_mapping_to_othergr);
refl_pos_in_othergr_inc = zeros(length_allbls_refl_list, 1);
refl_pos_in_othergr_inc(othergrs_ondet{ii_prevg}(othergrs_included{ii_prevg})) = (1:numel(othergrs_included{ii_prevg}));
twin_shared_refl_pos_in_prevgrs_inc{ii_prevg} = refl_pos_in_othergr_inc(twin_inc_import_mapping_to_othergr);
end
end
for ii_nextg = 1:num_grains_new - ii_g
ii_nextg_plus_ii_g_minus_one = ii_nextg + ii_g - 1;
twin_inc_import_mapping_to_othergr = matches_twin_to_othergrs(bool_twin_inc_import(:, ii_nextg_plus_ii_g_minus_one), ii_nextg_plus_ii_g_minus_one);
refl_pos_in_othergr_inc = zeros(length_allbls_refl_list, 1);
refl_pos_in_othergr_inc(othergrs_ondet{ii_nextg_plus_ii_g_minus_one}(othergrs_included{ii_nextg_plus_ii_g_minus_one})) = (1:numel(othergrs_included{ii_nextg_plus_ii_g_minus_one}));
twin_shared_refl_pos_in_nextgrs_inc{ii_nextg} = refl_pos_in_othergr_inc(twin_inc_import_mapping_to_othergr);
end
eproj(ii_g).refl_pos_inc_global = twin_refl_pos_inc_global;
eproj(ii_g).shared_refl_pos_in_prevgrs_inc = twin_shared_refl_pos_in_prevgrs_inc;
eproj(ii_g).shared_refl_pos_in_nextgrs_inc = twin_shared_refl_pos_in_nextgrs_inc;
eproj(ii_g).included = twin_inc;
eproj(ii_g).selected = twin_sel;
eproj(ii_g).allreflbool_included = bool_twin_inc;
eproj(ii_g).allreflbool_selected = bool_twin_sel;
end
end
function out = subsubfunc_bool_arrays_indices_mapping(indices_mapping_matrix, input_bool_arrays, arrays_size)
out = false(arrays_size);
for ii_otherg = 1:arrays_size(2)
out(indices_mapping_matrix(input_bool_arrays(:, ii_otherg), ii_otherg), ii_otherg) = true;
end
end
% function [problematic_bls, t2t_shr_bls] = detect_twin_only_shared_reflections(matches_ref_to_twin, matches_twin_to_ref, ii_g)
% % Suboptimal but it allows at least to detect and deselect them
%
% % The tables correspond only to the twins (so they have one fewer entry
% % than the number of grains), and inside the tables again, we only have
% % indices from the next grain and up. This means that the tables shrink
% % for successive grains, and the indeces have to adust for that
% shared_with_parent = matches_twin_to_ref{1}(:, ii_g-1) > 0;
%
% % There will be a 0 column (which corresponds to the considered twin
% % itself)
% t2t_shr_bls = false(size(shared_with_parent, 1), ii_g-1);
% for ii_g_prev = 2:ii_g-1
% shared_with_prev_twin = matches_twin_to_ref{ii_g_prev}(:, ii_g-ii_g_prev);
% t2t_shr_bls(:, ii_g_prev) = shared_with_prev_twin;
% end
%
% if (ii_g <= numel(matches_ref_to_twin))
% shared_with_next_twins = matches_ref_to_twin{ii_g};
% t2t_shr_bls = [t2t_shr_bls, shared_with_next_twins];
% end
%
% problematic_bls = any(t2t_shr_bls, 2) & ~shared_with_parent;
% end
function blob = load_blob(img_bboxes, img_sizes, stackUSize, stackVSize, p, conf)
blob = gtFwdSimBlobDefinition('blob');
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.mask = gtPlaceSubVolume( ...
false(blob_size_im), true(size(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
function merge_lists = determine_merge_list(estim_orient_bbox_rod, num_grains, conf)
merge_lists = cell(num_grains, 1);
bbsizes = (estim_orient_bbox_rod(:, 4:6) - estim_orient_bbox_rod(:, 1:3)) / 2;
centers = (estim_orient_bbox_rod(:, 4:6) + estim_orient_bbox_rod(:, 1:3)) / 2;
% Find overlapping regions (excluding master)
for ii_g = 1:num_grains
gr_inds = ii_g+1:num_grains;
overlaps = false(num_grains-ii_g, 1);
for ii_g_2 = gr_inds
distances = abs(centers(ii_g, :) - centers(ii_g_2, :)) ...
- (bbsizes(ii_g, :) + bbsizes(ii_g_2, :));
overlaps(ii_g_2-ii_g) = all(2 * atand(distances) <= conf.merge_dist_deg);
end
if (any(overlaps))
fprintf('Orientation box: %d overlaps with orientation boxes:%s. They will be merged...\n', ...
ii_g, sprintf(' %d', gr_inds(overlaps)));
merge_lists{ii_g} = gr_inds(overlaps);
end
end
merge_table = false(num_grains, num_grains);
for ii_g = 1:num_grains
merge_table([ii_g, merge_lists{ii_g}], [ii_g, merge_lists{ii_g}]) = true;
end
to_keep = true(num_grains, 1);
for ii_g = 1:num_grains
if to_keep(ii_g)
flag_incomplete_merge = true;
while(flag_incomplete_merge)
sub_table = any(merge_table(:, [ii_g, merge_lists{ii_g}]), 2);
sub_table([ii_g, merge_lists{ii_g}]) = false;
if any(sub_table)
merge_lists{ii_g} = [find(sub_table'), merge_lists{ii_g}];
merge_table([ii_g, merge_lists{ii_g}], [ii_g, merge_lists{ii_g}]) = true;
else
flag_incomplete_merge = false;
to_keep(merge_lists{ii_g}) = false;
merge_lists{ii_g} = [ii_g, sort(merge_lists{ii_g}, 'ascend')];
end
end
end
end
merge_lists = merge_lists(to_keep);
end