Tolerance for dw_scheme differences on header merge

This prevents the diffusion gradient table from being discarded when merging the key-value contents of two image headers (eg. mrcalc, mrmath) due to inconsequential differences in their string representations from floating-point imprecision.
Closes #3015.

cmd/dwiextract.cpp

@@ -19,6 +19,7 @@
 #include "phase_encoding.h"
 #include "progressbar.h"
 #include "dwi/gradient.h"
+#include "dwi/shells.h"
 #include "algo/loop.h"
 #include "adapter/extract.h"
 

cmd/mrinfo.cpp

@@ -24,6 +24,7 @@
 #include "file/json.h"
 #include "file/json_utils.h"
 #include "dwi/gradient.h"
+#include "dwi/shells.h"
 #include "image_io/pipe.h"
 
 

cmd/mrtransform.cpp

@@ -588,7 +588,7 @@ void run ()
 
     if (interp == 0)
       output_header.datatype() = DataType::from_command_line (input_header.datatype());
-    auto output = Image<float>::create (argument[1], output_header).with_direct_io();
+    auto output = Image<float>::create (argument[1], output_header);
 
     switch (interp) {
       case 0:
@@ -630,7 +630,7 @@ void run ()
       add_line (output_header.keyval()["comments"], std::string ("resliced using warp image \"" + warp.name() + "\""));
     }
 
-    auto output = Image<float>::create(argument[1], output_header).with_direct_io();
+    auto output = Image<float>::create(argument[1], output_header);
 
     if (warp.ndim() == 5) {
       Image<default_type> warp_deform;
@@ -684,7 +684,7 @@ void run ()
       output_header.transform() = linear_transform;
     else
       output_header.transform() = linear_transform.inverse() * output_header.transform();
-    auto output = Image<float>::create (argument[1], output_header).with_direct_io();
+    auto output = Image<float>::create (argument[1], output_header);
     copy_with_progress (input, output);
 
     if (fod_reorientation || modulate_jac) {

core/dwi/gradient.cpp

@@ -15,8 +15,8 @@
  */
 
 #include "dwi/gradient.h"
+#include "file/config.h"
 #include "file/nifti_utils.h"
-#include "dwi/shells.h"
 
 namespace MR
 {
@@ -85,6 +85,17 @@ namespace MR
 
 
 
+//CONF option: BZeroThreshold
+//CONF default: 10.0
+//CONF Specifies the b-value threshold for determining those image
+//CONF volumes that correspond to b=0.
+    default_type bzero_threshold () {
+      static const default_type value = File::Config::get_float ("BZeroThreshold", DWI_BZERO_THREHSOLD_DEFAULT);
+      return value;
+    }
+
+
+
     BValueScalingBehaviour get_cmdline_bvalue_scaling_behaviour ()
     {
       auto opt = App::get_options ("bvalue_scaling");

core/dwi/gradient.h

@@ -19,14 +19,16 @@
 
 #include "app.h"
 #include "header.h"
-#include "dwi/shells.h"
+#include "types.h"
 #include "file/config.h"
 #include "file/path.h"
 #include "math/condition_number.h"
 #include "math/sphere.h"
 #include "math/SH.h"
 
 
+#define DWI_BZERO_THREHSOLD_DEFAULT 10.0
+
 namespace MR
 {
   namespace App { class OptionGroup; }
@@ -40,6 +42,7 @@ namespace MR
     extern App::Option bvalue_scaling_option;
     extern const char* const bvalue_scaling_description;
 
+    default_type bzero_threshold ();
 
 
     //! check that the DW scheme matches the DWI data in \a header
@@ -168,6 +171,64 @@ namespace MR
         }
       }
 
+    //! store the DW gradient encoding matrix in a KeyValues structure
+    /*! this will store the DW gradient encoding matrix under the key 'dw_scheme'.
+     */
+    template <class MatrixType>
+      void set_DW_scheme (KeyValues& keyval, const MatrixType& G)
+      {
+        if (!G.rows()) {
+          auto it = keyval.find ("dw_scheme");
+          if (it != keyval.end())
+            keyval.erase (it);
+          return;
+        }
+        keyval["dw_scheme"] = scheme2str (G);
+      }
+
+    template <class MatrixType1, class MatrixType2>
+      Eigen::MatrixXd resolve_DW_scheme (const MatrixType1& one, const MatrixType2& two)
+      {
+        if (one.rows() != two.rows())
+          throw Exception ("Unequal numbers of rows between gradient tables");
+        if (one.cols() != two.cols())
+          throw Exception ("Unequal numbers of rows between gradient tables");
+        Eigen::MatrixXd result (one.rows(), one.cols());
+        // Don't know how to intellegently combine data beyond four columns;
+        //   therefore don't proceed if such data are present and are not precisely equivalent
+        if (one.cols() > 4) {
+          if (one.rightCols(one.cols()-4) != two.rightCols(two.cols()-4))
+            throw Exception ("Unequal dw_scheme contents beyond standard four columns");
+          result.rightCols(one.cols()-4) = one.rightCols(one.cols()-4);
+        }
+        for (ssize_t rowindex = 0; rowindex != one.rows(); ++rowindex) {
+          auto one_dir = one.template block<1,3>(rowindex, 0);
+          auto two_dir = two.template block<1,3>(rowindex, 0);
+          const default_type one_bvalue = one(rowindex, 3);
+          const default_type two_bvalue = two(rowindex, 3);
+          const bool is_bzero = std::max(one_bvalue, two_bvalue) <= bzero_threshold();
+          if (one_dir == two_dir) {
+            result.block<1,3>(rowindex, 0) = one_dir;
+          } else {
+            const Eigen::Vector3d mean_dir = (one_dir + two_dir).normalized();
+            if (!is_bzero && mean_dir.dot (one_dir) < 1.0 - 1e-3) {
+              throw Exception("Diffusion vector directions not equal within permissible imprecision "
+                              "(row " + str(rowindex) + ": " + str(one_dir.transpose()) + " <--> " + str(two_dir.transpose())
+                              + "; dot product " + str(mean_dir.dot(one_dir)) + ")");
+            }
+            result.block<1,3>(rowindex, 0) = mean_dir;
+          }
+          if (one_bvalue == two_bvalue) {
+            result(rowindex, 3) = one_bvalue;
+          } else if (is_bzero || std::abs(one_bvalue - two_bvalue) <= 1.0) {
+            result(rowindex, 3) = 0.5 * (one_bvalue + two_bvalue);
+          } else {
+            throw Exception("Diffusion gradient table b-values not equivalent");
+          }
+        }
+        return result;
+      }
+
 
 
     //! parse the DW gradient encoding matrix from a header

core/dwi/shells.h

@@ -23,7 +23,7 @@
 
 #include "app.h"
 #include "types.h"
-#include "file/config.h"
+#include "dwi/gradient.h"
 #include "misc/bitset.h"
 
 
@@ -37,16 +37,8 @@
 // Default number of volumes necessary for a shell to be retained
 //   (note: only applies if function reject_small_shells() is called explicitly)
 #define DWI_SHELLS_MIN_DIRECTIONS 6
-// Default b-value threshold for a shell to be classified as "b=0"
-#define DWI_SHELLS_BZERO_THREHSOLD 10.0
 
 
-
-//CONF option: BZeroThreshold
-//CONF default: 10.0
-//CONF Specifies the b-value threshold for determining those image
-//CONF volumes that correspond to b=0.
-
 //CONF option: BValueEpsilon
 //CONF default: 80.0
 //CONF Specifies the difference between b-values necessary for image
@@ -64,10 +56,6 @@ namespace MR
 
     extern const App::OptionGroup ShellsOption;
 
-    FORCE_INLINE default_type bzero_threshold () {
-      static const default_type value = File::Config::get_float ("BZeroThreshold", DWI_SHELLS_BZERO_THREHSOLD);
-      return value;
-    }
 
 
 
@@ -87,7 +75,7 @@ namespace MR
         default_type get_min()   const { return min; }
         default_type get_max()   const { return max; }
 
-        bool is_bzero()   const { return (mean < bzero_threshold()); }
+        bool is_bzero() const { return (mean <= MR::DWI::bzero_threshold()); }
 
 
         bool operator< (const Shell& rhs) const { return (mean < rhs.mean); }

core/filter/dwi_brain_mask.h

@@ -29,6 +29,7 @@
 #include "algo/copy.h"
 #include "algo/loop.h"
 #include "dwi/gradient.h"
+#include "dwi/shells.h"
 
 
 namespace MR

core/header.cpp

@@ -21,6 +21,7 @@
 
 #include "app.h"
 #include "axes.h"
+#include "math/math.h"
 #include "mrtrix.h"
 #include "phase_encoding.h"
 #include "stride.h"
@@ -128,12 +129,21 @@ namespace MR
         }
       } else {
         auto it = keyval().find (item.first);
-        if (it == keyval().end() || it->second == item.second)
+        if (it == keyval().end() || it->second == item.second) {
           new_keyval.insert (item);
-        else if (item.first == "SliceTiming")
+        } else if (item.first == "SliceTiming") {
           new_keyval["SliceTiming"] = resolve_slice_timing (item.second, it->second);
-        else
+        } else if (item.first == "dw_scheme") {
+          try {
+            auto scheme = DWI::resolve_DW_scheme (parse_matrix (item.second), parse_matrix (it->second));
+            DWI::set_DW_scheme (new_keyval, scheme);
+          } catch (Exception& e) {
+            WARN("Error merging DW gradient tables between headers");
+            new_keyval["dw_scheme"] = "variable";
+          }
+        } else {
           new_keyval[item.first] = "variable";
+        }
       }
     }
     std::swap (keyval_, new_keyval);

