gtReconstructGrainCluster.m

function gtReconstructGrainCluster(grain_ids, phase_id, parameters, varargin)
% gtReconstructGrainOrientation  6D reconstructions on a GPU machine
%     gtAstraReconstructGrain(grainIDs, phaseID, [parameters])
%     -------------------------------------------------------
    if (~exist('parameters', 'var') || isempty(parameters))
        parameters = gtLoadParameters();
    end
    parameters.fsim.check_spot = false;

    rec_opts = gtReconstruct6DGetParamenters(parameters);

    phase_dir = fullfile(parameters.acq.dir, '4_grains', ...
        sprintf('phase_%02d', phase_id));
    grains_str_ids = sprintf('_%04d', grain_ids);
    grain_cluster_file = fullfile(phase_dir, ...
        sprintf('grain_cluster%s.mat', grains_str_ids));
    if (~exist(grain_cluster_file, 'file'))
        grs = [];
        for ii_g = numel(grain_ids):-1:1
            grain_id = grain_ids(ii_g);
            grs(ii_g) = gtForwardSimulate_v2(grain_id, grain_id, ...
                pwd, phase_id, 'save_grain', false, 'display_figure', false);
        end

        [proj, gr, bb_ors] = gt6DCreateProjDataFromGrainCluster(grs, phase_id, varargin);
        gr.proj = proj;
        gr.bb_ors = bb_ors;
        fprintf('Saving the cluster file..')
        save(grain_cluster_file, '-struct', 'gr', '-v7.3');
        fprintf('\b\b: Done.\n')
    else
        fprintf('Loading the cluster file..')
        gr = load(grain_cluster_file);
        fprintf('\b\b: Done.\n')
        bb_ors = gr.bb_ors;
    end
    bb_gvdm = cat(1, bb_ors(:).R_vector)';

    sampler = GtOrientationSampling(parameters, gr);
%     sampler.make_simple_grid('cubic', rec_opts.grid_edge, bb_gvdm, 1);
    sampler.make_even_simple_grid('cubic', rec_opts.grid_edge, bb_gvdm, 1)
    if (rec_opts.super_sampling > 1)
        sampler.make_supersampling_simple_grid([1 2 3], rec_opts.super_sampling);
    end

    [algo, good_or] = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters);

    vols = algo.getCurrentSolution();

    % Restoring initial volume size (depending on the rounding)
    if (rec_opts.volume_downscaling > 1)
        fprintf('Expanding volumes..')
        c = tic();
        vols = gtMathsUpsampleVolume(vols, rec_opts.volume_downscaling);
        fprintf('\b\b: Done (%f seconds).\n', toc(c))
    end

    [avg_R_vecs, avg_R_vecs_int] = sampler.getAverageOrientations(vols, good_or);
    avg_R_vec = sampler.getAverageOrientation(vols, good_or);
    s_g_odf = reshape(sampler.getODF(vols, good_or), rec_opts.grid_edge([1 1 1]));

    vol_size = size(avg_R_vecs_int);
    shift = gtFwdSimComputeVolumeShifts(gr.proj, parameters, vol_size);

    [kam, gam] = gtDefComputeKernelAverageMisorientation(avg_R_vecs, avg_R_vecs_int);
    [igm, gos] = gtDefComputeIntraGranularMisorientation(avg_R_vecs, avg_R_vecs_int, 'R_vector', gr.R_vector);

    ODF6D = struct( ...
        'voxels_avg_R_vectors', {avg_R_vecs}, ...
        'intensity', {avg_R_vecs_int}, ...
        'shift', {shift}, ...
        'R_vectors', {sampler.get_R_vectors()}, ...
        'single_grain_ODF', {s_g_odf}, ...
        'single_grain_avg_R_vector', {avg_R_vec}, ...
        'kernel_average_misorientation', {kam}, ...
        'intra_granular_misorientation', {igm} ); %#ok<NASGU>

    cluster_rec_file = fullfile(phase_dir, ...
        sprintf('grain_cluster_details%s.mat', grains_str_ids));
    % Saving and cleaning at the same time
    if (exist(cluster_rec_file, 'file'))
        cl_rec = load(cluster_rec_file);
        cl_rec.ODF6D = ODF6D;
    else
        cl_rec = struct('ODF6D', ODF6D); %#ok<NASGU>
    end
    save(cluster_rec_file, '-struct', 'cl_rec', '-v7.3');

    if (algo.verbose)
        [proj_blobs, proj_spots] = algo.getProjectionOfCurrentSolution();

        ODF6D = struct( ...
            'orientation_volumes', {vols}, ...
            'fwd_projected_blobs', {proj_blobs}, ...
            'fwd_projected_spots', {proj_spots}, ...
            'grain_average_misorientation', {gam}, ...
            'grain_orientation_spread', {gos}, ...
            'compute_statistics', {algo.get_statistics()} ); %#ok<NASGU>
        grain_full_details_file = fullfile(phase_dir, ...
            sprintf('grain_cluster_full_details%s.mat', grains_str_ids));
        save(grain_full_details_file, 'ODF6D', '-v7.3');
    end
end