gt6DCreateProjDataFromTwinnedGrainFwdProj.m

function [cl, estim_space_bbox_pix, estim_orient_bbox_rod] = gt6DCreateProjDataFromTwinnedGrainFwdProj(parent_id, variants, phase_id, varargin)
% FUNCTION cl = gt6DCreateProjDataFromTwinnedGrainFwdProj(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'))
        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, ...
        'import_ref_included', true, ...
        'use_raw_images', false, ...
        'use_parent_mask', false, ...
        'include_all', false, ...
        'twin_only_shared_blobs', 'deselect', ... % {'deselect'} | 'ignore' | 'use'
        'save', false );
    conf = parse_pv_pairs(conf, varargin);

    if (~any(conf.use_raw_images) && conf.include_all)
        warning('gt6DCreateProjDataFromTwinnedGrainFwdProj:wrong_argument', ...
            '"include_all" cannot be enabled if "use_raw_images" is not enabled')
        conf.include_all = false;
    end

    p = gtLoadParameters();
    symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);

    if (numel(parent_id) > 1)
        error('gt6DCreateProjDataFromTwinnedGrainFwdProj:wrong_argument', ...
            'Only one grain should be passed as parent')
    end
    if (~isstruct(parent_id))
        grs(1) = gtLoadGrain(phase_id, parent_id);
        grs(1) = gtCalculateGrain(grs(1), p);
    else
        grs = parent_id;
    end

    % For compatibility
    if (isfield(grs(1), 'g_twins') && ~isempty(grs(1).g_twins))
        grs(1).twin_vars = grs(1).g_twins;
    end

    if (isfield(grs(1), 'twin_vars') && ~isempty(grs(1).twin_vars))
        theo_fwdproj_variants = cat(1, grs(1).twin_vars(:).R_vector);
    else
        theo_fwdproj_variants = [];
    end

    variants_num = numel(variants);
    fprintf('Loading variants: ')
    for ii_var = variants_num:-1:1
        num_chars = fprintf('%02d (type: %s)', variants_num-ii_var+1, variants(ii_var).type);

        switch (variants(ii_var).type)
            case {'gr_id', 'gr_str'}
                try
                    if (strcmpi(variants(ii_var).type, 'gr_id'))
                        gr_var_id = variants(ii_var).data;
                        grs(ii_var+1) = gtLoadGrain(phase_id, variants(ii_var).data);
                    else
                        gr_var_id = variants(ii_var).data.id;
                        grs(ii_var+1) = variants(ii_var).data;
                    end
                    grs(ii_var+1) = gtCalculateGrain(grs(ii_var+1), p);
                catch mexc
                    gtPrintException(mexc, 'Reforward-projecting grain because of:')
                    grs(ii_var+1) = gtForwardSimulate_v2(gr_var_id, gr_var_id, ...
                        pwd, phase_id, 'save_grain', false, 'display_figure', false);
                end

                if (~isempty(theo_fwdproj_variants))
                    dists = theo_fwdproj_variants - repmat(grs(ii_var+1).R_vector, [size(theo_fwdproj_variants, 1), 1]);
                    dists = 2 * atand(sqrt(sum(dists .^ 2, 2)));
                    indx = find(dists < 1);
                    if (~isempty(indx))
                        fprintf('\n - Corresponding to precomputed variant: %d\n', indx);
                        fprintf(['Loading variants: ' repmat(' ', [1 num_chars])]);
                    end
                end
            case 'var_num'
                grs(ii_var+1) = twin_fwd_sim_to_gr(grs(1), grs(1).twin_vars(variants(ii_var).data));
            otherwise
                error('gt6DCreateProjDataFromTwinnedGrainTheoFwdProj:wrong_argument', ...
                    'Variant type: %s cannot be handled', variants(ii_var).type)
        end

        fprintf(repmat('\b', [1 num_chars]));
        fprintf(repmat(' ', [1 num_chars]));
        fprintf(repmat('\b', [1 num_chars]));
    end
    fprintf('Done.\n')

    grs_list = cat(2, grs(:).id);
    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);
    end

    refgr_proj = refgr.proj(conf.det_index);

    if (numel(conf.use_raw_images) == 1)
        conf.use_raw_images = conf.use_raw_images(ones(num_grains, 1));
    end

    ref_ondet = refgr_proj.ondet;
    if (conf.use_raw_images(1) && conf.include_all)
        ref_included = 1:numel(ref_ondet);
        ref_selected = false(size(ref_included));
        ref_selected(refgr_proj.included(refgr_proj.selected)) = true;
    else
        ref_included = refgr_proj.included;
        ref_selected = refgr_proj.selected;
    end

    fprintf('Determining real-space bounding box..')
    c = tic();

    % Determining real-space bounding box. It is assumed to be contiguos.
    [gr_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(grs, conf.det_index, conf.rspace_oversize);
    gr_center_mm = gtGeoSam2Sam(gr_center_pix, p.recgeo, p.samgeo, false, false);

    estim_space_bbox_mm = [ ...
        gtGeoSam2Sam(estim_space_bbox_pix(1:3), p.recgeo, p.samgeo, false, false), ...
        gtGeoSam2Sam(estim_space_bbox_pix(4:6), p.recgeo, p.samgeo, false, false) ];
    bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3);

    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);

    max_img_sizes = zeros(num_grains, 2);

    extended_projs = get6DExtendedProjDefinition(num_grains);
    inconvenient_etas = cell(num_grains, 1);

    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 orienation 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

    samp_ors = cell(num_grains, 1);
    size_refl_list = size(grs(1).allblobs(conf.det_index).omega);

    % Determining the regions to be merged!
    merge_lists = determine_merge_list(estim_orient_bbox_rod, num_grains, conf);

    if (numel(merge_lists) < num_grains)
        num_grains = numel(merge_lists);
        % Merging regions
        new_grs = cell(num_grains, 1);
        for ii_m = 1:num_grains
            new_gr = grs(merge_lists{ii_m}(1));
            for ii_g = merge_lists{ii_m}(2:end)
                new_gr = gtGrainMerge(new_gr, grs(ii_g), p);
            end
            new_grs{ii_m} = new_gr;
        end
        grs = [new_grs{:}];
        % 'any' is ust a trick to deal with 1 and multiple entries
        conf.use_raw_images = cellfun(@(x)(any(conf.use_raw_images(x)) || numel(x) > 1), merge_lists);

        estim_orient_bbox_rod_old = estim_orient_bbox_rod;
        estim_orient_bbox_rod = zeros(num_grains, 6);
        for ii_g = 1:num_grains
            estim_orient_bbox_rod(ii_g, :) = [ ...
                min(estim_orient_bbox_rod_old(merge_lists{ii_g}, 1:3), [], 1), ...
                max(estim_orient_bbox_rod_old(merge_lists{ii_g}, 4:6), [], 1) ...
                ];
        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
            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
        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-1, 1);
    matches_twin_to_ref = cell(num_grains-1, 1);

    for ii_g = 1:num_grains-1
        % At the moment we only support the configuration: parent + twins
        % 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-ii_g);
        for ii_g_2 = (ii_g+1):num_grains
            twin_info = gtCrystTwinTest(samp_ors(ii_g).R_vector, samp_ors(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(ii_g), samp_ors((ii_g+1):end), angular_tolerance, conf.det_index);
        fprintf('\b\b: Done.\n')

        if (conf.verbose > 1)
            produce_matching_reflections_table(samp_ors(ii_g:end), conf, angular_tolerance, matches_ref_to_twin{ii_g})
        end
    end

    for ii_g = 1:num_grains
        % Populating the extended projection data-structures
        if (ii_g == 1)
            % We first deal with the parent grain, and we will create the
            % datastructures: refl_matches_ref and refl_matches_twin
            extended_projs(ii_g).ondet = ref_ondet;
            extended_projs(ii_g).included = ref_included;
            extended_projs(ii_g).selected = ref_selected;

            bool_ref_ondet = false(size_refl_list);
            bool_ref_ondet(ref_ondet) = true;
            bool_ref_inc = false(size_refl_list);
            bool_ref_inc(ref_ondet(ref_included)) = true;
            bool_ref_sel = false(size_refl_list);
            bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;

            extended_projs(ii_g).allreflbool_ondet = bool_ref_ondet;
            extended_projs(ii_g).allreflbool_included = bool_ref_inc;
            extended_projs(ii_g).allreflbool_selected = bool_ref_sel;
        else
            extended_projs(ii_g) = find_matches_with_refor(extended_projs(1), ...
                grs(ii_g), matches_ref_to_twin{1}(:, ii_g-1), matches_twin_to_ref{1}(:, ii_g-1), conf);

            extended_projs(1).selected = extended_projs(1).selected ...
                | extended_projs(ii_g).shared_to_be_selected_in_parent;

            extended_projs(1).allreflbool_selected(ref_ondet(ref_included)) = extended_projs(1).selected;
        end

        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

    % Following code fixes the problem of two twin variants might share the
    % same spot with the parent but the spot might not have been selected
    % in one of them - so this would lead to inconsistent projections
    for ii_g = 2 : num_grains
        extended_projs(ii_g) = find_matches_with_refor(extended_projs(1), ...
            grs(ii_g), matches_ref_to_twin{1}(:, ii_g-1), matches_twin_to_ref{1}(:, ii_g-1), conf);

        % (2) twin variants might share reflections which are not
        % shared with the parent - we currently cannot handle these, so
        % we deselect them
        switch (conf.twin_only_shared_blobs)
            case 'deselect'
                [twin_only_shared_bls, t2t_shr_bls_tbl] = detect_twin_only_shared_reflections(matches_ref_to_twin, matches_twin_to_ref, ii_g);
                twin_only_shared_bls_inc = twin_only_shared_bls(extended_projs(ii_g).ondet(extended_projs(ii_g).included));

                deselected = extended_projs(ii_g).selected & twin_only_shared_bls_inc;
                extended_projs(ii_g).selected = extended_projs(ii_g).selected & ~twin_only_shared_bls_inc;

                extended_projs(ii_g).allreflbool_selected(extended_projs(ii_g).ondet(extended_projs(ii_g).included)) = extended_projs(ii_g).selected;

                if (conf.verbose > 1 || true)
                    produce_deselected_reflections_table(samp_ors, extended_projs, ii_g, deselected, twin_only_shared_bls_inc, t2t_shr_bls_tbl, conf);
                end
            case 'ignore'
                % we do nothing
            otherwise
                error('gt6DCreateProjDataFromTwinnedGrainFwdProj:wrong_argument', ...
                    'Option: %s for twin_only_shared_blobs not implemented yet!', ...
                    conf.twin_only_shared_blobs)
        end

        if (any(extended_projs(ii_g).shared_refl_inc) && (conf.verbose || true))
            parent_ids = grs(1).fwd_sim(conf.det_index).spotid(ref_included(extended_projs(ii_g).shared_refl_pos_in_parent_inc))';

            twin_ids = zeros(1, numel(extended_projs(ii_g).shared_refl_inc));
            inds_of_replaced_blobs = grs(ii_g).proj(conf.det_index).included(~extended_projs(ii_g).included_old_to_keep);
            twin_ids(extended_projs(ii_g).included_replaced) = grs(ii_g).fwd_sim(conf.det_index).spotid(inds_of_replaced_blobs)';
            twin_ids = twin_ids(extended_projs(ii_g).shared_refl_inc);

            fprintf('Twin (%2d) blobs, shared blobs with reference.\n', ii_g)
            fprintf('Replaced by parents blobs (if 0, they were not included before):\n')
            fprintf(' - %02d: twin blob_id %5d -> parent blob_id %5d\n', ...
                [1:numel(parent_ids); twin_ids; parent_ids])
        end
    end

    fprintf('Determining UVW bounding boxes..')
    c = tic();

    if (any(conf.use_raw_images))
        uvw_tab = cell(num_grains, 1);
        img_bboxes = cell(num_grains, 1);
        img_sizes = cell(num_grains, 1);
    end

    for ii_g = 1:num_grains
        if (conf.use_raw_images(ii_g))
            local_ondet = extended_projs(ii_g).ondet;
            local_included = extended_projs(ii_g).included;
            local_inc_refl = local_ondet(local_included);

            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;

            max_img_sizes(ii_g, :) = max(img_sizes{ii_g}(:, 1:2), [], 1);

            if (conf.verbose || true)
                or_ns = samp_ors(ii_g).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{ii_g}, or_ns, or_ws * gtAcqGetOmegaStep(p, conf.det_index), ...
                    ~inconvenient_etas{ii_g}, ...
                    img_bboxes{ii_g}(:, [1 4 2 5 3 6]), ...
                    extended_projs(ii_g).selected ];

                fprintf('\nImages for orientation %d:\n', ii_g)
                if (ii_g == 1)
                    fprintf([tbl_out_format '\n'], ...
                        tbl_output');
                else
                    pos_of_shared_blobs = zeros(size(extended_projs(ii_g).shared_refl_inc));
                    pos_of_shared_blobs(extended_projs(ii_g).shared_refl_inc) = extended_projs(ii_g).shared_refl_pos_in_parent_inc;

                    fprintf([tbl_out_format ', shared: %d\n'], ...
                        [tbl_output,  pos_of_shared_blobs]');
                end
            end
        else
            max_img_sizes(ii_g, :) = [ ...
                size(grs(ii_g).proj(conf.det_index).stack, 1), ...
                size(grs(ii_g).proj(conf.det_index).stack, 3), ];
        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_g = 1:num_grains
        fprintf('Processing orientation-space region: %d/%d\n', ii_g, num_grains)

        num_blobs = numel(extended_projs(ii_g).included);

        c = tic();
        if (conf.use_raw_images(ii_g))
            fprintf(' - Loading raw images: ')

            blobs = gtFwdSimBlobDefinition('blob', num_blobs);
            for ii_b = 1:num_blobs
                num_chars = fprintf('%03d/%03d', ii_b, num_blobs);

                if (ii_g == 1 || ~extended_projs(ii_g).shared_refl_inc(ii_b))
                    blobs(ii_b) = load_blob( ...
                        img_bboxes{ii_g}(ii_b, :), img_sizes{ii_g}(ii_b, :), ...
                        stackUSize, stackVSize, p, conf);
                end

                fprintf(repmat('\b', [1, num_chars]));
            end
            % Using the same because it gives the same size/uvw-limits
            if (ii_g > 1)
                blobs(extended_projs(ii_g).shared_refl_inc) ...
                    = grs(1).proj(conf.det_index).bl(extended_projs(ii_g).shared_refl_pos_in_parent_inc);
            end
        else
            fprintf(' - Loading segmented blobs..')

            if (ii_g == 1)
                blobs = grs(ii_g).proj.bl;
            else
                blobs = gtFwdSimBlobDefinition('blob', num_blobs);
                blobs(~extended_projs(ii_g).shared_refl_inc) ...
                    = grs(ii_g).proj(conf.det_index).bl(extended_projs(ii_g).included_old_to_keep);
                % Using the same because it gives the same size/uvw-limits
                blobs(extended_projs(ii_g).shared_refl_inc) ...
                    = grs(1).proj(conf.det_index).bl(extended_projs(ii_g).shared_refl_pos_in_parent_inc);
            end
        end

        if (conf.use_raw_images(ii_g) && conf.use_parent_mask && ii_g > 1)
            fprintf('\b\b: Done in %g seconds.\n - Producing blob masks..', toc(c))
            c = tic();

            proj_bl_masks = project_masks(samp_ors(ii_g), extended_projs(ii_g), grs(1).fwd_sim(conf.det_index).gv_verts, p, conf);
            blobs = assign_masks(blobs, proj_bl_masks, ~extended_projs(ii_g).shared_refl_inc);
        end

        fprintf('\b\b: Done in %g seconds.\n - Equalizing blob sizes: ', toc(c))
        c = tic();

        for ii_b = 1:num_blobs
            num_chars = fprintf('%02d/%02d', ii_b, num_blobs);

            blobs(ii_b) = equalize_blob_size(blobs(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 - Producing proj data-structure..', toc(c))
        c = tic();

        spots = arrayfun(@(x){sum(x.intm, 3)}, blobs);
        spots = permute(cat(3, spots{:}), [1 3 2]);

        proj = gtFwdSimProjDefinition();

        proj.centerpix = gr_center_pix;
        proj.bl = blobs;
        proj.stack = spots;
        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_g).ondet;
        proj.included = extended_projs(ii_g).included;
        proj.selected = extended_projs(ii_g).selected;

        if (ii_g > 1)
            shared_inc = extended_projs(ii_g).shared_refl_inc;
            non_shared_blobs = blobs(~shared_inc);
            proj.global_pos = zeros(size(proj.included));
            proj.global_pos(~shared_inc) = numel(all_blobs_list) + (1:numel(non_shared_blobs))';
            proj.global_pos(shared_inc) = extended_projs(ii_g).shared_refl_pos_in_parent_inc;
            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(ii_g).proj(conf.det_index) = proj;

        fprintf('\b\b: Done in %g seconds.\n', toc(c))

        if (conf.check_spots && ii_g > 1)
            samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true, conf.det_index);
        end
    end

    if (conf.verbose)
        f = figure();
        ax = axes('parent', f);
        hold(ax, 'on');
        for ii_g = 1:num_grains
            gt6DPlotOrientationBBox(ax, cat(1, samp_ors(ii_g).bb_ors(:).R_vector));
            scatter3(ax, 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);

    if (conf.save)
        str_ids = sprintf('_%04d', grs_list);
        grain_filename = fullfile(p.acq.dir, '4_grains', ...
            sprintf('phase_%02d', phase_id), ...
            sprintf('grain_cluster%s.mat', str_ids));
        fprintf('Saving the cluster file "%s"..', grain_filename)
        save(grain_filename, '-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, 1);
        sample.phases{phase_id}.setClusterInfo(grs_list, cl_info);
        sample.saveToFile();
        fprintf('\b\b: Done.\n')
    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 gr = twin_fwd_sim_to_gr(ref_gr, twin)
    gr = ref_gr;
    gr.R_vector = twin.R_vector;
    gr.allblobs = twin.allblobs;
    gr.proj = twin.proj;
    gr.fwd_sim = twin.fwd_sim;

    % And the useless stuff
    gr.check     = twin.fwd_sim(1).check;
    gr.flag      = twin.fwd_sim(1).flag;
    gr.spotid    = twin.fwd_sim(1).spotid;

    % only included spots info
    gr.cands     = twin.fwd_sim(1).cands;
    gr.likelihood = twin.fwd_sim(1).likelihood;
    gr.difspotID = twin.fwd_sim(1).difspotID;
    gr.uv_exp    = twin.fwd_sim(1).uvw(:, 1:2);
    gr.om_exp    = twin.fwd_sim(1).uvw(:, 3);
    gr.bb_exp    = twin.fwd_sim(1).bb;
    gr.intensity = twin.fwd_sim(1).intensity;
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 produce_deselected_reflections_table(samp_ors, eproj, ii_g, deselected, twin_only_shared_bls_inc, t2t_shr_bls_tbl, conf)

    vis_inds = eproj(ii_g).ondet(eproj(ii_g).included(twin_only_shared_bls_inc));
    if (size(samp_ors(ii_g).allblobs(conf.det_index).hklsp, 2) == 4)
        hklsp = samp_ors(ii_g).allblobs(conf.det_index).hklsp(:, [1 2 4]);
    else
        hklsp = samp_ors(ii_g).allblobs(conf.det_index).hklsp;
    end

    [refl_inds, twin_ind] = find(t2t_shr_bls_tbl(vis_inds, :));

    fprintf('Twin (%2d) blobs identified due to overlap with only other twins:\n', ii_g)
    if (numel(vis_inds))
        fprintf(' - %02d: w %6.2f, n %6.2f, w_ind %d, hkl [%2d %2d %2d], deselected: %d\n', ...
            [(1:numel(vis_inds))', samp_ors(ii_g).allblobs(conf.det_index).omega(vis_inds), ...
            samp_ors(ii_g).allblobs(conf.det_index).eta(vis_inds), samp_ors(ii_g).allblobs(conf.det_index).omind(vis_inds), ...
            hklsp(vis_inds, :), deselected(twin_only_shared_bls_inc)]');
    else
        fprintf(' - None\n')
    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, ...
        'shared_refl_inc', repl_cell, ...
        'shared_refl_pos_in_parent_inc', repl_cell, ...
        'shared_to_be_selected_in_parent', repl_cell);
end

function eproj = find_matches_with_refor(ref_eproj, twingr, matches_ref_to_twin, matches_twin_to_ref, conf)
    % first translates from numel(gr.proj.ondet) into numel(gr.allblobs(conf.det_index).omega)
    % "indexing / numbering" frame, then identifies shared reflections and
    % translates back into numel(gr.proj.ondet) frame...

    eproj = get6DExtendedProjDefinition();

    size_allbls_refl_list = size(ref_eproj.allreflbool_ondet);

    ref_ondet = ref_eproj.ondet;
    ref_included = ref_eproj.included;

    % bool_* vectors are all of length numel(gr.allblobs(conf.det_index).omega)
    bool_ref_inc = ref_eproj.allreflbool_included;
    bool_ref_sel = ref_eproj.allreflbool_selected;

    %  the row index of refl_matches_ref corresponds to the hklsp index of
    %  shared reflection in the parent - the "value"points to the hklsp
    %  row index of the twin:  parent_id -> twin_id
    % list of parent reflections for which a shared spot has been
    % predicted with the current twin variant

    % bool_ref_shared  = false(size_refl_list);
    % bool_ref_shared(find(refl_matches_ref > 0)) = true;
    % list of twin reflections for which a shared spot has been
    % predicted
    bool_twin_shared = false(size_allbls_refl_list);
    bool_twin_shared(matches_twin_to_ref > 0) = true;

    twingr_proj = twingr.proj(conf.det_index);

    twin_ondet = twingr_proj.ondet;
    twin_inc = twingr_proj.included;
    twin_sel = twingr_proj.selected;

    bool_twin_ondet = false(size_allbls_refl_list);
    bool_twin_ondet(twin_ondet) = true;
    bool_twin_inc = false(size_allbls_refl_list);
    bool_twin_inc(twin_ondet(twin_inc)) = true;
    bool_twin_sel = false(size_allbls_refl_list);
    bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;

    % parent indices of shared reflections that are ondet for the twin
    ref_shared_ondet = false(size_allbls_refl_list);
    ref_shared_ondet(matches_twin_to_ref(bool_twin_shared & bool_twin_ondet)) = true;
    % ...and which were included for the parent
    ref_shared_ondet_inc = bool_ref_inc & ref_shared_ondet;
    ref_shared_selected = ref_shared_ondet_inc & bool_ref_sel;

    % translate this back into twin indices
    twin_shared_ref_inc = false(size_allbls_refl_list);
    twin_shared_ref_inc(matches_ref_to_twin(ref_shared_ondet_inc)) = true;

    % These will tell which ones to keep from the fwd-simulation
    bool_twin_inc_keep = bool_twin_inc & ~twin_shared_ref_inc;

    % this is used for tracking which blobs were previously included and
    % then replaced in favor of reference's blobs
    bool_twin_inc_replaced = bool_twin_inc & twin_shared_ref_inc;

    eproj.included_old_to_keep = bool_twin_inc_keep(twin_ondet(twin_inc));

    % These wll be the ones that will be included from the parent
    if (conf.import_ref_included)
        bool_twin_inc_import = bool_twin_ondet & twin_shared_ref_inc;

        bool_twin_inc(bool_twin_inc_import) = true;
        twin_inc = find(bool_twin_inc(twin_ondet));
    else
        bool_twin_inc_import = bool_twin_inc & twin_shared_ref_inc;
    end
    % These are the ones that were included and replaced. In case the
    % import option is false, these are identical to shared
    eproj.included_replaced = bool_twin_inc_replaced(twin_ondet(twin_inc));

    % If a shared & twin included spot is selected in parent, select
    % it for the twin, as well
    bool_twin_sel(matches_ref_to_twin(ref_shared_selected)) = bool_twin_inc(matches_ref_to_twin(ref_shared_selected));
    twin_sel = bool_twin_sel(twin_ondet(twin_inc));

    bool_ref_matched_twin_inc = false(size_allbls_refl_list);
    bool_ref_matched_twin_inc(matches_twin_to_ref(bool_twin_inc_import)) = true;
    shared_to_be_selected_in_parent = bool_ref_matched_twin_inc & ~bool_ref_sel;

    eproj.shared_to_be_selected_in_parent ...
        = shared_to_be_selected_in_parent(ref_ondet(ref_included));

    % We now find the blobs in the old twin bl structure that are
    % shared with the parent
    eproj.shared_refl_inc = twin_shared_ref_inc(twin_ondet(twin_inc));

    % We have to find the included IDs of the corresponding shared
    % reflections
    refl_pos_in_parent_inc = zeros(size_allbls_refl_list);
    refl_pos_in_parent_inc(ref_ondet(ref_included)) = (1:numel(ref_included));

    eproj.shared_refl_pos_in_parent_inc = refl_pos_in_parent_inc(matches_twin_to_ref(bool_twin_inc_import));

    eproj.ondet = twin_ondet;
    eproj.included = twin_inc;
    eproj.selected = twin_sel;
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();

    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 = true(blob_size_im);

    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 proj_bl_masks = project_masks(refor, eproj, parent_verts, p, conf)
    sel = eproj.ondet(eproj.included);
    proj_bl_masks = gt6DSpreadProjectVertices2Det(refor, parent_verts, sel, p, conf.det_index);
end

function blobs = assign_masks(blobs, proj_mask_blobs, sel)
    blobs_u_lims = cat(1, blobs(sel).bbuim);
    blobs_v_lims = cat(1, blobs(sel).bbvim);

    proj_masks_u_lims = cat(1, proj_mask_blobs(sel).bbuim);
    proj_masks_v_lims = cat(1, proj_mask_blobs(sel).bbvim);

    shifts = [ ...
        (proj_masks_u_lims(:, 1) - blobs_u_lims(:, 1)), ...
        (proj_masks_v_lims(:, 1) - blobs_v_lims(:, 1)) ];
    shifts(:, 3) = 0;

    inds = find(sel);
    for ii = 1:numel(inds)
        ii_b = inds(ii);

        mask = uint8(proj_mask_blobs(ii_b).mask);
        mask = mask(:, :, ones(blobs(ii_b).bbsize(3), 1));

        blobs(ii_b).mask = gtPlaceSubVolume( ...
            zeros(blobs(ii_b).bbsize, 'uint8'), mask, shifts(ii, :) );

        blobs(ii_b).intm(~blobs(ii_b).mask) = 0;
        blob_int = gtMathsSumNDVol(blobs(ii_b).intm);
        blobs(ii_b).intensity = blob_int;
    end
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 = 2: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

    to_clean_up = false(num_grains, 1);
    for ii_m = 1:num_grains
        if (numel(merge_lists{ii_m}) == 1 && merge_lists{ii_m} == 0)
            to_clean_up(ii_m) = true;
        else
            merge_lists{ii_m} = explore_list(merge_lists, merge_lists{ii_m});
            for ii_g = merge_lists{ii_m}
                % Removing reference to other, so to avoid
                % double merging
                merge_lists{ii_g} = 0;
            end
            merge_lists{ii_m} = [ii_m, merge_lists{ii_m}];
        end
    end
    merge_lists(to_clean_up) = [];

    function partial_list_out = explore_list(merge_lists, partial_list_in)
        partial_list_out = partial_list_in;
        for ii = 1:numel(partial_list_in)
            entries_list = merge_lists{partial_list_in(ii)};
            partial_list_out = [partial_list_out, explore_list(merge_lists, entries_list)]; %#ok<AGROW>
        end
    end
end