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

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

    mem_cons = [ ...
        algo.get_peak_memory_consumption(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(algorithm))
        case '6DL1'
            algo.cp_ls(30); % To get close to the feasible region
            algo.cp_l1(num_iter, lambda);
        case '6DLS'
            algo.cp_ls(num_iter);
        case '6DTV'
            algo.cp_ls(30); % To get close to the feasible region
            algo.cp_tv(num_iter);
        case '6DTVL1'
            algo.cp_ls(30); % To get close to the feasible region
            algo.cp_tvl1(num_iter, lambda);
        otherwise
            error('gtReconstructGrainOrientation:wrong_parameter', ...
                'Unknown 6D algorithm: "%s"', algorithm)
    end

    vols = algo.getCurrentSolution();

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