gtGetODFFromGvdm.m

function odf = gtGetODFFromGvdm(gvdm, ospace_grid_points, gvint, mode)
% FUNCTION odf = gtGetODFFromGvdm(gvdm, ospace_grid_points, gvint, mode)
%
% ddm is a <n_voxels x 3> vector, where each row is the Rodriguez vector of
% each voxel

    if (~exist('mode', 'var') || isempty(mode))
        mode = 'linear';
    end

    if (size(gvdm, 1) == 3)
        gvdm = gvdm';
    end

    if (~exist('gvint', 'var'))
        gvint = ones(size(gvdm, 1), 1);
    else
        gvint = reshape(gvint, [], 1);
    end

    num_voxels = numel(gvint);

    if (iscell(ospace_grid_points))
        o_voxel_size = ospace_grid_points{2, 2, 2}.R_vector - ospace_grid_points{1, 1, 1}.R_vector;

        min_o_voxel_center = ospace_grid_points{1, 1, 1}.R_vector;
        max_o_voxel_center = ospace_grid_points{end, end, end}.R_vector;
    else
        % We expect a <1x9> vector with <1x6> o-space bb, plus <1x3> the
        % resolution
        o_voxel_size = ospace_grid_points(7:9);

        min_o_voxel_center = ospace_grid_points(1:3);
        max_o_voxel_center = ospace_grid_points(4:6);
    end
    half_o_voxel_size = o_voxel_size / 2;
    inv_o_voxel_size = 1 ./ o_voxel_size;

    min_o_corner = min_o_voxel_center - half_o_voxel_size;
    max_o_corner = max_o_voxel_center + half_o_voxel_size;

    size_odf = round((max_o_corner - min_o_corner) .* inv_o_voxel_size);

    gvdm_rescaled = gvdm - min_o_voxel_center(ones(num_voxels, 1), :);
    gvdm_rescaled = gvdm_rescaled .* inv_o_voxel_size(ones(num_voxels, 1), :) + 1;

    fprintf('Producing ODF..');
    c = tic();
    switch (mode)
        case 'nearest'
            inds = round(gvdm_rescaled);
            ints = gvint;
        case 'linear'
            [inds, ints] = gtMathsGetInterpolationIndices(gvdm_rescaled, gvint);
    end

    valid = ints > eps('single') & all(inds > 0, 2) ...
        & inds(:, 1) <= size_odf(1) & inds(:, 2) <= size_odf(2) & inds(:, 3) <= size_odf(3);
    inds = inds(valid, :);
    ints = ints(valid, :);

    odf = accumarray(inds, ints, size_odf);

    fprintf('\b\b (%3.1f s) Done. Total intensity: %g, included: %g\n', ...
        toc(c), sum(gvint), sum(ints));
end