Def: Added function to compute a 6D ODF from distribution of orientations

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Deformation/gtGetODF6DFromGvdm.m

+function odf_6D = gtGetODF6DFromGvdm(gvdm, gvcs, ospace_grid_points, rspace_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 (size(gvcs, 1) == 3)
+        gvcs = gvcs';
+    end
+
+    if (~exist('gvint', 'var'))
+        gvint = ones(size(gvdm, 1), 1);
+    else
+        gvint = reshape(gvint, [], 1);
+    end
+
+    if (iscell(ospace_grid_points))
+        o_voxel_size = ospace_grid_points{2, 2, 2}.R_vector - ospace_grid_points{1, 1, 1}.R_vector;
+
+        size_odf = size(ospace_grid_points);
+
+        o_z_edges = linspace(ospace_grid_points{1, 1, 1}.R_vector(3)-o_voxel_size(3)/2, ...
+            ospace_grid_points{1, 1, end}.R_vector(3)+o_voxel_size(3)/2, size_odf(3)+1);
+        o_y_edges = linspace(ospace_grid_points{1, 1, 1}.R_vector(2)-o_voxel_size(2)/2, ...
+            ospace_grid_points{1, end, 1}.R_vector(2)+o_voxel_size(2)/2, size_odf(2)+1);
+        o_x_edges = linspace(ospace_grid_points{1, 1, 1}.R_vector(1)-o_voxel_size(1)/2, ...
+            ospace_grid_points{end, 1, 1}.R_vector(1)+o_voxel_size(1)/2, size_odf(1)+1);
+    else
+        % We expect a <1x9> vector with <1x6> o-space bb, plus <1x3> the
+        % resolution
+        o_voxel_size = ospace_grid_points(7:9);
+
+        size_odf = round((ospace_grid_points(4:6) - ospace_grid_points(1:3)) ./ o_voxel_size);
+
+        o_z_edges = linspace(ospace_grid_points(3) - o_voxel_size(3) / 2, ...
+            ospace_grid_points(6) + o_voxel_size(3) / 2, size_odf(3) + 1);
+        o_y_edges = linspace(ospace_grid_points(2) - o_voxel_size(2) / 2, ...
+            ospace_grid_points(5) + o_voxel_size(2) / 2, size_odf(2) + 1);
+        o_x_edges = linspace(ospace_grid_points(1) - o_voxel_size(1) / 2, ...
+            ospace_grid_points(4) + o_voxel_size(1) / 2, size_odf(1) + 1);
+    end
+
+    % We expect a <1x9> vector with <1x6> r-space bb, plus <1x3> the resolution
+    r_voxel_size = rspace_grid_points(7:9);
+
+    size_vol = round((rspace_grid_points(4:6) - rspace_grid_points(1:3)) ./ r_voxel_size);
+
+    r_z_edges = linspace(rspace_grid_points(3) - r_voxel_size(3) / 2, ...
+        rspace_grid_points(6) + r_voxel_size(3) / 2, size_vol(3) + 1);
+    r_y_edges = linspace(rspace_grid_points(2) - r_voxel_size(2) / 2, ...
+        rspace_grid_points(5) + r_voxel_size(2) / 2, size_vol(2) + 1);
+    r_x_edges = linspace(rspace_grid_points(1) - r_voxel_size(1) / 2, ...
+        rspace_grid_points(4) + r_voxel_size(1) / 2, size_vol(1) + 1);
+
+    odf_6D = zeros([size_odf, size_vol]);
+
+    fprintf('Producing ODF..');
+    c = tic();
+    for rz_ii = 1:size_vol(3)
+        included_rz = (r_z_edges(rz_ii) <= gvcs(:, 3)) & (gvcs(:, 3) <= r_z_edges(rz_ii+1));
+        gvcs_inc_rz = gvcs(included_rz, :);
+        gvdm_inc_rz = gvdm(included_rz, :);
+        gvint_inc_rz = gvint(included_rz);
+
+        for ry_ii = 1:size_vol(2)
+            included_ry = (r_y_edges(ry_ii) <= gvcs_inc_rz(:, 2)) & (gvcs_inc_rz(:, 2) <= r_y_edges(ry_ii+1));
+            gvcs_inc_ry = gvcs_inc_rz(included_ry, :);
+            gvdm_inc_ry = gvdm_inc_rz(included_ry, :);
+            gvint_inc_ry = gvint_inc_rz(included_ry);
+
+            for rx_ii = 1:size_vol(1)
+                included_rx = (r_x_edges(rx_ii) <= gvcs_inc_ry(:, 1)) & (gvcs_inc_ry(:, 1) <= r_x_edges(rx_ii+1));
+                gvdm_inc_rx = gvdm_inc_ry(included_rx, :);
+                gvint_inc_rx = gvint_inc_ry(included_rx);
+
+                for oz_ii = 1:size_odf(3)
+                    included_oz = (o_z_edges(oz_ii) <= gvdm_inc_rx(:, 3)) & (gvdm_inc_rx(:, 3) <= o_z_edges(oz_ii+1));
+                    gvdm_inc_oz = gvdm_inc_rx(included_oz, :);
+                    gvint_inc_oz = gvint_inc_rx(included_oz);
+
+                    for oy_ii = 1:size_odf(2)
+                        included_oy = (o_y_edges(oy_ii) <= gvdm_inc_oz(:, 2)) & (gvdm_inc_oz(:, 2) <= o_y_edges(oy_ii+1));
+                        gvdm_inc_oy = gvdm_inc_oz(included_oy, :);
+                        gvint_inc_oy = gvint_inc_oz(included_oy);
+
+                        for ox_ii = 1:size_odf(1)
+                            included_ox = (o_x_edges(ox_ii) <= gvdm_inc_oy(:, 1)) & (gvdm_inc_oy(:, 1) <= o_x_edges(ox_ii+1));
+                            gvint_inc_ox = gvint_inc_oy(included_ox);
+
+                            odf_6D(ox_ii, oy_ii, oz_ii, rx_ii, ry_ii, rz_ii) = sum(gvint_inc_ox);
+                        end
+                    end
+                end
+            end
+        end
+    end
+    fprintf('\b\b (%3.1f s) Done.\n', toc(c));
+end
\ No newline at end of file