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, ...
        '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)
        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);
        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), ...
            grs(1).allblobs.hklsp(vis_sel, [1 2 4]), ...
            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

    if (conf.verbose && false)
        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(phase_id), 'int_exp'))
            predicted_ints = p.cryst(phase_id).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(phase_id).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

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

    if (conf.verbose)
        fprintf('\n');
        fprintf('Estimated spatial voxel BBox: [%3d, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix);
        fprintf('                   BBox size: %3d, %3d, %3d (%f, %f, %f mm)\n', bbox_size_pix, bbox_size_mm);
    end

    img_sizes = zeros(num_grains, 2);

    for ii_g = num_grains:-1:1
        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);

        if (conf.verbose)
            fprintf('%d) Estimated orientation BBox: [%3.3f, %3.3f, %3.3f] -> [%3.3f, %3.3f, %3.3f]\n', ii_g, estim_orient_bbox);
            fprintf('                    BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
            fprintf('\n');
        end

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

        img_sizes(ii_g, :) = [ ...
            size(grs(ii_g).proj(conf.det_index).stack, 1), ...
            size(grs(ii_g).proj(conf.det_index).stack, 3), ];

        samp_ors(ii_g) = refor;
    end

    [stackUSize, stackVSize] = get_stack_UV_size(cat(1, img_sizes), p, conf);

    bool_ref_inc = false(size(refgr.allblobs.omega));
    bool_ref_inc(ref_ondet(ref_included)) = true;
    bool_ref_sel = false(size(refgr.allblobs.omega));
    bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;

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

    for ii_g = 1:num_grains
        fprintf('Loading raw images for grain %d: ', ii_g)

        if (ii_g > 1)
            twingr = grs(ii_g);

            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(refgr.allblobs.omega));
            bool_twin_ondet(twin_ondet) = true;
            bool_twin_inc = false(size(refgr.allblobs.omega));
            bool_twin_inc(twin_ondet(twin_inc)) = true;
            bool_twin_sel = false(size(refgr.allblobs.omega));
            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;

            if (false)
                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);
            else
                test_shared = refl_matches(bool_test_inc, ii_g-1);
            end

            % let's add the shared as included
            shared_parent = false(size(refgr.allblobs.omega));
            bool_test_inc_indx = find(bool_test_inc);
            shared_parent(bool_test_inc_indx(test_shared)) = true;

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

            % We now find the blobs in the old twin bl structure that are
            % shared with the parent
            shared_parent_inc_redundant = 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_parent_pos = bool_ref_inc & shared_parent;
            shared_parent_pos = shared_parent_pos(ref_ondet(ref_included));

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

            local_ondet = twin_ondet;
            local_included = twin_inc;
            local_selected = twin_sel;
        else
            local_ondet = ref_ondet;
            local_included = ref_included;
            local_selected = ref_selected;
        end

        num_blobs = numel(local_included);

        if (ii_g == 1)
            blobs = grs(ii_g).proj.bl;
        else
            blobs = gtFwdSimBlobDefinition('blob', num_blobs);
            blobs(~shared_parent) = grs(ii_g).proj.bl(~shared_parent_inc_redundant);
            % Using the same because it gives the same size/uvw-limits
            blobs(shared_parent) = samp_ors(1).proj.bl(shared_parent_pos);
        end

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

            blobs(ii_b) = include_blob(blobs(ii_b), stackUSize, stackVSize);

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

        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 + 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 = local_ondet;
        proj.included = local_included;
        proj.selected = local_selected;

        if (ii_g > 1)
            proj.shared_bls_parent = zeros(size(shared_parent));
            proj.shared_bls_parent(shared_parent) = find(shared_parent_pos);

            samp_ors(1).proj(conf.det_index).shared_bls_twins(shared_parent_pos) = true;
        else
            proj.shared_bls_twins = false(num_blobs, 1);
        end

        samp_ors(ii_g).proj(conf.det_index) = proj;

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

        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

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

    max_img_sizes = [max(img_sizes(:, 1)), max(img_sizes(:, 2))];
    stack_imgs_oversize = min(p.fsim.oversize, conf.stack_oversize);
    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 = include_blob(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;

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

    % only included spots info
    gr.cands     = twin.fwd_sim.cands;
    gr.likelihood = twin.fwd_sim.likelihood;
    gr.difspotID = twin.fwd_sim.difspotID;
    gr.uv_exp    = twin.fwd_sim.uvw(:, 1:2);
    gr.om_exp    = twin.fwd_sim.uvw(:, 3);
    gr.bb_exp    = twin.fwd_sim.bb;
    gr.intensity = twin.fwd_sim.intensity;
end