gt6DCreateProjDataFromTwinnedGrainFwdProj.m

function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFromTwinnedGrainFwdProj(parent_id, variants, phase_id, varargin)
% FUNCTION samp_ors = 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_toll', 2, ...
        'rspace_oversize', 1.1, ...
        'ospace_oversize', 1.1, ...
        'check_spots', false, ...
        'stack_oversize', 1.4, ...
        'use_raw_images', false, ...
        'save', false );
    conf = parse_pv_pairs(conf, varargin);

    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.omind;
    for ii_g = num_grains:-1:2
        grs(ii_g) = gtGrainAllblobsFilterOrder(grs(ii_g), refgr_omind);
    end

    % At the moment we only support the configuration: parent + twins
    refl_matches = gt6DGetMatchingReflections(refgr, grs(2:end), conf.angular_toll);

    refgr_proj = refgr.proj(conf.det_index);

    ref_ondet = refgr_proj.ondet;
    ref_included = refgr_proj.included;
    ref_selected = refgr_proj.selected;

    if (conf.verbose > 1)
        produce_matching_reflections_table(grs, conf, refl_matches, ref_ondet, ref_included)
    end

    if (conf.verbose > 1 && false)
        produce_beam_intensity_renorm_figure(refgr, p, ref_ondet, ref_included, ref_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);

    samp_ors = cell(num_grains, 1);

    % Determining orientation-space bounding boxes. They are assumed to be
    % detached from each others.
    for ii_g = 1:num_grains
        fprintf('Determining orientation-space region: %d/%d\n', ii_g, num_grains)

        sampler = GtOrientationSampling(p, grs(ii_g), ...
            'detector_index', conf.det_index, 'verbose', conf.verbose);
        r_vecs = sampler.guess_ODF_BB()';

        estim_orient_bbox = [min(r_vecs, [], 1), max(r_vecs, [], 1)];
        bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);

        % oversizing the orienation a bit
        delta_bbox_size_rod = bbox_size_rod * (conf.ospace_oversize - 1) / 2;
        estim_orient_bbox = estim_orient_bbox + [-delta_bbox_size_rod, delta_bbox_size_rod];
        bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);

        bbox_size_deg = 2 * atand(bbox_size_rod);

        fprintf(' - Estimated extent: [%g, %g, %g] -> [%g, %g, %g]\n', estim_orient_bbox);
        fprintf('          BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);

        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, bbox_size_rod, refgr, p);

        samp_ors{ii_g} = refor;

        % Populating the extended projection data-structures
        if (ii_g > 1)
            extended_projs(ii_g) = find_matches_with_refor(refgr, ...
                grs(ii_g), refl_matches(:, ii_g-1), conf);
        else
            extended_projs(ii_g).ondet = ref_ondet;
            extended_projs(ii_g).included = ref_included;
            extended_projs(ii_g).selected = ref_selected;
        end

        local_ondet = extended_projs(ii_g).ondet;
        local_included = extended_projs(ii_g).included;
        refor_ns = refor.allblobs.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

    samp_ors = [samp_ors{:}];

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

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

    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.detector(conf.det_index).uvw(local_inc_refl, :);
            end

            % Let's treat those blobs at the w edge 360->0
            % (from the sampled orientations perspective)
            refor_ws = refor.allblobs.omega(local_inc_refl) / gtGetOmegaStepDeg(p, conf.det_index);
            uvw_tab{ii_g}(:, :, 3) = gtGrainAnglesTabularFix360deg(uvw_tab{ii_g}(:, :, 3), refor_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)
                refor_ns = refor.allblobs.eta(local_inc_refl);

                fprintf('\nImages for orientation %d:\n', ii_g)
                fprintf('%02d)du %8d, dv %8d, dw %8d, eta: %7.3f <- used: %d, u: [%4d %4d], v: [%4d %4d], w: [%4d %4d]\n', ...
                    [(1:numel(refor_ns))', img_sizes{ii_g}, refor_ns, ~inconvenient_etas{ii_g}, img_bboxes{ii_g}(:, [1 4 2 5 3 6]) ]');
            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(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) = samp_ors(1).proj(conf.det_index).bl(extended_projs(ii_g).shared_refl_pos_in_parent);
            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) = grs(ii_g).proj(conf.det_index).bl(~extended_projs(ii_g).shared_refl_included_pos);
                % Using the same because it gives the same size/uvw-limits
                blobs(extended_projs(ii_g).shared_refl) = samp_ors(1).proj(conf.det_index).bl(extended_projs(ii_g).shared_refl_pos_in_parent);

                if (conf.verbose)
                    parent_ids = grs(1).fwd_sim(conf.det_index).spotid(ref_included(extended_projs(ii_g).shared_refl_pos_in_parent))';
                    twin_ids = grs(ii_g).fwd_sim(conf.det_index).spotid(twin_inc(extended_projs(ii_g).shared_refl))';
                    fprintf('Shared included blobs:\n')
                    fprintf(' - %d: parent blob_id %d, twin blob_id %d\n', ...
                        [1:numel(find(extended_projs(ii_g).shared_refl)); parent_ids; twin_ids])
                end
            end
        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;

        % Temporary solution that allows the parent + twin configuration to
        % work, but it should be replaced by a more advanced global
        % solution for handling shared blobs by different orientation boxes
        if (ii_g > 1)
            proj.shared_bls_parent = zeros(size(extended_projs(ii_g).shared_refl));
            proj.shared_bls_parent(extended_projs(ii_g).shared_refl) = find(extended_projs(ii_g).shared_refl_pos_in_parent);

            samp_ors(1).proj(conf.det_index).shared_bls_twins(extended_projs(ii_g).shared_refl_pos_in_parent) = true;
        else
            proj.shared_bls_twins = false(num_blobs, 1);
        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, grs(ii_g).R_vector(1), grs(ii_g).R_vector(2), grs(ii_g).R_vector(3), 30);
        end
        hold(ax, 'off');

        drawnow();
    end

    if (conf.save > 0)
        str_ids = sprintf('_%04d', grs_list);
        grain_filename = fullfile(p.acq.dir, '4_grains', ...
            sprintf('phase_%02d', phase_id), ...
            sprintf('grain_cluster%s.mat', str_ids));
        fprintf('Saving the cluster file "%s"..', grain_filename)
        save(grain_filename, 'samp_ors', '-v7.3');
        fprintf('\b\b: Done.\n')

        if (conf.save == 1)
            fprintf('Saving to sample.mat..')
            sample = GtSample.loadFromFile();
            sample.phases{phase_id}.clusters(end+1) = ...
                GtPhase.makeCluster(grs_list, cat(1, samp_ors(:).R_vector), estim_space_bbox_pix, [], 1);
            sample.saveToFile();
            fprintf('\b\b: Done.\n')
        end
    end
