Fix CI tests for #2602

bin/dwifslpreproc

@@ -1140,31 +1140,51 @@ def execute(): #pylint: disable=unused-variable
     bvecs = matrix.load_matrix(bvecs_path)
     bvecs_combined_transpose = [ ]
     bvals_combined = [ ]
+    nans_present_bvecs = 0
 
     for pair in volume_pairs:
       bvec_mean = [ 0.5*(bvecs[0][pair[0]] + bvecs[0][pair[1]]),
                     0.5*(bvecs[1][pair[0]] + bvecs[1][pair[1]]),
                     0.5*(bvecs[2][pair[0]] + bvecs[2][pair[1]]) ]
-      norm2 = matrix.dot(bvec_mean, bvec_mean)
-
-      # If one diffusion sensitisation gradient direction is reversed with respect to
-      #   the other, still want to enable their recombination; but need to explicitly
-      #   account for this when averaging the two directions
-      if norm2 < 0.5:
-        bvec_mean = [ 0.5*(bvecs[0][pair[0]] - bvecs[0][pair[1]]),
-                      0.5*(bvecs[1][pair[0]] - bvecs[1][pair[1]]),
-                      0.5*(bvecs[2][pair[0]] - bvecs[2][pair[1]]) ]
-        norm2 = matrix.dot(bvec_mean, bvec_mean)
-
-      # Occasionally a b=0 volume can have a zero vector
-      if norm2:
-        factor = 1.0 / math.sqrt(norm2)
-        new_vec = [ bvec_mean[0]*factor, bvec_mean[1]*factor, bvec_mean[2]*factor ]
+      if any(math.isnan(value) for value in bvec_mean):
+        nans_present_bvecs += 1
+        new_bvec = [0.0, 0.0, 0.0]
+        new_bval = 0.0
       else:
-        new_vec = [ 0.0, 0.0, 0.0 ]
-      bvecs_combined_transpose.append(new_vec)
-      bvals_combined.append(0.5 * (grad[pair[0]][3] + grad[pair[1]][3]))
+        norm2 = matrix.dot(bvec_mean, bvec_mean)
 
+        # If one diffusion sensitisation gradient direction is reversed with respect to
+        #   the other, still want to enable their recombination; but need to explicitly
+        #   account for this when averaging the two directions
+        if norm2 < 0.5:
+          bvec_mean = [ 0.5*(bvecs[0][pair[0]] - bvecs[0][pair[1]]),
+                        0.5*(bvecs[1][pair[0]] - bvecs[1][pair[1]]),
+                        0.5*(bvecs[2][pair[0]] - bvecs[2][pair[1]]) ]
+          norm2 = matrix.dot(bvec_mean, bvec_mean)
+
+        # Occasionally a b=0 volume can have a zero vector
+        if norm2:
+          factor = 1.0 / math.sqrt(norm2)
+          new_bvec = [ bvec_mean[0]*factor, bvec_mean[1]*factor, bvec_mean[2]*factor ]
+        else:
+          new_bvec = [ 0.0, 0.0, 0.0 ]
+
+        new_bval = 0.5 * (grad[pair[0]][3] + grad[pair[1]][3])
+      # End check for NaN values in bvec
+
+      bvecs_combined_transpose.append(new_bvec)
+      bvals_combined.append(new_bval)
+
+    if nans_present_bvecs:
+      app.warn('FSL "eddy"\'s rotated bvecs contained '
+               + (str(nans_present_bvecs) + ' entries with NaN values; '
+                      'these have been read as zero-length vectors, '
+                      'and will be saved as b=0 volumes'
+                  if nans_present_bvecs > 1
+                  else 'an entry with NaN values; '
+                      'this has been read as a zero-length vector, '
+                      'and will be saved as a b=0 volume')
+               + ' in the output series')
     bvecs_combined = matrix.transpose(bvecs_combined_transpose)
     matrix.save_matrix('bvecs_combined', bvecs_combined, add_to_command_history=False)
     matrix.save_vector('bvals_combined', bvals_combined, add_to_command_history=False)

core/dwi/dwi.h

@@ -1,4 +1,4 @@
-/* Copyright (c) 2008-2023 the MRtrix3 contributors.
+/* Copyright (c) 2008-2025 the MRtrix3 contributors.
  *
  * This Source Code Form is subject to the terms of the Mozilla Public
  * License, v. 2.0. If a copy of the MPL was not distributed with this

core/dwi/gradient.cpp

@@ -167,6 +167,38 @@ namespace MR
       grad.leftCols<3>().transpose() = header.transform().rotation() * G.transpose();
       grad.col(3) = bvals.row(0);
 
+      // Substitute NaNs with b=0 volumes
+      ssize_t nans_present_bvecs = false;
+      ssize_t nans_present_bvals = false;
+      ssize_t nan_linecount = 0;
+      for (ssize_t n = 0; n != grad.rows(); ++n) {
+        bool zero_row = false;
+        if (std::isnan(grad(n, 3))) {
+          if (grad.block<1,3>(n, 0).squaredNorm() > 0.0)
+            throw Exception("Corrupt content in bvecs/bvals data (" + bvecs_path + " & " + bvals_path + ") "
+                            "(NaN present in bval but valid direction in bvec)");
+          nans_present_bvals = true;
+          zero_row = true;
+        }
+        if (grad.block<1,3>(n, 0).hasNaN()) {
+          if (grad(n, 3) > 0.0)
+            throw Exception("Corrupt content in bvecs/bvals data (" + bvecs_path + " & " + bvals_path + ") "
+                            "(NaN bvec direction but non-zero value in bval)");
+          nans_present_bvecs = true;
+          zero_row = true;
+        }
+        if (zero_row) {
+          grad.block<1,4>(n, 0).setZero();
+          ++nan_linecount;
+        }
+      }
+      if (nan_linecount > 0) {
+        WARN(str(nan_linecount) + " row" + (nan_linecount > 1 ? "s" : "") + " with NaN values detected in "
+             + (nans_present_bvecs ? "bvecs file " + bvecs_path + (nans_present_bvals ? " and" : "") : "")
+             + (nans_present_bvals ? "bvals file " + bvals_path : "")
+             + "; these have been interpreted as b=0 volumes by MRtrix");
+      }
+
       return grad;
     }