gt6DCreateProjDataFromTwinnedGrain.m

function [cl, estim_space_bbox_pix, estim_orient_bbox_rod] = gt6DCreateProjDataFromTwinnedGrain(grains_ids, 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( 'parent_ind', <ind>, 'var_num', <ind> )
% 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, ...
        'ospace_lims', [], ...
        'check_spots', false, ...
        'stack_oversize', 1.4, ...
        'merge_dist_pixels', 0, ...
        'match_tol', [0, 0, 1], ...
        'use_raw_images', false, ...
        'use_volume_mask', false, ...
        'prefer_raw_shared', false, ...
        'include_all', false, ...
        'always_select_shared', false, ...
        '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 (~isstruct(grains_ids))
        grs = gtLoadGrain(phase_id, grains_ids);
        for ii_g = 1:numel(grs)
            grs(ii_g) = gtCalculateGrain(grs(ii_g), p);
        end
    else
        grs = grains_ids;
    end

    variants_num = numel(variants);
    if (variants_num > 0)
        fprintf('Loading variants: ')
        for ii_var = variants_num:-1:1
            n_p = variants(ii_var).parent_ind;
            n_v = variants(ii_var).var_num;

            num_chars = fprintf('%02d (parent indx: %d, var_num: %d)', variants_num-ii_var+1, n_p, n_v);

            grs(end+1) = twin_fwd_sim_to_gr(grs(n_p), grs(n_p).twin_vars(n_v)); %#ok<AGROW>

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

    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

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

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

    % Determining real-space bounding box. It is assumed to be contiguos.
    [cl_center_pix, bbox_size_pix, estim_space_bbox_pix] = gt6DMergeRealSpaceVolumes(grs, conf.det_index, conf.rspace_oversize);
    cl_center_mm = gtGeoSam2Sam(cl_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);

    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)

        if (size(conf.ospace_lims, 1) >= ii_g)
            diff_r_vecs = tand(conf.ospace_lims(ii_g, :) / 2);
            r_vecs = [grs(ii_g).R_vector + diff_r_vecs(1:3); grs(ii_g).R_vector + diff_r_vecs(4:6)];
        else
            sampler = GtOrientationSampling(p, grs(ii_g), ...
                'detector_index', conf.det_index, 'verbose', conf.verbose);
            r_vecs = sampler.guess_ODF_BB()';
        end

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

    % Computing the new grain structures
    for ii_g = 1:num_grains
        refor = struct( ...
            'id', grs(ii_g).id, 'phaseid', grs(ii_g).phaseid, ...
            'center', cl_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{:}];

    num_images = gtAcqTotNumberOfImages(p);

    % Let's now list all blobs, and the corresponding UVW bounding boxes
    % from the orientation-bbs projections (not segmentation of the bls)
    all_blobs_list = grs(1).proj(conf.det_index).bl;
    for ii_g = 2:num_grains
        all_blobs_list = cat(2, all_blobs_list, grs(ii_g).proj(conf.det_index).bl);
    end
    num_blobs = numel(all_blobs_list);
    for ii_b = 1:num_blobs
        all_blobs_list(ii_b).intm = all_blobs_list(ii_b).intm ...
            / gtMathsSumNDVol(all_blobs_list(ii_b).intm) ...
            * all_blobs_list(ii_b).intensity;
    end

    blobs_uvw_boundaries = zeros(num_blobs, 6);
    blobs_ids = zeros(num_blobs, 1);
    base_counter = 0;
    is_selected = false(num_blobs, 1);
    blob_use_raw_images = false(num_blobs, 1);
    blobs_etas = zeros(num_blobs, 1);

    for ii_g = 1:num_grains
        num_bl_grain = numel(grs(ii_g).proj(conf.det_index).bl);
        inds_bls = (base_counter+1):(base_counter+num_bl_grain);

        inc_bls = grs(ii_g).proj(conf.det_index).ondet(grs(ii_g).proj(conf.det_index).included);

        num_bb_ors = numel(samp_ors(ii_g).bb_ors);
        if (conf.use_raw_images(ii_g))
            uvw_tab = zeros(numel(inc_bls), 3, num_bb_ors);
            for ii_ab = 1:num_bb_ors
                uvw_tab(:, :, ii_ab) = samp_ors(ii_g).bb_ors(ii_ab).allblobs(conf.det_index).detector.uvw(inc_bls, :);
            end
        else
            bls_u_lims = cat(1, grs(ii_g).proj(conf.det_index).bl(:).mbbu);
            bls_v_lims = cat(1, grs(ii_g).proj(conf.det_index).bl(:).mbbv);
            bls_w_lims = cat(1, grs(ii_g).proj(conf.det_index).bl(:).mbbw);
            uvw_tab = cat(3, ...
                [bls_u_lims(:, 1), bls_v_lims(:, 1), bls_w_lims(:, 1)], ...
                [bls_u_lims(:, 2), bls_v_lims(:, 2), bls_w_lims(:, 2)] );
        end
        blobs_uvw_boundaries(inds_bls, 1:3) = min(floor(uvw_tab), [], 3);
        blobs_uvw_boundaries(inds_bls, 4:6) = max(ceil(uvw_tab), [], 3);
        blobs_uvw_boundaries(inds_bls, [3 6]) = gtGrainAnglesTabularFix360deg( ...
            blobs_uvw_boundaries(inds_bls, [3 6]), samp_ors(ii_g).allblobs(conf.det_index).detector.uvw(inc_bls, 3), p);

        blobs_etas(inds_bls) = samp_ors(ii_g).allblobs(conf.det_index).eta(inc_bls);
        blobs_ids(inds_bls) = grs(ii_g).fwd_sim(conf.det_index).difspotID;
        is_selected(inds_bls) = grs(ii_g).proj(conf.det_index).selected;
        blob_use_raw_images(inds_bls) = conf.use_raw_images(ii_g);

        base_counter = base_counter + num_bl_grain;
    end

    % Transforming blob boundaries to blob center + halfwidth
    blobs_uvw_center_width = [ ...
        blobs_uvw_boundaries(:, 4:6) + blobs_uvw_boundaries(:, 1:3), ...
        blobs_uvw_boundaries(:, 4:6) - blobs_uvw_boundaries(:, 1:3) ] / 2;

    is_merged = false(size(blobs_ids));
    % We will now scan for blobs which have the same ID and merge them
    for ii_b = 1:num_blobs
        ind_b_others = (ii_b+1:num_blobs)';

        matches = (blobs_ids(ii_b) == blobs_ids(ind_b_others)) & ~is_merged(ind_b_others);

        if (any(matches))
            merge_list = [ii_b; ind_b_others(matches)];
            % making sure that every blobs is on the same side:
            temp_bounds = blobs_uvw_boundaries(merge_list, :);
            temp_bounds(:, [3 6]) = gtGrainAnglesTabularFix360deg( ...
                temp_bounds(:, [3 6]), blobs_uvw_center_width(ii_b, 3), p);

            blobs_uvw_boundaries(ii_b, 1:3) = min(temp_bounds(:, 1:3), [], 1);
            blobs_uvw_boundaries(ii_b, 4:6) = max(temp_bounds(:, 4:6), [], 1);

            blobs_etas(ii_b) = mean(blobs_etas(merge_list));

            is_selected(ii_b) = any(is_selected(merge_list)) || conf.always_select_shared;
            if (conf.prefer_raw_shared)
                blob_use_raw_images(ii_b) = any(blob_use_raw_images(merge_list));
            else
                blob_use_raw_images(ii_b) = all(blob_use_raw_images(merge_list));
            end

            is_merged(ind_b_others(matches)) = true;
        end
    end
    % Updating structures, to remove merged blobs
    all_blobs_list(is_merged) = [];
    blobs_uvw_boundaries(is_merged, :) = [];
    blob_use_raw_images(is_merged, :) = [];
    blobs_etas(is_merged) = [];
    is_selected(is_merged) = [];
    blobs_uvw_center_width(is_merged, :) = [];
    num_blobs = numel(all_blobs_list);

    fprintf('Determining which blobs are close in UVW space: ')
    c = tic();

    % and now we scan for overlaps of blobs
    is_merged = false(size(all_blobs_list));
    % We will now scan for blobs which have the same ID and merge them
    for ii_b = 1:num_blobs
        num_chars = fprintf('%03d/%03d', ii_b, num_blobs);

        ind_b_others = (ii_b+1:num_blobs)';

        dist_centers_w = gtMathsCircularDistance( ...
            blobs_uvw_center_width(ii_b, 3), ...
            blobs_uvw_center_width(ind_b_others, 3), num_images);
        dist_centers_uv = abs( bsxfun(@minus, ...
            blobs_uvw_center_width(ii_b, 1:2), ...
            blobs_uvw_center_width(ind_b_others, 1:2)) );

        dist_centers = [dist_centers_uv, dist_centers_w];

        sizes_sum = bsxfun(@plus, ...
            blobs_uvw_center_width(ii_b, 4:6), ...
            blobs_uvw_center_width(ind_b_others, 4:6));
        sizes_sum = bsxfun(@plus, sizes_sum, conf.merge_dist_pixels);

        matches = all(dist_centers < sizes_sum, 2) ...
            & reshape(~is_merged(ind_b_others), [], 1);

        if (any(matches))
            merge_list = [ii_b; ind_b_others(matches)];
            % making sure that every blobs is on the same side:
            temp_bounds = blobs_uvw_boundaries(merge_list, :);
            temp_bounds(:, [3 6]) = gtGrainAnglesTabularFix360deg( ...
                temp_bounds(:, [3 6]), blobs_uvw_center_width(ii_b, 3), p);

            blobs_uvw_boundaries(ii_b, 1:3) = min(temp_bounds(:, 1:3), [], 1);
            blobs_uvw_boundaries(ii_b, 4:6) = max(temp_bounds(:, 4:6), [], 1);

            all_blobs_list(ii_b) = merge_blobs(all_blobs_list(merge_list), blobs_uvw_center_width(ii_b, 3), p);

            blobs_etas(ii_b) = mean(blobs_etas(merge_list));

            is_selected(ii_b) = any(is_selected(merge_list)) || conf.always_select_shared;
            if (conf.prefer_raw_shared)
                blob_use_raw_images(ii_b) = any(blob_use_raw_images(merge_list));
            else
                blob_use_raw_images(ii_b) = all(blob_use_raw_images(merge_list));
            end

            is_merged(ind_b_others(matches)) = true;
        end

        fprintf(repmat('\b', [1, num_chars]));
    end
    % Updating structures, to remove merged blobs
    all_blobs_list(is_merged) = [];
    blobs_uvw_boundaries(is_merged, :) = [];
    blob_use_raw_images(is_merged, :) = [];
    is_selected(is_merged) = [];
    blobs_etas(is_merged) = [];
    num_blobs = numel(all_blobs_list);

    % updating center and halfwidth
    blobs_uvw_center_width = [ ...
        blobs_uvw_boundaries(:, 4:6) + blobs_uvw_boundaries(:, 1:3), ...
        blobs_uvw_boundaries(:, 4:6) - blobs_uvw_boundaries(:, 1:3) ] / 2;

    % Fixing blobs' centers to be between 0-360
    blobs_w_c = mod(blobs_uvw_center_width(:, 3), num_images);
    blobs_uvw_boundaries(:, [3 6]) = gtGrainAnglesTabularFix360deg(blobs_uvw_boundaries(:, [3 6]), blobs_w_c, p);

    for ii_b = 1:num_blobs
        all_blobs_list(ii_b).bbwim = gtGrainAnglesTabularFix360deg(all_blobs_list(ii_b).bbwim, blobs_w_c(ii_b), p);
        all_blobs_list(ii_b).mbbw = gtGrainAnglesTabularFix360deg(all_blobs_list(ii_b).mbbw, blobs_w_c(ii_b), p);
    end

    % updating center and halfwidth
    blobs_uvw_center_width = [ ...
        blobs_uvw_boundaries(:, 4:6) + blobs_uvw_boundaries(:, 1:3), ...
        blobs_uvw_boundaries(:, 4:6) - blobs_uvw_boundaries(:, 1:3) ] / 2;

    fprintf('\b\b: Done in %g seconds.\nDetermining which orientation regions project to the blobs: ', toc(c))
    c = tic();

    % Matching all ondet for all the regions, to include and select (if
    % conf.always_select_shared is true) the ones that match
    or_inc_pos = cell(num_grains, 1);

    blobs_props = struct( ...
        'ors', cell(num_blobs, 1), ...
        'ors_ondet', cell(num_blobs, 1) );
    temp_proj = gtFwdSimProjDefinition([], num_grains);

    % allow for some additional tolerance in omega position
    blobs_uvw_match_tol = bsxfun(@plus, blobs_uvw_center_width(:, 4:6), conf.match_tol);

    for ii_g = 1:num_grains
        num_chars = fprintf('%03d/%03d', ii_g, num_grains);

        ondet_bls = grs(ii_g).proj(conf.det_index).ondet;
        ondet_uvw = samp_ors(ii_g).allblobs(conf.det_index).detector.uvw(ondet_bls, :);

        or_inc_pos{ii_g} = zeros(size(ondet_bls));

        for ii_o = 1:numel(ondet_bls)
            dist_centers_w = gtMathsCircularDistance( ...
                ondet_uvw(ii_o, 3), ...
                blobs_uvw_center_width(:, 3), num_images);
            dist_centers_uv = abs( bsxfun(@minus, ...
                ondet_uvw(ii_o, 1:2), ...
                blobs_uvw_center_width(:, 1:2)) );

            dist_centers = [dist_centers_uv, dist_centers_w];

            matches = dist_centers <= blobs_uvw_match_tol;
            matches = all(matches, 2);

            if (any(matches)) % TODO: find better logic (fwd-sim type)
                ii_o
                cands = find(matches)
                dist_centers_scalar = sqrt(sum(dist_centers .^ 2, 2));
                [~, ii_c] = min(dist_centers_scalar(cands));

                ii_b = cands(ii_c);

                or_inc_pos{ii_g}(ii_o) = ii_b;
                blobs_props(ii_b).ors = [blobs_props(ii_b).ors, ii_g];
                blobs_props(ii_b).ors_ondet = [blobs_props(ii_b).ors_ondet, ii_o];
            end
        end

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

    if (conf.always_select_shared)
        for ii_b = 1:num_blobs
            is_selected(ii_b) = is_selected(ii_b) || (numel(blobs_props(ii_b).ors) > 1);
        end
    end

    % Filtering vertical blobs
    eta_ok = acosd(abs(cosd(blobs_etas))) > conf.min_eta;
    is_selected = is_selected & eta_ok;

    for ii_g = 1:num_grains
        ondet_bls = grs(ii_g).proj(conf.det_index).ondet;
        valid_inc_pos = or_inc_pos{ii_g} > 0;
        temp_proj(ii_g).ondet = ondet_bls;
        temp_proj(ii_g).included = find(valid_inc_pos);
        temp_proj(ii_g).selected = is_selected(or_inc_pos{ii_g}(valid_inc_pos));
        temp_proj(ii_g).global_pos = or_inc_pos{ii_g}(valid_inc_pos);
    end

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

    for ii_b = 1:num_blobs
        fprintf(' - Blob %03d [u: %4d-%4d, v: %4d-%4d, w %4d-%4d], (is sel: %d, eta ok: %d, use raw: %d) orientation BBs:%s\n', ...
            ii_b, blobs_uvw_boundaries(ii_b, [1 4 2 5 3 6]), is_selected(ii_b), eta_ok(ii_b), blob_use_raw_images(ii_b), sprintf(' %d', blobs_props(ii_b).ors))
    end

    % If conf.use_raw_images is selected, it is now time to act!
    if (any(conf.use_raw_images))
        fprintf('Loading raw images: ')
        c = tic();

        inds_raw_blobs = reshape(find(blob_use_raw_images), 1, []);
        for ii_b = inds_raw_blobs
            num_chars = fprintf('%03d/%03d', ii_b, num_blobs);

            all_blobs_list(ii_b) = load_blob_simple(blobs_uvw_boundaries(ii_b, :), p, conf);

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

    if (any(conf.use_raw_images) && conf.use_volume_mask)
        volume_verts = zeros(0, 3);
        for ii_g = 1:num_grains
            center_shift = grs(ii_g).center - cl_center_mm;
            volume_verts = cat(1, volume_verts, ...
                bsxfun(@plus, grs(ii_g).fwd_sim.gv_verts, center_shift));
        end

        for ii_g = 1:num_grains
            local_global_pos = temp_proj(ii_g).global_pos;
            local_ondet = temp_proj(ii_g).ondet;
            local_included = temp_proj(ii_g).included;

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

                sel = local_ondet(local_included);
                proj_bl_masks = gt6DSpreadProjectVertices2Det(samp_ors(ii_g), volume_verts, sel, p, conf.det_index);

                all_blobs_list = assign_masks(all_blobs_list, proj_bl_masks, local_global_pos);
            end
        end

        for ii_b = 1:numel(all_blobs_list)
            all_blobs_list(ii_b).intensity = sum(all_blobs_list(ii_b).intm(all_blobs_list(ii_b).mask));
        end
    end

    fprintf('\b\b: Done in %g seconds.\nDetermining UVW bounding boxes..', toc(c))
    c = tic();

    bls_bb_sizes = cat(1, all_blobs_list(:).bbsize);
    [stackUSize, stackVSize] = get_stack_UV_size(bls_bb_sizes, blob_use_raw_images, conf);

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

        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_g = 1:num_grains
        fprintf('Processing orientation-space region: %d/%d\n - Producing blobs and stacks..', ii_g, num_grains)

        temp_proj(ii_g).bl = all_blobs_list(temp_proj(ii_g).global_pos);
        temp_proj(ii_g).difstack = zeros(stackUSize, numel(temp_proj(ii_g).bl), stackVSize, 'single');

        for ii_s = 1:numel(temp_proj(ii_g).bl)
            % Essentially zero pads the spot to make it fit into ASTRA's
            % diffraction stack
            spot_pad = sum(temp_proj(ii_g).bl(ii_s).intm, 3);
            temp_proj(ii_g).stack(:, ii_s, :) = spot_pad;
        end

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

    for ii_g = 1:num_grains
        fprintf('Processing orientation-space region: %d/%d\n - Producing proj data-structure..', ii_g, num_grains)
        c = tic();

        proj = gtFwdSimProjDefinition();

        proj.centerpix = cl_center_pix;
        proj.bl = temp_proj(ii_g).bl;
        proj.stack = temp_proj(ii_g).stack;
        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 = temp_proj(ii_g).ondet;
        proj.included = temp_proj(ii_g).included;
        proj.selected = temp_proj(ii_g).selected;

        % Indices that indicate where the blobs are in the list
        proj.global_pos = temp_proj(ii_g).global_pos;

        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}, ...
        'blobs_props', {blobs_props}, ...
        '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, blob_use_raw_images, conf)

    max_img_sizes = max(img_sizes(~blob_use_raw_images, 1:2), 1);
    max_img_sizes = max([ceil(max_img_sizes * conf.stack_oversize); img_sizes(blob_use_raw_images, 1:2)], [], 1);
    u_size = round(max_img_sizes(1));
    v_size = round(max_img_sizes(2));

    fprintf('\n');
    fprintf('               Maximum images size: [%3d, %3d]\n', max_img_sizes);
    fprintf('Stack images size (oversize: %1.2f): [%3d, %3d]\n', conf.stack_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 blob = load_blob_simple(img_bboxes, p, conf)
    blob = gtFwdSimBlobDefinition();

    img_sizes = img_bboxes(1, 4:6) - img_bboxes(1, 1:3) + 1;

    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 = img_sizes;
    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.intm = single(blob_vol);
    if (conf.use_volume_mask)
        blob.mask = false(img_sizes);
    else
        blob.mask = true(img_sizes);
    end

    blob.bbsize = img_sizes;
    blob.bbuim = blob.mbbu;
    blob.bbvim = blob.mbbv;
    blob.bbwim = blob.mbbw;

    blob.intensity = gtMathsSumNDVol(blob.intm);
end

function blob = merge_blobs(bls, ref_w, p)
    blob = gtFwdSimBlobDefinition();

    bls_bb_u = cat(1, bls(:).bbuim);
    bls_bb_v = cat(1, bls(:).bbvim);
    bls_bb_w = cat(1, bls(:).bbwim);
    bls_bb_w = gtGrainAnglesTabularFix360deg(bls_bb_w, ref_w, p);

    blob.bbuim = [min(bls_bb_u(:, 1)), max(bls_bb_u(:, 2))];
    blob.bbvim = [min(bls_bb_v(:, 1)), max(bls_bb_v(:, 2))];
    blob.bbwim = [min(bls_bb_w(:, 1)), max(bls_bb_w(:, 2))];

    blob.bbsize = [ ...
        blob.bbuim(2) - blob.bbuim(1) + 1, ...
        blob.bbvim(2) - blob.bbvim(1) + 1, ...
        blob.bbwim(2) - blob.bbwim(1) + 1 ];

    bls_mbb_u = cat(1, bls(:).mbbu);
    bls_mbb_v = cat(1, bls(:).mbbv);
    bls_mbb_w = cat(1, bls(:).mbbw);
    bls_mbb_w = gtGrainAnglesTabularFix360deg(bls_mbb_w, ref_w, p);

    blob.mbbu = [min(bls_mbb_u(:, 1)), max(bls_mbb_u(:, 2))];
    blob.mbbv = [min(bls_mbb_v(:, 1)), max(bls_mbb_v(:, 2))];
    blob.mbbw = [min(bls_mbb_w(:, 1)), max(bls_mbb_w(:, 2))];

    blob.mbbsize = [ ...
        blob.mbbu(2) - blob.mbbu(1) + 1, ...
        blob.mbbv(2) - blob.mbbv(1) + 1, ...
        blob.mbbw(2) - blob.mbbw(1) + 1 ];

    shifts = [ ...
        bls_bb_u(:, 1) - blob.bbuim(1), ...
        bls_bb_v(:, 1) - blob.bbvim(1), ...
        bls_bb_w(:, 1) - blob.bbwim(1) ];
    blob.intm = zeros(blob.bbsize, 'like', bls(1).intm);
    blob.mask = zeros(blob.bbsize, 'like', bls(1).mask);

    for ii_b = 1:numel(bls)
        blob.intm = gtPlaceSubVolume(blob.intm, bls(ii_b).intm, shifts(ii_b, :));
        blob.mask = gtPlaceSubVolume(blob.mask, bls(ii_b).mask, shifts(ii_b, :));
    end

    blob.intensity = sum([bls(:).intensity]);
end

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

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

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

    for ii = 1:numel(global_pos)
        ii_b = global_pos(ii);

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

        blobs(ii_b).mask = gtPlaceSubVolume(blobs(ii_b).mask, mask, shifts(ii, :));
    end
end