gtReconstructGrainCluster.m

function varargout = 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

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

    rec_opts = gtReconstruct6DGetParamenters(parameters);
    if (~isempty(conf.ospace_resolution))
        rec_opts.ospace_resolution = conf.ospace_resolution;
    end

    % handling input, and making sure that things are called by their name
    if (isstruct(grain_ids))
        cl = grain_ids;
    else
        fprintf('Loading the cluster file..')
        c = tic();
        cl = gtLoadCluster(phase_id, grain_ids);
        fprintf('\b\b: Done in %g seconds.\n', toc(c))
    end

    if (isfield(cl, 'samp_ors'))
        % this means it is a twin cluster, but this is the only place where
        % it matters
        if (isfield(cl, 'blobs'))
            blobs = cl.blobs;
        else
            blobs = [];
        end
        cl = cl.samp_ors;
        grain_ids = [cl(:).id];
    else
        grain_ids = cl.gids;
        blobs = [];
    end

    num_sub_ors = numel(cl);
    if ((num_sub_ors > 1) && (numel(rec_opts.ospace_resolution) == 1))
        rec_opts.ospace_resolution(2:num_sub_ors) = rec_opts.ospace_resolution;
    end

    for ii_g = 1:num_sub_ors
        sampler(ii_g) = GtOrientationSampling(parameters, cl(ii_g)); %#ok<AGROW>

        if (~isempty(rec_opts.ospace_resolution))
            if (isempty(conf.ospace_lims) || (size(conf.ospace_lims, 1) < ii_g))
                bb_gvdm = cat(1, cl(ii_g).bb_ors(:).R_vector);
            else
                diff_r_vecs = tand(conf.ospace_lims(ii_g, :) / 2);
                bb_gvdm = [cl(ii_g).R_vector + diff_r_vecs(1:3); cl(ii_g).R_vector + diff_r_vecs(4:6)];
            end
            sampler(ii_g).make_grid_resolution('cubic', ...
                rec_opts.ospace_resolution(ii_g), bb_gvdm', ...
                rec_opts.ospace_oversize, rec_opts.max_grid_edge_points);
        else
            bb_gvdm = cat(1, cl(ii_g).bb_ors(:).R_vector);
            sampler(ii_g).make_grid_numpoints('cubic', ...
                rec_opts.max_grid_edge_points, bb_gvdm', true, ...
                rec_opts.ospace_oversize)
        end

        if (rec_opts.ospace_super_sampling > 1)
            sampler(ii_g).make_supersampling_simple_grid([1 2 3], ...
                rec_opts.ospace_super_sampling);
        end
    end

    algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters, 'blobs', blobs);

    vols = algo.getCurrentSolution();

    ors_ranges = algo.orientation_groups;

    gam = zeros(num_sub_ors, 1);
    gos = zeros(num_sub_ors, 1);

    ODF6D = gtReconstructDataStructureDefinition('VOL6D', num_sub_ors);

    fprintf('Producing output data-structures: ')
    c = tic();
    for ii_g = 1:num_sub_ors
        num_chars = fprintf('%03d/%03d', ii_g, num_sub_ors);

        ODF6D(ii_g).options = rec_opts;
        ODF6D(ii_g).compute_statistics = algo.get_statistics();

        or_sizes = sampler(ii_g).get_orientation_sampling_size();
        gr_vols = vols(ors_ranges(ii_g, 1):ors_ranges(ii_g, 2));

        [avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = sampler(ii_g).getAverageOrientations(gr_vols);
        avg_R_vec = sampler(ii_g).getAverageOrientation(gr_vols);
        s_g_odf = reshape(sampler(ii_g).getODF(gr_vols), or_sizes{1});

        [kam, gam(ii_g)] = gtDefComputeKernelAverageMisorientation(avg_R_vecs, avg_R_vecs_int);
        [igm, gos(ii_g)] = gtDefComputeIntraGranularMisorientation(avg_R_vecs, avg_R_vecs_int, 'R_vector', cl(ii_g).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(cl(ii_g).proj, parameters, vol_size);

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

        fprintf(repmat('\b', [1, num_chars]));
    end
    fprintf('Done (%d) in %g seconds.\n', num_sub_ors, toc(c))

    if (conf.save)
        fprintf('Saving the cluster reconstruction file..')
        c = tic();
        % Saving and cleaning at the same time
        cl_rec = struct('ODF6D', ODF6D);
        gtSaveClusterRec(phase_id, grain_ids, cl_rec, 'clean', true);
        fprintf('\b\b: Done in %g seconds.\n', toc(c))
    end

    if (nargout > 0)
        varargout{1} = ODF6D;
    end

    if (conf.extra_output)
        clear('ODF6D')
        [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

        for ii_g = num_sub_ors:-1:1
            gr_vols = vols(ors_ranges(ii_g, 1):ors_ranges(ii_g, 2));

            ODF6D(ii_g) = struct( ...
                'orientation_volumes', {gr_vols}, ...
                'fwd_projected_blobs', {proj_blobs}, ...
                'fwd_projected_spots', {proj_spots}, ...
                'grain_average_misorientation', {gam(ii_g)}, ...
                'grain_orientation_spread', {gos(ii_g)} );
        end

        if (conf.save)
            phase_dir = fullfile(parameters.acq.dir, '4_grains', ...
                sprintf('phase_%02d', phase_id));
            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
        if (nargout > 1)
            varargout{2} = ODF6D;
        end
    end
end