function vols = 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 = parameters.rec;
    if (isfield(rec, 'grains') && isfield(rec.grains, 'options') ...
            && ~isempty(rec.grains.options))
        rec_opts = rec.grains.options;
        if (~isfield(rec_opts, 'volume_downscaling') ...
                || isempty(rec_opts.volume_downscaling))
            rec_opts.volume_downscaling = 1;
        end
    else
        warning('gtReconstructGrainOrientation:wrong_parameters', ...
            'The rec.grains structure doesn''t seem to be valid. Falling back to defaults')
        rec_opts = struct( ...
            'grid_edge', 7, 'num_interp', 1, 'lambda', 1e-1, ...
            'super_sampling', 1, 'volume_downscaling', 1 );
    end
    if (isfield(rec, 'grains'))
        rec_opts.num_iter = rec.grains.num_iter;
        rec_opts.algorithm = rec.grains.algorithm;
    else
        rec_opts.num_iter = rec.num_iter;
        rec_opts.algorithm = '6DL1';
    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);
    if (rec_opts.super_sampling > 1)
        sampler.make_supersampling_simple_grid([1 2 3], rec_opts.super_sampling);
    end

    rec_factory = Gt6DReconstructionAlgorithmFactory(parameters, ...
        'volume_downscaling', rec_opts.volume_downscaling);
    [algo, good_or] = rec_factory.getSubBlobReconstructionAlgo(sampler, rec_opts.num_interp);

    mem_cons = [ ...
        algo.get_peak_memory_consumption(rec_opts.algorithm), ...
        algo.get_memory_consumption_volumes(), ...
        algo.get_memory_consumption_blobs(), ...
        algo.get_memory_consumption_sinograms(), ...
        algo.get_memory_consumption_astra() ];
    fprintf(['Estimated memory consumption: %d MB (Volumes %d MB, Blobs %d MB, Sinograms %d MB)'...
        ' of which ASTRA''s copies: %d MB\n'], ...
        round(mem_cons / (2^20)))

    switch (upper(rec_opts.algorithm))
        case '6DL1'
            algo.cp_ls(30); % To get close to the feasible region
            algo.cp_l1(rec_opts.num_iter, rec_opts.lambda);
        case '6DLS'
            algo.cp_ls(rec_opts.num_iter);
        case '6DTV'
            algo.cp_ls(30); % To get close to the feasible region
            algo.cp_tv(rec_opts.num_iter);
        case '6DTVL1'
            algo.cp_ls(30); % To get close to the feasible region
            algo.cp_tvl1(rec_opts.num_iter, rec_opts.lambda);
        otherwise
            error('gtReconstructGrainOrientation:wrong_parameter', ...
                'Unknown 6D algorithm: "%s"', rec_opts.algorithm)
    end

    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>

    grain_details_file = fullfile(phase_dir, ...
        sprintf('grain_cluster_details%s.mat', grains_str_ids));
    save(grain_details_file, 'ODF6D', '-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