Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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