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);
    if (isfield(parameters.rec, 'grains') ...
            && isfield(parameters.rec.grains, 'options') ...
            && ~isempty(parameters.rec.grains.options) ...
            && (~isfield(parameters.rec.grains.options, 'use_predicted_scatter_ints') ...
                || isempty(parameters.rec.grains.options.use_predicted_scatter_ints)) )
        rec_opts.use_predicted_scatter_ints = true;
    end

    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.ospace_super_sampling > 1)
        sampler.make_supersampling_simple_grid([1 2 3], rec_opts.ospace_super_sampling);
    end

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

    vols = algo.getCurrentSolution();

    or_sizes = sampler.get_orientation_sampling_size();

    [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), or_sizes{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);

    % Restoring initial volume size (depending on the rounding)
    if (rec_opts.volume_downscaling > 1)
        avg_R_vecs_int = gtMathsUpsampleVolume(avg_R_vecs_int, rec_opts.volume_downscaling);
        avg_R_vecs = gtMathsUpsampleVolume(avg_R_vecs, rec_opts.volume_downscaling);
        kam = gtMathsUpsampleVolume(kam, rec_opts.volume_downscaling);
        igm = gtMathsUpsampleVolume(igm, rec_opts.volume_downscaling);
    end

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

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

        % 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

        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