src/dwi/sdeconv/csd.h

@@ -20,6 +20,7 @@
 #include "app.h"
 #include "header.h"
 #include "dwi/gradient.h"
+#include "dwi/shells.h"
 #include "math/SH.h"
 #include "math/ZSH.h"
 #include "dwi/directions/predefined.h"

testing/binaries/tests/mrcalc

@@ -2,3 +2,5 @@ mrcalc mrcalc/in.mif 2 -mult -neg -exp 10 -add - | testing_diff_image - mrcalc/o
 mrcalc mrcalc/in.mif 1.224 -div -cos mrcalc/in.mif -abs -sqrt -log -atanh -sub - | testing_diff_image - mrcalc/out2.mif -frac 1e-5
 mrcalc mrcalc/in.mif 0.2 -gt mrcalc/in.mif mrcalc/in.mif -1.123 -mult 0.9324 -add -exp -neg -if - | testing_diff_image - mrcalc/out3.mif -frac 1e-5
 mrcalc mrcalc/in.mif 0+1j -mult -exp mrcalc/in.mif -mult 1.34+5.12j -mult - | testing_diff_image - mrcalc/out4.mif -frac 1e-5
+mrinfo dwi.mif -dwgrad > tmp.txt && truncate -s-1 tmp.txt && mrconvert dwi.mif -grad tmp.txt - | mrcalc - dwi.mif -add 0.5 -mult - | mrinfo - -dwgrad
+if [ mrtransform dwi.mif -flip 0 - | mrcalc - dwi.mif -add 0.5 -mult - | mrinfo - -property dw_scheme ]; then exit 1; else exit 0; fi