From fc1d01092b3983ec102cf87dd82c235a325b4d2b Mon Sep 17 00:00:00 2001
From: Zheheng Liu <zheheng.liu@esrf.fr>
Date: Thu, 30 Jan 2025 17:50:43 +0100
Subject: [PATCH 1/3] create forward prjections for twin cluster.
Signed-off-by: Zheheng Liu <zheheng.liu@esrf.fr>
---
...gt6DCreateProjDataFromTwinClusterFwdProj.m | 940 ++++++++++++++++++
.../gtTwinClusterFwdProjUpdateSelected.m | 38 +
2 files changed, 978 insertions(+)
create mode 100755 zUtil_Twins/gt6DCreateProjDataFromTwinClusterFwdProj.m
create mode 100755 zUtil_Twins/gtTwinClusterFwdProjUpdateSelected.m
diff --git a/zUtil_Twins/gt6DCreateProjDataFromTwinClusterFwdProj.m b/zUtil_Twins/gt6DCreateProjDataFromTwinClusterFwdProj.m
new file mode 100755
index 00000000..5d87dc93
--- /dev/null
+++ b/zUtil_Twins/gt6DCreateProjDataFromTwinClusterFwdProj.m
@@ -0,0 +1,940 @@
+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
diff --git a/zUtil_Twins/gtTwinClusterFwdProjUpdateSelected.m b/zUtil_Twins/gtTwinClusterFwdProjUpdateSelected.m
new file mode 100755
index 00000000..b94e0535
--- /dev/null
+++ b/zUtil_Twins/gtTwinClusterFwdProjUpdateSelected.m
@@ -0,0 +1,38 @@
+function cl = gtTwinClusterFwdProjUpdateSelected(cl, det_ind, global_pos_to_selected, global_pos_to_deselected, varargin)
+ if (nargin < 2 || isempty(det_ind))
+ det_ind = 1;
+ end
+ if (nargin < 3)
+ global_pos_to_selected = [];
+ end
+ if (nargin < 4)
+ global_pos_to_deselected = [];
+ end
+ conf = struct(...
+ 'flag_selected_all', false, ...
+ 'deselect_eta_tol', 0);
+ conf = parse_pv_pairs(conf, varargin);
+ if conf.flag_selected_all
+ for ii_g = 1:numel(cl.samp_ors)
+ cl.samp_ors(ii_g).proj(det_ind).selected = true(size(cl.samp_ors(ii_g).proj(det_ind).selected));
+ end
+ end
+ if conf.deselect_eta_tol
+ for ii_g = 1:numel(cl.samp_ors)
+ tmp_proj = cl.samp_ors(ii_g).proj(det_ind);
+ tmp_eta = cl.samp_ors(ii_g).allblobs.eta(tmp_proj.ondet(tmp_proj.included));
+ tmp_eta = 90 - abs(abs(tmp_eta - 180) - 90);
+ tmp_ind = (tmp_eta < conf.deselect_eta_tol);
+ tmp_ind = tmp_proj.global_pos(tmp_ind);
+ global_pos_to_deselected = [global_pos_to_deselected(:); tmp_ind(:)];
+ end
+ end
+ global_pos_to_selected = unique(global_pos_to_selected);
+ global_pos_to_deselected = unique(global_pos_to_deselected);
+ for ii_g = 1:numel(cl.samp_ors)
+ ind_bool = any(bsxfun(@eq, cl.samp_ors(ii_g).proj(det_ind).global_pos(:), global_pos_to_selected(:)'), 2);
+ cl.samp_ors(ii_g).proj(det_ind).selected(ind_bool) = true;
+ ind_bool = any(bsxfun(@eq, cl.samp_ors(ii_g).proj(det_ind).global_pos(:), global_pos_to_deselected(:)'), 2);
+ cl.samp_ors(ii_g).proj(det_ind).selected(ind_bool) = false;
+ end
+end
\ No newline at end of file
--
GitLab
From 73b5cb463e3655066ddf6f1c3f33a85bb6544c0d Mon Sep 17 00:00:00 2001
From: Zheheng Liu <zheheng.liu@esrf.fr>
Date: Thu, 30 Jan 2025 17:51:03 +0100
Subject: [PATCH 2/3] some changes for gtGrainMerge
Signed-off-by: Zheheng Liu <zheheng.liu@esrf.fr>
---
4_grains/gtGrainMerge.m | 18 ++++++++-----
4_grains/gtGrainMergeProjBl.m | 49 ++++++++++++-----------------------
2 files changed, 28 insertions(+), 39 deletions(-)
diff --git a/4_grains/gtGrainMerge.m b/4_grains/gtGrainMerge.m
index 9723d20c..5d1fc0be 100644
--- a/4_grains/gtGrainMerge.m
+++ b/4_grains/gtGrainMerge.m
@@ -79,15 +79,19 @@ function grn = gtGrainMerge(gr1, gr2, p, mode, flag_merge_bl)
gr2_bl_to_merge = gr2.proj(ii_d).bl(merge_blob_ind_on_gr2_bl_list);
if ~isempty(gr1_bl_to_merge) && ~isempty(gr2_bl_to_merge)
tmp_bl_merge = gtGrainMergeProjBl([gr1_bl_to_merge, gr2_bl_to_merge]);
- elseif ~isempty(gr1_bl_to_merge)
- tmp_bl_merge = gr1_bl_to_merge;
- elseif ~isempty(gr2_bl_to_merge)
- tmp_bl_merge = gr2_bl_to_merge;
+ elseif flag_merge_bl == 1
+ if ~isempty(gr1_bl_to_merge)
+ tmp_bl_merge = gr1_bl_to_merge;
+ elseif ~isempty(gr2_bl_to_merge)
+ tmp_bl_merge = gr2_bl_to_merge;
+ else
+ error('gtGrainMerge:missing_blob', ...
+ 'Both gr1.proj(%d).bl(%d) and gr2.proj(%d).bl(%d) are empty', ii_d, merge_blob_ind_on_gr1_bl_list, ii_d, merge_blob_ind_on_gr2_bl_list)
+ end
else
- error('gtGrainMerge:missing_blob', ...
- 'Both gr1.proj(%d).bl(%d) and gr2.proj(%d).bl(%d) are empty', ii_d, merge_blob_ind_on_gr1_bl_list, ii_d, merge_blob_ind_on_gr2_bl_list)
+ tmp_bl_merge = grn.proj(ii_d).bl(ii_merge_bl);
end
- grn.proj(ii_d).bl(ii_merge_bl) = tmp_bl_merge;
+ grn.proj(ii_d).bl(ii_merge_bl) = gtAddMatFile(grn.proj(ii_d).bl(ii_merge_bl), tmp_bl_merge, true, false, false, false, false);
end
end
diff --git a/4_grains/gtGrainMergeProjBl.m b/4_grains/gtGrainMergeProjBl.m
index 348f40b0..70de07eb 100644
--- a/4_grains/gtGrainMergeProjBl.m
+++ b/4_grains/gtGrainMergeProjBl.m
@@ -14,7 +14,11 @@ function bl_merge = gtGrainMergeProjBl(bl_array, difspotID_array)
end
end
% Initialize the merged bl
- bl_merge = bl_array(1);
+ if isempty(bl_array)
+ bl_merge = gtFwdSimBlobDefinition('blob');
+ else
+ bl_merge = bl_array(1);
+ end
if numel(bl_array(:)) > 1
% Merge bounds and sizes
@@ -39,38 +43,19 @@ function bl_merge = gtGrainMergeProjBl(bl_array, difspotID_array)
if exist('difspotID_array','var') && numel(difspotID_array(:)) == numel(bl_array(:))
[~, tmp_ind] = unique(difspotID_array(:));
bl_array = bl_array(tmp_ind);
-
- % Merge intensity
- intensity_array = cat(1, bl_array(:).intensity);
- tmp_intensity = sum(intensity_array);
- intensity_array = intensity_array/tmp_intensity;
-
- % Merge mask and intm
- tmp_mask = false(bl_merge.bbsize);
- tmp_intm = single(zeros(bl_merge.bbsize));
- for i_bl = 1:numel(bl_array(:))
- tmp_vol_ind = sfVolInd(bound(:, :, i_bl), bound_merge, bl_merge.bbsize);
- tmp_mask(tmp_vol_ind) = tmp_mask(tmp_vol_ind) + bl_array(i_bl).mask(:);
- tmp_intm(tmp_vol_ind) = tmp_intm(tmp_vol_ind) + intensity_array(i_bl) * bl_array(i_bl).intm(:);
- end
- else
- % Merge mask, intm and intensity.
- intensity_array = cat(1, bl_array(:).intensity);
- tmp_mask = false(bl_merge.bbsize);
- tmp_intm = single(zeros(bl_merge.bbsize));
- for i_bl = 1:numel(bl_array(:))
- tmp_vol_ind = sfVolInd(bound(:, :, i_bl), bound_merge, bl_merge.bbsize);
- % Ignore the bl whose mask overlaps the ones of previous bls
- if any(bl_array(i_bl).mask(:) & tmp_mask(tmp_vol_ind))
- intensity_array(i_bl) = 0;
- else
- tmp_mask(tmp_vol_ind) = tmp_mask(tmp_vol_ind) + bl_array(i_bl).mask(:);
- tmp_intm(tmp_vol_ind) = tmp_intm(tmp_vol_ind) + intensity_array(i_bl) * bl_array(i_bl).intm(:);
- end
- end
- tmp_intensity = sum(intensity_array(:));
- tmp_intm = tmp_intm/tmp_intensity;
end
+ % Merge mask, intm and intensity.
+ tmp_mask = false(bl_merge.bbsize);
+ tmp_intm = single(zeros(bl_merge.bbsize));
+ for i_bl = 1:numel(bl_array(:))
+ tmp_vol_ind = sfVolInd(bound(:, :, i_bl), bound_merge, bl_merge.bbsize);
+ % Remove the overlaps
+ tmp_not_overlapped_vol = (bl_array(i_bl).mask(:) & ~tmp_mask(tmp_vol_ind)) .* bl_array(i_bl).intm(:) .* bl_array(i_bl).intensity;
+ tmp_mask(tmp_vol_ind) = tmp_mask(tmp_vol_ind) | bl_array(i_bl).mask(:);
+ tmp_intm(tmp_vol_ind) = tmp_intm(tmp_vol_ind) + tmp_not_overlapped_vol;
+ end
+ tmp_intensity = sum(tmp_intm(:));
+ tmp_intm = tmp_intm/tmp_intensity;
bl_merge.intensity = tmp_intensity;
bl_merge.mask = tmp_mask;
--
GitLab
From eaced50124b7924b1777915c46da3e969b881b1b Mon Sep 17 00:00:00 2001
From: Zheheng Liu <zheheng.liu@esrf.fr>
Date: Thu, 30 Jan 2025 18:23:02 +0100
Subject: [PATCH 3/3] Add noise variances to cluster forward projection
structure
Signed-off-by: Zheheng Liu <zheheng.liu@esrf.fr>
---
zUtil_Twins/GtGrainClusterCalcBlobNoiseVars.m | 540 ++++++++++++++++++
1 file changed, 540 insertions(+)
create mode 100644 zUtil_Twins/GtGrainClusterCalcBlobNoiseVars.m
diff --git a/zUtil_Twins/GtGrainClusterCalcBlobNoiseVars.m b/zUtil_Twins/GtGrainClusterCalcBlobNoiseVars.m
new file mode 100644
index 00000000..9c7bb010
--- /dev/null
+++ b/zUtil_Twins/GtGrainClusterCalcBlobNoiseVars.m
@@ -0,0 +1,540 @@
+classdef GtGrainClusterCalcBlobNoiseVars < handle
+ properties
+ cluster;
+ parameters;
+ invalid_pixel_nvar = 3e38;
+ detector_margin = 1;
+ use_detector_mask = [];
+ end
+
+ properties (Access = protected)
+ bbuims = [];
+ bbvims = [];
+ bbwims = [];
+ data_type;
+ det_ind;
+ detector_mask;
+ end
+
+ methods (Access = public)
+ function [self, cl] = GtGrainClusterCalcBlobNoiseVars(cl_ids, phase_id, p, detinds, varargin)
+ % GtGrainClusterCalcBlobNoiseVars(cl, phase_id, p, detinds)
+ % The function reads the blobs as the estimate of
+ % Poisson noise variance. Next, upon the parameter
+ % detector_weights.nvar_images, the background images are read
+ % to get more accurate Poisson noise. An alternative method is
+ % to load raw images as noise variances. The parameters are
+ % located at parameters.rec.grains.options.detector_weights.
+ % The optional input 'detector_margin' will inactivate the
+ % pixels at the edges of detector and default number is 1; the
+ % default value of 'use_detector_mask' is true, which uses the
+ % detector mask saved in 'detector%d_mask.edf' (same convention
+ % as dark.edf).
+ % If we have got the object det_weights =
+ % GtGrainClusterCalcBlobNoiseVars(gr, phase_id, p, detinds_1),
+ % we can read the noise variance again by
+ % det_weights.loadNoiseVarianceForIndicatedDet(detinds_2)
+
+ conf = struct(...
+ 'detector_margin', 1, ...
+ 'use_detector_mask', true);
+ conf = parse_pv_pairs(conf, varargin);
+
+ self.detector_margin = conf.detector_margin;
+ self.use_detector_mask = conf.use_detector_mask;
+ if (~exist('phase_id', 'var') || isempty(phase_id))
+ phase_id = 1;
+ end
+ if (~exist('p', 'var') || isempty(p))
+ self.parameters = gtLoadParameters;
+ else
+ self.parameters = p;
+ end
+ if isnumeric(cl_ids)
+ cl_ids = sort(unique(cl_ids), 'ascend');
+ self.cluster = gtLoadCluster(phase_id, cl_ids);
+ else
+ self.cluster = cl_ids;
+ cl_ids = sort(unique(cat(2, cl_ids.samp_ors.id)), 'ascend');
+ end
+ self.set_unused_bl_to_empty();
+ if isstruct(self.cluster.blobs)
+ self.cluster.blobs = {self.cluster.blobs};
+ end
+ self.normalize_blobs();
+ if (~exist('detinds', 'var') || isempty(detinds))
+ detinds = 1:numel(self.cluster.blobs);
+ end
+ self.loadNoiseVarianceForIndicatedDet(detinds)
+ if nargout > 1
+ cl = self.cluster;
+ end
+ end
+
+ function gr = loadNoiseVarianceForIndicatedDet(self, detinds)
+ % this function is used to read the noise variance for
+ % indicated detector indices.
+ % For example, as we have got the object det_weights =
+ % GtGrainClusterCalcBlobNoiseVars(gr, phase_id, p, detinds_1),
+ % we can read the noise variance again by
+ % det_weights.loadNoiseVarianceForIndicatedDet(detinds_2)
+ for ii_d = 1:numel(detinds)
+ if ~isempty(self.cluster.blobs{detinds(ii_d)})
+ self.det_ind = detinds(ii_d);
+ self.get_bboxes_and_datatype();
+ self.get_noise_var_into_grain();
+ end
+ end
+ if nargout
+ gr = self.grain;
+ end
+ end
+
+ function saveClusterToFile(self, phase_id)
+ cl_ids = sort(unique(cat(2, self.cluster.samp_ors.id)), 'ascend');
+ gtSaveCluster(phase_id, cl_ids, self.cluster)
+ end
+
+ function applyDetectorMask(self, detinds)
+ if (~exist('detinds', 'var')) || isempty(detinds)
+ detinds = 1;
+ end
+
+ % fill the missing use_detector_mask for the detectors.
+ nof_missing_flags = max(detinds) - numel(self.use_detector_mask);
+ if (nof_missing_flags > 0)
+ self.use_detector_mask = [self.use_detector_mask(:)', self.use_detector_mask(ones(1, nof_missing_flags))];
+ end
+
+ for ii_d = 1:numel(detinds)
+ if isempty(self.cluster.blobs{detinds(ii_d)})
+ self.use_detector_mask(detinds(ii_d)) = false;
+ else
+ if (self.use_detector_mask(detind))
+ self.load_detector_mask(detinds(ii_d)); % load detector mask from detector%d_mask.edf
+ end
+ if (self.use_detector_mask(detind))
+ self.apply_detector_mask(detinds(ii_d)); % apply the detector mask
+ end
+ end
+ end
+ end
+
+ function inactivatePixelsOutsideDetector(self, detinds)
+ if (~exist('detinds', 'var')) || isempty(detinds)
+ detinds = 1;
+ end
+ for ii_d = 1:numel(detinds)
+ detbbmin = 1 + self.detector_margin;
+ detbbumax = self.parameters.acq(detinds(ii_d)).xdet - self.detector_margin;
+ detbbvmax = self.parameters.acq(detinds(ii_d)).ydet - self.detector_margin;
+ for ii_b = 1:numel(self.cluster.blobs{detinds(ii_d)})
+ bbuim = self.cluster.blobs{detinds(ii_d)}(ii_b).bbuim;
+ bbvim = self.cluster.blobs{detinds(ii_d)}(ii_b).bbvim;
+ bbuim = bbuim(1):bbuim(2);
+ bbvim = bbvim(1):bbvim(2);
+ bbuim = (bbuim > detbbumax) & (bbuim < detbbmin);
+ bbvim = (bbvim > detbbvmax) & (bbvim < detbbmin);
+ if any(bbuim)
+ self.cluster.blobs{detinds(ii_d)}(ii_b).nvar(bbuim, :, :) = self.invalid_pixel_nvar;
+ end
+ if any(bbvim)
+ self.cluster.blobs{detinds(ii_d)}(ii_b).nvar(:, bbvim, :) = self.invalid_pixel_nvar;
+ end
+ end
+ end
+ end
+
+ function deselectBlobCoveredByInvalidPixels(self, detinds)
+ nof_ors = numel(self.cluster.samp_ors);
+ if iscell(self.cluster.blobs)
+ if (~exist('detinds', 'var')) || isempty(detinds)
+ detinds = 1;
+ end
+ for ii_d = 1:numel(detinds)
+ for ii_ors = 1:nof_ors
+ global_pos = self.cluster.samp_ors(ii_ors).proj(detinds(ii_d)).global_pos;
+ for ii_b = 1:numel(global_pos)
+ exist_valid_pixel = self.cluster.blobs{detinds(ii_d)}(global_pos(ii_b)).mask > 0.1;
+ exist_valid_pixel = self.cluster.blobs{detinds(ii_d)}(global_pos(ii_b)).nvar(exist_valid_pixel);
+ exist_valid_pixel = exist_valid_pixel < (self.invalid_pixel_nvar * 0.9);
+ exist_valid_pixel = any(exist_valid_pixel(:));
+ if ~exist_valid_pixel
+ self.cluster.samp_ors(ii_ors).proj(detinds(ii_d)).selected(ii_b) = false;
+ end
+ end
+ end
+ end
+ else
+ for ii_ors = 1:nof_ors
+ global_pos = self.cluster.samp_ors(ii_ors).proj(1).global_pos;
+ for ii_b = 1:numel(global_pos)
+ exist_valid_pixel = self.cluster.blobs(global_pos(ii_b)).mask > 0.1;
+ exist_valid_pixel = self.cluster.blobs(global_pos(ii_b)).nvar(exist_valid_pixel);
+ exist_valid_pixel = exist_valid_pixel < (self.invalid_pixel_nvar * 0.9);
+ exist_valid_pixel = any(exist_valid_pixel(:));
+ if ~exist_valid_pixel
+ self.cluster.samp_ors(ii_ors).proj(1).selected(ii_b) = false;
+ end
+ end
+ end
+ end
+ end
+
+ function showNoiseVariance2DStack(self, detind, lim_max)
+ if (~exist('detind', 'var')) || isempty(detind)
+ detind = 1;
+ end
+ if (~exist('lim_max', 'var')) || isempty(lim_max)
+ lim_max = 1e4;
+ end
+ stack = arrayfun(@(x) sum(x.nvar, 3), self.cluster.blobs{detind}, 'UniformOutput', false);
+ stack = cat(3, stack{:});
+ GtVolView(min(stack, lim_max));
+ end
+ end
+
+ methods (Access = protected)
+ function set_unused_bl_to_empty(self)
+ for ii_ors = 1:numel(self.cluster.samp_ors)
+ for ii_det = 1:numel(self.cluster.samp_ors(ii_ors).proj)
+ for ii_b = 1:numel(self.cluster.samp_ors(ii_ors).proj(ii_det).bl)
+ self.cluster.samp_ors(ii_ors).proj(ii_det).bl(ii_b).intm = [];
+ self.cluster.samp_ors(ii_ors).proj(ii_det).bl(ii_b).mask = [];
+ end
+ end
+ end
+ end
+
+ function normalize_blobs(self)
+ num_det = numel(self.cluster.blobs);
+ for ii_d = 1:num_det
+ nof_bl = numel(self.cluster.blobs{ii_d});
+ for ii_b = 1:nof_bl
+ sum_intm = sum(abs(self.cluster.blobs{ii_d}(ii_b).intm(:)));
+ if ((~isempty(sum_intm)) && (abs(sum_intm) > 1e-5) && (abs(sum_intm - 1) > 1e-5))
+ self.cluster.blobs{ii_d}(ii_b).intm = self.cluster.blobs{ii_d}(ii_b).intm / sum_intm;
+ end
+ end
+ end
+ end
+
+ function get_bboxes_and_datatype(self)
+ detind = self.det_ind;
+ self.data_type = class(self.cluster.blobs{detind}(1).intm);
+ self.bbuims = cat(1, self.cluster.blobs{detind}(:).bbuim);
+ self.bbvims = cat(1, self.cluster.blobs{detind}(:).bbvim);
+
+ % The bbwim is only used to read the median images, so TT
+ % does not need it as there is no median images for TT.
+ if isfield(self.cluster.blobs{detind}, 'bbwim')
+ self.bbwims = cat(1, self.cluster.blobs{detind}(:).bbwim);
+ end
+ end
+
+ function load_detector_mask(self, detind)
+ % load detector mask from detector_mask.edf under
+ % p.acq(detind).dir (currently only work for dct detector 1)
+ detmaskdir = sprintf('detector%d_mask.edf', detind);
+ detmaskdir = fullfile(self.parameters.acq(detind).dir, detmaskdir);
+ if exist(detmaskdir, 'file')
+ det_mask = edf_read(detmaskdir);
+ self.detector_mask{detind} = det_mask' > 0;
+ else
+ self.use_detector_mask(detind) = false;
+ end
+ end
+
+ function get_noise_var_into_grain(self)
+ % get parameters
+ nvar_param = self.get_parameters_of_nvar_generation();
+
+ % get the Poisson noise variance
+ if isempty(nvar_param.nvar_images)
+ nvar_param.nvar_images = 'blob';
+ end
+ fprintf('Loading %s images ...\n', image_type)
+ var_dark = [];
+ dark = [];
+ c = tic;
+ switch nvar_param.nvar_images
+ case 'median' % read background noise from median images. The total noise variance is blob + median image
+ self.copy_bl_int_to_nvar()
+ self.load_bg_from_median_imgs();
+ case 'raw' % directly use raw images as noise variance
+ [var_dark, dark] = GtGrainCalcBlobNoiseVars.calcDarkImageVariance(self.parameters.acq(self.det_ind));
+ self.load_nvar_from_raw_imgs(dark);
+ otherwise % only use blob as noise variance
+ self.copy_bl_int_to_nvar()
+ end
+ fprintf('Done in %.2f s.\n', toc(c));
+
+ % get the noise variance of dark images
+ if (nvar_param.use_dark_image_variance)
+ c = tic;
+ fprintf('Calculating background noise variance from dark images ... ')
+ self.add_noise_of_dark_images(var_dark, dark) % Calculate the noise in dark images
+ fprintf('Done in %.2f s.\n', toc(c));
+ end
+ end
+
+ function nvar_param = get_parameters_of_nvar_generation(self)
+ nvar_param = struct(...
+ 'nvar_images', [], ...
+ 'use_dark_image_variance', false, ...
+ 'save_cluster', false);
+ rec_opt = self.parameters.rec.grains.options;
+ if (isfield(rec_opt, 'detector_weights') && ...
+ ~isempty(rec_opt.detector_weights))
+ weight_opt = rec_opt.detector_weights;
+ if (isfield(weight_opt, 'nvar_images') && ...
+ ischar(weight_opt.nvar_images))
+ nvar_param.nvar_images = lower(weight_opt.nvar_images);
+ end
+ if (isfield(weight_opt, 'use_dark_image_variance') && ...
+ ~isempty(weight_opt.use_dark_image_variance))
+ nvar_param.use_dark_image_variance = weight_opt.use_dark_image_variance;
+ end
+ if (isfield(weight_opt, 'save_cluster') && ...
+ ~isempty(weight_opt.save_cluster))
+ nvar_param.save_cluster = weight_opt.save_cluster;
+ end
+ end
+
+ % As for TT, we can only use the blob as the noise variance.
+ if (~strcmpi(self.parameters.acq(self.det_ind).type, '360degree'))
+ nvar_param.nvar_images = [];
+ nvar_param.use_dark_image_variance = false;
+ end
+ end
+
+ function copy_bl_int_to_nvar(self)
+ switch lower(self.parameters.acq(self.det_ind).type)
+ case {'180degree', '360degree'}
+ blob_type = 'blob';
+ case 'topotomo'
+ blob_type = 'blob_topo';
+ end
+ nof_bl = numel(self.cluster.blobs{self.det_ind});
+ if ~isfield(self.cluster.blobs{self.det_ind}(1), 'nvar')
+ tmpbl = gtFwdSimBlobDefinition(blob_type, nof_bl);
+ for ii_b = nof_bl:-1:1
+ tmpbl(ii_b) = gtAddMatFile(tmpbl(ii_b), self.cluster.blobs{self.det_ind}(ii_b), true, false, false, false, false);
+ end
+ self.cluster.blobs{self.det_ind} = tmpbl;
+ end
+ for ii_b = nof_bl:-1:1
+ self.cluster.blobs{self.det_ind}(ii_b).nvar = abs(self.cluster.blobs{self.det_ind}(ii_b).intm) * self.cluster.blobs{self.det_ind}(ii_b).intensity;
+ end
+ end
+
+ function load_bg_from_median_imgs(self)
+ % According to gtPreprocessing, the median image is regarded as
+ % the background of next parameters.prep.fullint projection
+ % images.
+ acq = self.parameters.acq(self.det_ind);
+ outdir = fullfile(acq.dir, '1_preprocessing', 'full');
+ fullint = self.parameters.prep.fullint;
+ totproj = 2 * acq.nproj;
+ nof_med = ceil(totproj/fullint);
+ gauge = GtGauge(nof_med, 'loading background noise from median images ... ');
+ % loop for median images
+ for fullint_lb = 0:fullint:(fullint * (nof_med - 1))
+ fullint_ub = min(fullint_lb + fullint, totproj);
+ fname = fullfile(outdir, sprintf('med%04d.edf', fullint_lb));
+ gauge.incrementAndDisplay();
+ self.add_median_image_to_spots(fname, fullint_lb, fullint_ub, totproj, acq.xdet, acq.ydet);
+ end
+ end
+
+ function add_median_image_to_spots(self, fname, fullint_lb, fullint_ub, totproj, xdet, ydet)
+ % This function finds the blobs (partially) use this median
+ % image as background and adds this median image to (the part of)
+ % the noise variance of the blobs.
+
+ % truncate the blobs that are partially outside the range of
+ % the median image, and finds the blobs (partially) in the range
+ [bbwims_modified, bool_inint] = self.modify_bbwims_and_check_bl_exist_med(fullint_lb, fullint_ub, totproj);
+
+ if any(bool_inint)
+ % find the blob regions inside the detector.
+ bbuim_inint = self.bbuims(bool_inint, :);
+ bbuim_inint_ondet = GtGrainCalcBlobNoiseVars.calcOndetBbim(bbuim_inint, 1 + self.detector_margin, xdet - self.detector_margin);
+ subbbu = bsxfun(@minus, bbuim_inint_ondet, bbuim_inint(:, 1) - 1);
+ bbvim_inint = self.bbvims(bool_inint, :);
+ bbvim_inint_ondet = GtGrainCalcBlobNoiseVars.calcOndetBbim(bbvim_inint, 1 + self.detector_margin, ydet - self.detector_margin);
+ subbbv = bsxfun(@minus, bbvim_inint_ondet, bbvim_inint(:, 1) - 1);
+ inds_inint = find(bool_inint);
+
+ % The information helping to read the median image
+ info = edf_info(fname);
+
+ % loop for the blobs: read the median image and add it to the blob
+ for ii_inint = 1:numel(inds_inint)
+ ind_b = inds_inint(ii_inint);
+
+ % determine the available bounding box
+ bbwim = bbwims_modified(ind_b, :);
+ bbwlow = max(fullint_lb, bbwim(1)) - bbwim(1) + 1;
+ bbwup = min(fullint_ub - 1, bbwim(2)) - bbwim(1) + 1;
+ bb_ii = [bbuim_inint_ondet(ii_inint, :); bbvim_inint_ondet(ii_inint, :)];
+ bb_ii(:, 2) = bb_ii(:, 2) - bb_ii(:, 1) + 1;
+
+ % read median image in the bounding box
+ medimg_crop = cast(edf_read(fname, bb_ii(:)', [], info), self.data_type)';
+ medimg_crop(medimg_crop <= 0) = self.invalid_pixel_nvar;
+
+ % add the median image to the blob noise variance
+ medbg = medimg_crop(:, :, ones(1, bbwup - bbwlow + 1));
+ subu = subbbu(ii_inint, :);
+ subv = subbbv(ii_inint, :);
+ spotaddbg = self.cluster.blobs{self.det_ind}(ind_b).nvar(subu(1):subu(2), subv(1):subv(2), bbwlow:bbwup) + medbg;
+ self.cluster.blobs{self.det_ind}(ind_b).nvar(:, :, bbwlow:bbwup) = cast(self.invalid_pixel_nvar, self.data_type);
+ self.cluster.blobs{self.det_ind}(ind_b).nvar(subu(1):subu(2), subv(1):subv(2), bbwlow:bbwup) = cast(spotaddbg, self.data_type);
+ end
+ end
+ end
+
+ function [bbwims_modified, bool_inint] = modify_bbwims_and_check_bl_exist_med(self, fullint_lb, fullint_ub, totproj)
+ % truncate the blobs that are partially outside the range of
+ % the median image, and finds the blobs (partially) in the range
+ bbwims_modified = self.bbwims;
+ bool_inint = any(bbwims_modified < (fullint_ub - totproj), 2) & any(bbwims_modified >= (fullint_lb - totproj), 2);
+ if any(bool_inint)
+ bbwims_modified(bool_inint, :) = bbwims_modified(bool_inint, :) - totproj;
+ end
+ bool_inint = any(bbwims_modified < (fullint_ub + totproj), 2) & any(bbwims_modified >= (fullint_lb + totproj), 2);
+ if any(bool_inint)
+ bbwims_modified(bool_inint, :) = bbwims_modified(bool_inint, :) + totproj;
+ end
+ bool_inint = any(bbwims_modified < fullint_ub, 2) & any(bbwims_modified >= fullint_lb, 2);
+ end
+
+ function load_nvar_from_raw_imgs(self, dark)
+ acq = self.parameters.acq(self.det_ind);
+ dir = fullfile(acq.dir, '0_rawdata', acq.name);
+ if (~exist('dark', 'var') || isempty(dark))
+ dark = edf_read(fullfile(dir, 'dark.edf'));
+ end
+ dark = cast(dark', self.data_type); % dark image needs substracting from raw images
+ totproj = 2 * acq.nproj;
+ gauge = GtGauge(totproj, 'loading raw images as noise ... ');
+ self.set_bl_nvar_to_inf(); % initialize the noise variance to a very large value
+ fname = fullfile(acq.dir, '0_rawdata', acq.name, sprintf('%s0000.edf', acq.name));
+ info = edf_info(fname); % The information helping to read the median image
+
+ % find the blob regions inside the detector.
+ bbuims_ondet = GtGrainCalcBlobNoiseVars.calcOndetBbim(self.bbuims, 1 + self.detector_margin, acq.xdet - self.detector_margin);
+ subbbu = bsxfun(@minus, bbuims_ondet, self.bbuims(:, 1) - 1);
+ bbvims_ondet = GtGrainCalcBlobNoiseVars.calcOndetBbim(self.bbvims, 1 + self.detector_margin, acq.ydet - self.detector_margin);
+ subbbv = bsxfun(@minus, bbvims_ondet, self.bbvims(:, 1) - 1);
+
+ % loop for raw images
+ for ii_proj = 0:(totproj-1)
+ gauge.incrementAndDisplay();
+ [bbwims_modified, bool_inproj] = self.modify_bbwims_and_check_bl_exist_raw(ii_proj, totproj);
+ if any(bool_inproj)
+ fname = fullfile(acq.dir, '0_rawdata', acq.name, sprintf('%s%04d.edf', acq.name, ii_proj));
+ inds_inproj = find(bool_inproj);
+
+ % loop for the blobs: read the raw image and add it to the blob
+ for ind_b = inds_inproj(:)'
+ bbwim = bbwims_modified(ind_b, :);
+ bbw = ii_proj - bbwim(1) + 1;
+ bb_ii = [bbuims_ondet(ind_b, :); bbvims_ondet(ind_b, :)];
+ dark_crop = dark(bb_ii(1, 1):bb_ii(1, 2), bb_ii(2, 1):bb_ii(2, 2)); % dark image needs substracting from raw images
+ bb_ii(:, 2) = bb_ii(:, 2) - bb_ii(:, 1) + 1;
+ rawimg_crop = cast(edf_read(fname, bb_ii(:)', [], info), self.data_type)'; % read raw image
+ rawimg_crop(rawimg_crop <= 0) = self.invalid_pixel_nvar;
+ rawimg_crop = rawimg_crop - dark_crop; % dark image needs substracting from raw images
+ subu = subbbu(ind_b, :);
+ subv = subbbv(ind_b, :);
+ self.cluster.blobs{self.det_ind}(ind_b).nvar(subu(1):subu(2), subv(1):subv(2), bbw) = cast(rawimg_crop, self.data_type);
+ end
+ end
+ end
+ end
+
+ function [bbwims_modified, bool_inproj] = modify_bbwims_and_check_bl_exist_raw(self, ii_proj, totproj)
+ bbwims_modified = self.bbwims;
+ bool_inproj = any(bbwims_modified <= (ii_proj - totproj), 2) & any(bbwims_modified >= (ii_proj - totproj), 2);
+ if any(bool_inproj)
+ bbwims_modified(bool_inproj, :) = bbwims_modified(bool_inproj, :) - totproj;
+ end
+ bool_inproj = any(bbwims_modified <= (ii_proj + totproj), 2) & any(bbwims_modified >= (ii_proj + totproj), 2);
+ if any(bool_inproj)
+ bbwims_modified(bool_inproj, :) = bbwims_modified(bool_inproj, :) + totproj;
+ end
+ bool_inproj = any(bbwims_modified <= ii_proj, 2) & any(bbwims_modified >= ii_proj, 2);
+ end
+
+ function set_bl_nvar_to_inf(self)
+ % initialize the noise variance to a very large value
+ switch lower(self.parameters.acq(self.det_ind).type)
+ case {'180degree', '360degree'}
+ blob_type = 'blob';
+ case 'topotomo'
+ blob_type = 'blob_topo';
+ end
+ nof_bl = numel(self.cluster.blobs{self.det_ind});
+ if ~isfield(self.cluster.blobs{self.det_ind}(1), 'nvar')
+ tmpbl = gtFwdSimBlobDefinition(blob_type, nof_bl);
+ for ii_b = nof_bl:-1:1
+ tmpbl(ii_b) = gtAddMatFile(tmpbl(ii_b), self.cluster.blobs{self.det_ind}(ii_b), true, false, false, false, false);
+ end
+ self.cluster.blobs{self.det_ind} = tmpbl;
+ end
+ for ii_b = nof_bl:-1:1
+ self.cluster.blobs{self.det_ind}(ii_b).nvar = cast(ones(size(self.cluster.blobs{self.det_ind}(ii_b).intm)) * self.invalid_pixel_nvar, self.data_type);
+ end
+ end
+
+ function add_noise_of_dark_images(self, var_n_dark, dark)
+ bbuim = self.bbuims;
+ bbvim = self.bbvims;
+
+ % find the range of the detector
+ bbuim_ondet = GtGrainCalcBlobNoiseVars.calcOndetBbim(bbuim, 1 + self.detector_margin, self.parameters.acq(self.det_ind).xdet - self.detector_margin);
+ subbbu = bsxfun(@minus, bbuim_ondet, bbuim(:, 1) - 1);
+ bbvim_ondet = GtGrainCalcBlobNoiseVars.calcOndetBbim(bbvim, 1 + self.detector_margin, self.parameters.acq(self.det_ind).ydet - self.detector_margin);
+ subbbv = bsxfun(@minus, bbvim_ondet, bbvim(:, 1) - 1);
+ bbwsize = abs(self.bbwims(:,2) - self.bbwims(:,1)) + 1;
+
+ detind = self.det_ind;
+ if ~exist('var_n_dark', 'var') || isempty(var_n_dark)
+ [var_n_dark, dark] = GtGrainCalcBlobNoiseVars.calcDarkImageVariance(self.parameters.acq(detind));
+ end
+ var_n_dark(~dark) = self.invalid_pixel_nvar; % set the noise variance of bad pixels to a very large number.
+ var_n_dark = cast(var_n_dark', self.data_type);
+ nof_blobs = numel(self.cluster.blobs{detind});
+ for n = 1:nof_blobs
+ tmp_dark_var = var_n_dark(bbuim_ondet(n, 1):bbuim_ondet(n, 2), bbvim_ondet(n, 1):bbvim_ondet(n, 2));
+ tmp_dark_var = tmp_dark_var(:, :, ones(1, bbwsize(n)));
+ spotaddbg = self.cluster.blobs{detind}(n).nvar(subbbu(n, 1):subbbu(n, 2), subbbv(n, 1):subbbv(n, 2), :) + tmp_dark_var;
+
+ % set the noise variance out of the detector to a very large number
+ self.cluster.blobs{detind}(n).nvar = cast(ones(size(self.cluster.blobs{detind}(n).nvar)) * self.invalid_pixel_nvar, self.data_type);
+
+ self.cluster.blobs{detind}(n).nvar(subbbu(n, 1):subbbu(n, 2), subbbv(n, 1):subbbv(n, 2), :) = cast(spotaddbg, self.data_type);
+ end
+ end
+
+ function apply_detector_mask(self, detind)
+ nof_bl = numel(self.cluster.blobs{detind});
+ size_detmask = size(self.detector_mask);
+ for ii_b = 1:nof_bl
+ size_nvar = size(self.cluster.blobs{detind}(ii_b).nvar);
+ nvarinactmask = false(size_nvar(1), size_nvar(2));
+ bbim = self.cluster.blobs{detind}(ii_b).bbuim;
+ bbuinds = bbim(1):bbim(2);
+ bool_bbuindsondet = ~((bbuinds < 1) | (bbuinds > size_detmask(1)));
+ bbim = self.cluster.blobs{detind}(ii_b).bbvim;
+ bbvinds = bbim(1):bbim(2);
+ bool_bbvindsondet = ~((bbvinds < 1) | (bbvinds > size_detmask(2)));
+ nvarinactmask(bool_bbuindsondet, bool_bbvindsondet) = ~self.detector_mask(bbuinds(bool_bbuindsondet), bbvinds(bool_bbvindsondet));
+ nvarinactmask = nvarinactmask(:, :, ones(1, size_nvar(3)));
+ self.cluster.blobs{detind}(ii_b).nvar(nvarinactmask) = self.invalid_pixel_nvar;
+ end
+ end
+ end
+end
\ No newline at end of file
--
GitLab