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

    size_odf = size(grid_points);
    odf = zeros(size_odf);

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

    voxel_size = grid_points{2, 2, 2}.R_vector - grid_points{1, 1, 1}.R_vector;

    z_edges = linspace(grid_points{1, 1, 1}.R_vector(3)-voxel_size(3)/2, ...
        grid_points{1, 1, end}.R_vector(3)+voxel_size(3)/2, size_odf(3)+1);
    y_edges = linspace(grid_points{1, 1, 1}.R_vector(2)-voxel_size(2)/2, ...
        grid_points{1, end, 1}.R_vector(2)+voxel_size(2)/2, size_odf(2)+1);
    x_edges = linspace(grid_points{1, 1, 1}.R_vector(1)-voxel_size(1)/2, ...
        grid_points{end, 1, 1}.R_vector(1)+voxel_size(1)/2, size_odf(1)+1);

    fprintf('Producing ODF..');
    c = tic();
    for l = 1:size_odf(3)
        included3 = (z_edges(l) <= gvdm(:, 3)) & (gvdm(:, 3) <= z_edges(l+1));
        included3 = gvdm(included3, :);
        for n = 1:size_odf(2)
            included2 = (y_edges(n) <= included3(:, 2)) & (included3(:, 2) <= y_edges(n+1));
            included2 = included3(included2, :);
            for m = 1:size_odf(1)
                included1 = (x_edges(m) <= included2(:, 1)) & (included2(:, 1) <= x_edges(m+1));

                odf(m, n, l) = numel(find(included1));
            end
        end
    end
    fprintf('\b\b (%3.1f s) Done.\n', toc(c));
end