end

function [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf)

    max_img_sizes = [max(img_sizes(:, 1)), max(img_sizes(:, 2))];
    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);

    if (conf.verbose)
        fprintf('\n');
        fprintf('               Maximum images size: [%3d, %3d]\n', max_img_sizes);
        fprintf('Stack images size (oversize: %1.2f): [%3d, %3d]\n', stack_imgs_oversize, u_size, v_size);
        fprintf('\n');
    end
end

function blob = 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, refl_matches, ref_ondet, ref_included)
    vis_sel = ref_ondet(ref_included);
    w_diffs = abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180);
    if (size(grs(1).allblobs.hklsp, 2) == 4)
        hklsp = grs(1).allblobs.hklsp(vis_sel, [1 2 4]);
    else
        hklsp = grs(1).allblobs.hklsp(vis_sel, :);
    end
    fprintf('Reflections included by the parent:\n')
    fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_ind %d, hkl [%2d %2d %2d] <- matches %d, diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
        [(1:numel(vis_sel))', ...
        grs(1).allblobs.omega(vis_sel), grs(2).allblobs.omega(vis_sel), ...
        grs(1).allblobs.eta(vis_sel), grs(2).allblobs.eta(vis_sel), ...
        grs(1).allblobs.omind(vis_sel), hklsp, refl_matches(vis_sel, 1), ...
        w_diffs, w_diffs < grs(1).allblobs.lorentzfac(vis_sel), ...
        grs(1).allblobs.lorentzfac(vis_sel), ...
        acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
        acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
        ]');

    vis_sel = find(refl_matches(:, 1));
    w_diffs = abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180);
    if (size(grs(1).allblobs.hklsp, 2) == 4)
        hklsp = grs(1).allblobs.hklsp(vis_sel, [1 2 4]);
    else
        hklsp = grs(1).allblobs.hklsp(vis_sel, :);
    end
    fprintf('Reflections shared by parent and first twin:\n')
    fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_ind %d, hkl [%2d %2d %2d] <- matches %d, diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
        [(1:numel(vis_sel))', ...
        grs(1).allblobs.omega(vis_sel), grs(2).allblobs.omega(vis_sel), ...
        grs(1).allblobs.eta(vis_sel), grs(2).allblobs.eta(vis_sel), ...
        grs(1).allblobs.omind(vis_sel), hklsp, refl_matches(vis_sel, 1), ...
        w_diffs, w_diffs < grs(1).allblobs.lorentzfac(vis_sel), ...
        grs(1).allblobs.lorentzfac(vis_sel), ...
        acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
        acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
        ]');
end

