gtDefComputeGrainODFNearField.m

function [odfw, odfw_R_vectors] = gtDefComputeGrainODFNearField(phase_id, grain_id, parameters, varargin)
    if (~exist('parameters', 'var') || isempty(parameters))
        parameters = gtLoadParameters();
    end

    conf = struct( ...
        'save', true, ...
        'compare_ODF6D', false, ...
        'verbose', false);

    conf = parse_pv_pairs(conf, varargin);

    sol = GtGrainODFwSolver(parameters, 'verbose', conf.verbose);

    phase_dir = fullfile(parameters.acq.dir, '4_grains', ...
        sprintf('phase_%02d', phase_id));
    gr = gtLoadGrain(phase_id, grain_id);

    odfw = sol.solve(gr);

    odfw_R_vectors = sol.get_R_vectors();

    if (conf.compare_ODF6D || conf.save)
        grain_det_file = fullfile(phase_dir, sprintf('grain_details_%04d.mat', grain_id));
        gr_det = load(grain_det_file);
    end

    if (conf.save)
        gr_det.ODFw = struct(...
            'single_grain_ODF', {odfw}, ...
            'R_vectors', {odfw_R_vectors} );

        save(grain_det_file, '-struct', 'gr_det', '-v7.3');
    end

    if (conf.compare_ODF6D)
        if (~conf.verbose)
            fprintf([ 'Overriding verbose=false (but not too much) ' ...
                'because a comparison with ODF6D was asked.\n' ])
        end

        % The orientation information comes as a DMVOL structure, and we want
        % it to be in GVDM format.
        odfw_orientations = sol.get_orientations();
        gvdm6D = gtDefDmvol2Gvdm(gr_det.ODF6D.voxels_avg_R_vectors);
        odf6D = gtGetODFFromGvdm(gvdm6D, odfw_orientations, gr_det.ODF6D.intensity);

        odfw_coeffs = reshape(odfw, [], 1);
        odfw_coeffs = odfw_coeffs(:, [1 1 1]);

        avg_orient_ODFw = sum(odfw_R_vectors .* odfw_coeffs, 1) ./ sum(odfw_coeffs, 1);

        fprintf([ 'Average orientations:\n' ...
            ' - INDXT : (%f, %f, %f)\n' ...
            ' - ODFw  : (%f, %f, %f)\n' ...
            ' - ODF6D : (%f, %f, %f)\n' ], ...
            gr.R_vector, avg_orient_ODFw, gr_det.ODF6D.single_grain_avg_R_vector);

        if (conf.verbose)
            f = figure();
            ax = axes('parent', f);
            hold(ax, 'on')
            scatter3(ax, gr.R_vector(1), gr.R_vector(2), gr.R_vector(3), 30);
            scatter3(ax, avg_orient_ODFw(1), avg_orient_ODFw(2), avg_orient_ODFw(3), 30, 'r');
            scatter3(ax, gr_det.ODF6D.single_grain_avg_R_vector(1), ...
                gr_det.ODF6D.single_grain_avg_R_vector(2), ...
                gr_det.ODF6D.single_grain_avg_R_vector(3), 30, 'g');
            hold(ax, 'off')

            GtVolView.compareVolumes({odf6D, odfw})
        end
    end
end