gtGetODFFromGvdm.m

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

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

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

    if (iscell(grid_points))
        voxel_size = grid_points{2, 2, 2}.R_vector - grid_points{1, 1, 1}.R_vector;

        size_odf = size(grid_points);

        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);
    else
        % We expect a <1x9> vector with <1x6> o-space bb, plus <1x3> the
        % resolution
        voxel_size = grid_points(7:9);

        size_odf = round((grid_points(4:6) - grid_points(1:3)) ./ voxel_size);

        z_edges = linspace(grid_points(3) - voxel_size(3) / 2, ...
            grid_points(6) + voxel_size(3) / 2, size_odf(3) + 1);
        y_edges = linspace(grid_points(2) - voxel_size(2) / 2, ...
            grid_points(5) + voxel_size(2) / 2, size_odf(2) + 1);
        x_edges = linspace(grid_points(1) - voxel_size(1) / 2, ...
            grid_points(4) + voxel_size(1) / 2, size_odf(1) + 1);
    end
    odf = zeros(size_odf);

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

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