function gtReconstructGrainOrientation(grain_id, phase_id, parameters, varargin)
% gtReconstructGrainOrientation  6D reconstructions on a GPU machine
%     gtAstraReconstructGrain(grainID, phaseID, [parameters])
%     -------------------------------------------------------
    if (~exist('parameters','var') || isempty(parameters))
        parameters = gtLoadParameters();
    end

    conf = struct( ...
        'ospace_resolution', [], ...
        'ospace_lims', [] );
    [conf, ~] = parse_pv_pairs(conf, varargin);

    rec_opts = gtReconstruct6DGetParamenters(parameters);

    sample = GtSample.loadFromFile();
    is_extended = sample.phases{phase_id}.getUseExtended(grain_id);

    fprintf('Loading the grain file..')
    c = tic();
    gr = gtLoadGrain(phase_id, grain_id, 'is_extended', is_extended);

    if (~isfield(gr.proj, 'bl'))
        gr_extra = gtLoadGrain(phase_id, gr.id, 'fields', {'bl'});
        gr.proj.bl = gr_extra.bl;
    end
    if (~isfield(gr.proj, 'selected'))
        fields_to_load = {'ondet', 'included', 'selected'};
        gr_extra = gtLoadGrain(phase_id, gr.id, 'fields', fields_to_load);
        gr.proj.ondet = gr_extra.ondet;
        gr.proj.included = gr_extra.included;
        gr.proj.selected = gr_extra.selected;
    end
    fprintf('\b\b: Done in %g seconds.\n', toc(c))

    if (is_extended)
        bb_ors = gr.bb_ors;
        ospace_bb = cat(1, bb_ors(:).R_vector);
    elseif (~isempty(conf.ospace_lims))
        diff_r_vecs = tand(conf.ospace_lims / 2);
        ospace_bb = [gr.R_vector + diff_r_vecs(1:3); gr.R_vector + diff_r_vecs(4:6)];
    end

    sampler = GtOrientationSampling(parameters, gr);
    if (exist('ospace_bb', 'var') && ~isempty(ospace_bb))
%         sampler.make_simple_grid('cubic', rec_opts.grid_edge, ospace_bb', 1)
        if (~isempty(conf.ospace_resolution))
            sampler.make_res_simple_grid('cubic', conf.ospace_resolution, ospace_bb', rec_opts.ospace_oversize);
        else
            sampler.make_even_simple_grid('cubic', rec_opts.grid_edge, ospace_bb', 1)
        end
    else
        sampler.make_simple_grid_estim_ODF('cubic', rec_opts.grid_edge, false, rec_opts.ospace_oversize);
    end
    if (rec_opts.ospace_super_sampling > 1)
        sampler.make_supersampling_simple_grid([1 2 3], rec_opts.ospace_super_sampling);
    end

    algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters);

    vols = algo.getCurrentSolution();

    or_sizes = sampler.get_orientation_sampling_size();

    fprintf('Producing output data-structure..')
    c = tic();
    [avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = sampler.getAverageOrientations(vols);
    avg_R_vec = sampler.getAverageOrientation(vols);
    s_g_odf = reshape(sampler.getODF(vols), or_sizes{1});

    [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);
        stddev_R_vecs = gtMathsUpsampleVolume(stddev_R_vecs, rec_opts.volume_downscaling);
        kam = gtMathsUpsampleVolume(kam, rec_opts.volume_downscaling);
        igm = gtMathsUpsampleVolume(igm, rec_opts.volume_downscaling);
    end

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

    ODF6D = gtReconstructDataStructureDefinition('VOL6D');

    ODF6D.voxels_avg_R_vectors = avg_R_vecs;
    ODF6D.intensity = avg_R_vecs_int;
    ODF6D.shift = shift;
    ODF6D.R_vectors = sampler.get_R_vectors();
    ODF6D.voxels_stddev_R_vectors = stddev_R_vecs;
    ODF6D.single_grain_ODF = s_g_odf;
    ODF6D.single_grain_avg_R_vector = avg_R_vec;
    ODF6D.kernel_average_misorientation = kam;
    ODF6D.intra_granular_misorientation = igm;

    fprintf('\b\b: Done in %g seconds.\n', toc(c))

    fprintf('Saving the reconstruction file..')
    c = tic();
    % Saving and cleaning at the same time
    gr_rec = struct('ODF6D', ODF6D);
    gtSaveGrainRec(phase_id, gr.id, gr_rec, ...
        'is_extended', is_extended, 'clean', true);
    fprintf('\b\b: Done in %g seconds.\n', toc(c))

    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>

        phase_dir = fullfile(parameters.acq.dir, '4_grains', sprintf('phase_%02d', phase_id));
        if (is_extended)
            grain_full_details_file = fullfile(phase_dir, ...
                sprintf('grain_extended_full_details_%04d.mat', gr.id));
        else
            grain_full_details_file = fullfile(phase_dir, ...
                sprintf('grain_full_details_%04d.mat', gr.id));
        end
        save(grain_full_details_file, 'ODF6D', '-v7.3');
    end
end