function produce_beam_intensity_renorm_figure(refgr, p, ref_ondet, ref_included, ref_selected)
    try
        refgr = gtComputeIncomingBeamIntensity(refgr, p);
        beam_ints = refgr.beam_intensity(ref_selected);
        beam_ints = beam_ints / mean(beam_ints);
    catch mexc
        gtPrintException(mexc, 'Skipping beam intensity renormalization')
        beam_ints = 1;
    end
    blob_ints = refgr.intensity(ref_selected);

    vis_spots = ref_ondet(ref_included(ref_selected));
    thetatypes = refgr.allblobs.thetatype(vis_spots);

    if (isfield(p.cryst(refgr.phaseid), 'int_exp'))
        predicted_ints = p.cryst(refgr.phaseid).int_exp(thetatypes)';

        f = figure();
        ax = axes('parent', f);
        plot(ax, blob_ints)
        hold(ax, 'on')
        plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
        predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
        plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
        predicted_ints = predicted_ints .* beam_ints;
        plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
    end

    predicted_ints = p.cryst(refgr.phaseid).int(thetatypes)';

    f = figure();
    ax = axes('parent', f);
    plot(ax, blob_ints)
    hold(ax, 'on')
    plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
    predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
    plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
    predicted_ints = predicted_ints .* beam_ints;
    plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')

    fprintf('%d) %d, %g\n', ...
        [(1:numel(blob_ints))', refgr.allblobs.thetatype(vis_spots), ...
        blob_ints - predicted_ints / mean(predicted_ints) * mean(blob_ints)]')
    hold(ax, 'off')
end

function 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, ...
        'shared_refl', repl_cell, ...
        'shared_refl_pos_in_parent', repl_cell, ...
        'shared_refl_included_pos', repl_cell );
end

function eproj = find_matches_with_refor(refor, twingr, refl_matches, conf)
    size_refl_list = size(refor.allblobs.omega);

    refor_proj = refor.proj(conf.det_index);

    ref_ondet = refor_proj.ondet;
    ref_included = refor_proj.included;
    ref_selected = refor_proj.selected;

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

    twingr_proj = twingr.proj(conf.det_index);

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

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

    % Indices of reflections that exist on the detector for the
    % twin, and were included for the parent
    bool_test_inc = bool_ref_inc & bool_twin_ondet;

    % Alternative code to the use of refl_matches
%     inc_ref_omegas = samp_ors(1).allblobs.omega(bool_test_inc);
%     inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(bool_test_inc);
%     inc_ref_etas = samp_ors(1).allblobs.eta(bool_test_inc);
%
%     inc_twin_omegas = samp_ors(ii_g).allblobs.omega(bool_test_inc);
%     inc_twin_etas = samp_ors(ii_g).allblobs.eta(bool_test_inc);
%
%     w_diff = abs(mod(inc_ref_omegas - inc_twin_omegas + 180, 360) - 180);
%     n_diff = abs(mod(inc_ref_etas   - inc_twin_etas   + 180, 360) - 180);
%
%     test_shared = ...
%         (w_diff < conf.angular_toll * inc_ref_Lfac) ...
%         & (n_diff < conf.angular_toll);

    test_shared = refl_matches(bool_test_inc);

    % let's add the shared as included
    shared_refl = false(size_refl_list);
    bool_test_inc_indx = find(bool_test_inc);
    shared_refl(bool_test_inc_indx(test_shared)) = true;

    bool_twin_inc_redundant = bool_twin_inc & shared_refl;
    % So if it is shared, we use the parent's blob
    bool_twin_inc = bool_twin_inc | shared_refl;
    % If the parent enables it, we enable it as well
    bool_twin_sel = bool_twin_sel | (shared_refl & bool_ref_sel);

    eproj = get6DExtendedProjDefinition();

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

    twin_inc = find(bool_twin_inc(twin_ondet));
    twin_sel = bool_twin_sel(twin_ondet(twin_inc));

    shared_refl_pos_in_parent = bool_ref_inc & shared_refl;
    eproj.shared_refl_pos_in_parent = shared_refl_pos_in_parent(ref_ondet(ref_included));

    % We cut it down to the size of the included
    eproj.shared_refl = shared_refl(twin_ondet(twin_inc));

    eproj.ondet = twin_ondet;
    eproj.included = twin_inc;
    eproj.selected = twin_sel;
end

function blob = load_blob(img_bboxes, img_sizes, stackUSize, stackVSize, p, conf)
    blob = gtFwdSimBlobDefinition();

    bb = [img_bboxes(1, 1:2), img_sizes(1, 1:2)];
    blob_vol = gtGetRawRoi(img_bboxes(3), img_bboxes(6), p.acq, bb, conf.det_index);
    blob_vol(blob_vol < 0) = 0;
    blob_vol(isnan(blob_vol)) = 0;
    blob_bb = [img_bboxes(1, 1:3), img_sizes(1, :)];

    % Transposing to keep the same convention as spots
    blob_vol = permute(blob_vol, [2 1 3]);

    blob.mbbsize = blob_bb(4:6);
    blob.mbbu = [blob_bb(1), blob_bb(1) + blob_bb(4) + 1];
    blob.mbbv = [blob_bb(2), blob_bb(2) + blob_bb(5) + 1];
    blob.mbbw = [blob_bb(3), blob_bb(3) + blob_bb(6) + 1];

    blob_size_im = [stackUSize, stackVSize, blob_bb(6)+2];

    shifts_blob = gtFwdSimGetStackShifts(stackUSize, stackVSize, blob_bb, false);
    shifts_blob = [shifts_blob.u, shifts_blob.v, 1];

    blob.intm = gtPlaceSubVolume( ...
        zeros(blob_size_im, 'single'), single(blob_vol), shifts_blob);

    blob.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