ODF creation: fixed a few bugs

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

zUtil_Deformation/gtGetODF6DFromGvdm.m

@@ -27,37 +27,42 @@ function odf_6D = gtGetODF6DFromGvdm(gvdm, gvcs, ospace_grid_points, rspace_grid
     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_R_vector = ospace_grid_points{1, 1, 1}.R_vector;
-        max_R_vector = ospace_grid_points{end, end, end}.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_R_vector = ospace_grid_points(1:3);
-        max_R_vector = ospace_grid_points(4:6);
+        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_R_vector - o_voxel_size/2;
-    max_o_corner = max_R_vector + o_voxel_size/2;
+    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) ./ 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;
 
     % We expect a <1x9> vector with <1x6> r-space bb, plus <1x3> the resolution
     r_voxel_size = rspace_grid_points(7:9);
+    half_r_voxel_size = r_voxel_size / 2;
+    inv_r_voxel_size = 1 ./ r_voxel_size;
 
-    min_r_corner = rspace_grid_points(1:3) - r_voxel_size/2;
-    max_r_corner = rspace_grid_points(4:6) + r_voxel_size/2;
+    min_r_voxel_center = rspace_grid_points(1:3);
+    max_r_voxel_center = rspace_grid_points(4:6);
 
-    size_vol = round((max_r_corner - min_r_corner) ./ r_voxel_size);
+    min_r_corner = min_r_voxel_center - half_r_voxel_size;
+    max_r_corner = max_r_voxel_center + half_r_voxel_size;
 
-    range_o_corners = (size_odf - 1) ./ (max_o_corner - min_o_corner);
-    gvdm_rescaled = gvdm - min_o_corner(ones(num_voxels, 1), :);
-    gvdm_rescaled = gvdm_rescaled .* range_o_corners(ones(num_voxels, 1), :) + 1;
+    size_vol = round((max_r_corner - min_r_corner) .* inv_r_voxel_size);
 
-    range_r_corners = (size_vol - 1) ./ (max_r_corner - min_r_corner);
-    gvcs_rescaled = gvcs - min_r_corner(ones(num_voxels, 1), :);
-    gvcs_rescaled = gvcs_rescaled .* range_r_corners(ones(num_voxels, 1), :) + 1;
+    gvcs_rescaled = gvcs - min_r_voxel_center(ones(num_voxels, 1), :);
+    gvcs_rescaled = gvcs_rescaled .* inv_r_voxel_size(ones(num_voxels, 1), :) + 1;
 
     fprintf('Producing ODF..');
     c = tic();
@@ -68,26 +73,32 @@ function odf_6D = gtGetODF6DFromGvdm(gvdm, gvcs, ospace_grid_points, rspace_grid
 
             odf_6D = accumarray([inds_o inds_r], gvint, [size_odf, size_vol]);
         case 'linear'
-            [inds8_o, ints8_o] = gtMathsGetInterpolationIndices(gvdm_rescaled, gvint);
-            [inds8_r, ints8_r] = gtMathsGetInterpolationIndices(gvcs_rescaled, ones(size(gvint)));
+            ones_gvint = gtMathsGetSameSizeOnes(gvint);
+            [inds8_o, ints8_o] = gtMathsGetInterpolationIndices(gvdm_rescaled, ones_gvint);
+            [inds8_r, ints8_r] = gtMathsGetInterpolationIndices(gvcs_rescaled, ones_gvint);
 
             % inds8 is <(nx8) x 3> while ints8 is <(nx8) x 1>
 
-            ints8_r = permute(reshape(ints8_r, [], 8), [2 3 1]);
             ints8_o = permute(reshape(ints8_o, [], 8), [3 2 1]);
+            ints8_r = permute(reshape(ints8_r, [], 8), [2 3 1]);
             ints64 = ints8_o(ones(8, 1), :, :) .* ints8_r(:, ones(8, 1), :);
-            ints64 = reshape(ints64, [], 1);
+            gvint = reshape(gvint, 1, []);
+            gvint = gvint(ones(64, 1), :);
+            ints64 = reshape(ints64, [], 1) .* reshape(gvint, [], 1);
 
             inds8_o = reshape(inds8_o, [], 8, 3);
-            inds8_o = permute(inds8_o, [3 4 2 1]);
-
             inds8_r = reshape(inds8_r, [], 8, 3);
+
+            inds8_o = permute(inds8_o, [3 4 2 1]);
             inds8_r = permute(inds8_r, [3 2 4 1]);
+            inds8_o = inds8_o(:, ones(8, 1), :, :);
+            inds8_r = inds8_r(:, :, ones(8, 1), :);
+            inds8_o = reshape(inds8_o, 3, [])';
+            inds8_r = reshape(inds8_r, 3, [])';
 
-            inds64 = [inds8_o(:, ones(8, 1), :, :); inds8_r(:, :, ones(8, 1), :)];
-            inds64 = reshape(inds64, 6, [])';
+            inds64 = [inds8_o, inds8_r];
 
-            valid = ints64 > 0;
+            valid = ints64 > eps('single');
             inds64 = inds64(valid, :);
             ints64 = ints64(valid, :);
 

zUtil_Deformation/gtGetODFFromGvdm.m

@@ -23,25 +23,26 @@ function odf = gtGetODFFromGvdm(gvdm, ospace_grid_points, gvint, mode)
     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_R_vector = ospace_grid_points{1, 1, 1}.R_vector;
-        max_R_vector = ospace_grid_points{end, end, end}.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_R_vector = ospace_grid_points(1:3);
-        max_R_vector = ospace_grid_points(4:6);
+        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_R_vector - o_voxel_size/2;
-    max_o_corner = max_R_vector + o_voxel_size/2;
+    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) ./ o_voxel_size);
+    size_odf = round((max_o_corner - min_o_corner) .* inv_o_voxel_size);
 
-    range_corners = (size_odf - 1) ./ (max_o_corner - min_o_corner);
-    gvdm_rescaled = gvdm - min_o_corner(ones(num_voxels, 1), :);
-    gvdm_rescaled = gvdm_rescaled .* range_corners(ones(num_voxels, 1), :) + 1;
+    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();
@@ -53,7 +54,7 @@ function odf = gtGetODFFromGvdm(gvdm, ospace_grid_points, gvint, mode)
         case 'linear'
             [inds8, ints8] = gtMathsGetInterpolationIndices(gvdm_rescaled, gvint);
 
-            valid = ints8 > 0;
+            valid = ints8 > eps('single');
             inds8 = inds8(valid, :);
             ints8 = ints8(valid, :);