Fixes to bvecs import

In order for a bvecs file to be properly read in, the information about how the image was transformed upon image load needs to be available. However if an attempt is made to perform such a read based on a Header that it has itself been constructed from an arbitrary ImageType, this information is lost. Therefore, reading from a bvecs file needs to be done using the Header class in which the corresponding NIfTI was originally opened. In the process, some functionalities of the tracking of that transformation on image load have been made more explicit. Also fixes some syntax errors relating to the manual back-porting of #2917 in #3011.
MRtrix3 · Sep 22, 2024 · 178f1cc · 178f1cc
1 parent 01a09b8
commit 178f1cc
Show file tree

Hide file tree

Showing 23 changed files with 225 additions and 193 deletions.
diff --git a/cmd/amp2sh.cpp b/cmd/amp2sh.cpp
@@ -201,8 +201,9 @@ class Amp2SH { MEMALIGN(Amp2SH)
 
 void run ()
 {
-  auto amp = Image<value_type>::open (argument[0]).with_direct_io (3);
-  Header header (amp);
+  Header header_in (Header::open (argument[0]));
+  Header header_out (header_in);
+  header_out.datatype() = DataType::Float32;
 
   vector<size_t> bzeros, dwis;
   Eigen::MatrixXd dirs;
@@ -213,8 +214,8 @@ void run ()
       dirs = Math::Sphere::cartesian2spherical (dirs);
   }
   else {
-    auto hit = header.keyval().find ("directions");
-    if (hit != header.keyval().end()) {
+    auto hit = header_in.keyval().find ("directions");
+    if (hit != header_in.keyval().end()) {
       vector<default_type> dir_vector;
       for (auto line : split_lines (hit->second)) {
         auto v = parse_floats (line);
@@ -225,35 +226,34 @@ void run ()
         dirs(i/2, 0) = dir_vector[i];
         dirs(i/2, 1) = dir_vector[i+1];
       }
-      header.keyval()["basis_directions"] = hit->second;
-      header.keyval().erase (hit);
+      header_out.keyval()["basis_directions"] = hit->second;
+      header_out.keyval().erase (hit);
     }
     else {
-      auto grad = DWI::get_DW_scheme (amp);
+      auto grad = DWI::get_DW_scheme (header_in);
       DWI::Shells shells (grad);
       shells.select_shells (true, false, false);
       if (shells.smallest().is_bzero())
         bzeros = shells.smallest().get_volumes();
       dwis = shells.largest().get_volumes();
       dirs = DWI::gen_direction_matrix (grad, dwis);
-      DWI::stash_DW_scheme (header, grad);
+      DWI::stash_DW_scheme (header_out, grad);
     }
   }
-  Metadata::PhaseEncoding::clear_scheme (header.keyval());
+  Metadata::PhaseEncoding::clear_scheme (header_out.keyval());
 
   auto sh2amp = DWI::compute_SH2amp_mapping (dirs, true, 8);
 
-
   bool normalise = get_options ("normalise").size();
   if (normalise && !bzeros.size())
     throw Exception ("the normalise option is only available if the input data contains b=0 images.");
 
+  header_out.size (3) = sh2amp.cols();
+  Stride::set_from_command_line (header_out);
 
-  header.size (3) = sh2amp.cols();
-  Stride::set_from_command_line (header);
-  auto SH = Image<value_type>::create (argument[1], header);
-
+  auto amp = header_in.get_image<value_type>().with_direct_io (3);
   Amp2SHCommon common (sh2amp, bzeros, dwis, normalise);
+  auto SH = Image<value_type>::create (argument[1], header_out);
 
   opt = get_options ("rician");
   if (opt.size()) {

diff --git a/cmd/dwi2adc.cpp b/cmd/dwi2adc.cpp
@@ -80,11 +80,11 @@ class DWI2ADC { MEMALIGN(DWI2ADC)
 
 
 void run () {
-  auto dwi = Header::open (argument[0]).get_image<value_type>();
-  auto grad = DWI::get_DW_scheme (dwi);
+  auto header_in = Header::open (argument[0]);
+  auto grad = DWI::get_DW_scheme (header_in);
 
   size_t dwi_axis = 3;
-  while (dwi.size (dwi_axis) < 2)
+  while (header_in.size (dwi_axis) < 2)
     ++dwi_axis;
   INFO ("assuming DW images are stored along axis " + str (dwi_axis));
 
@@ -96,14 +96,15 @@ void run () {
 
   auto binv = Math::pinv (b);
 
-  Header header (dwi);
-  header.datatype() = DataType::Float32;
-  DWI::stash_DW_scheme (header, grad);
-  Metadata::PhaseEncoding::clear_scheme (header.keyval());
-  header.ndim() = 4;
-  header.size(3) = 2;
+  Header header_out (header_in);
+  header_out.datatype() = DataType::Float32;
+  DWI::stash_DW_scheme (header_out, grad);
+  Metadata::PhaseEncoding::clear_scheme (header_out.keyval());
+  header_out.ndim() = 4;
+  header_out.size(3) = 2;
 
-  auto adc = Image<value_type>::create (argument[1], header);
+  auto dwi = header_in.get_image<value_type>();
+  auto adc = Image<value_type>::create (argument[1], header_out);
 
   ThreadedLoop ("computing ADC values", dwi, 0, 3)
     .run (DWI2ADC (binv, dwi_axis), dwi, adc);

diff --git a/cmd/dwi2mask.cpp b/cmd/dwi2mask.cpp
@@ -66,13 +66,13 @@ OPTIONS
 
 void run () {
 
-  auto input = Image<float>::open (argument[0]).with_direct_io (3);
-  auto grad = DWI::get_DW_scheme (input);
-
-  if (input.ndim() != 4)
+  auto header = Header::open(argument[0]);
+  if (header.ndim() != 4)
     throw Exception ("input DWI image must be 4D");
+  auto grad = DWI::get_DW_scheme (header);
+  auto input = header.get_image<float>().with_direct_io (3);
 
-  Filter::DWIBrainMask dwi_brain_mask_filter (input, grad);
+  Filter::DWIBrainMask dwi_brain_mask_filter (header, grad);
   dwi_brain_mask_filter.set_message ("computing dwi brain mask");
   auto temp_mask = Image<bool>::scratch (dwi_brain_mask_filter, "brain mask");
   dwi_brain_mask_filter (input, temp_mask);

diff --git a/cmd/dwi2tensor.cpp b/cmd/dwi2tensor.cpp
@@ -202,52 +202,53 @@ inline Processor<MASKType, B0Type, DKTType, PredictType> processor (const Eigen:
 
 void run ()
 {
-  auto dwi = Header::open (argument[0]).get_image<value_type>();
-  auto grad = DWI::get_DW_scheme (dwi);
+  auto header_in = Header::open (argument[0]);
+  auto grad = DWI::get_DW_scheme (header_in);
 
   Image<bool> mask;
   auto opt = get_options ("mask");
   if (opt.size()) {
     mask = Image<bool>::open (opt[0][0]);
-    check_dimensions (dwi, mask, 0, 3);
+    check_dimensions (header_in, mask, 0, 3);
   }
 
   bool ols = get_options ("ols").size();
 
   // depending on whether first (initialisation) loop should be considered an iteration
   auto iter = get_option_value ("iter", DEFAULT_NITER);
 
-  Header header (dwi);
-  header.datatype() = DataType::Float32;
-  header.ndim() = 4;
-  DWI::stash_DW_scheme (header, grad);
-  Metadata::PhaseEncoding::clear_scheme (header.keyval());
+  Header header_out (header_in);
+  header_out.datatype() = DataType::Float32;
+  header_out.ndim() = 4;
+  DWI::stash_DW_scheme (header_out, grad);
+  Metadata::PhaseEncoding::clear_scheme (header_out.keyval());
 
   Image<value_type> predict;
   opt = get_options ("predicted_signal");
   if (opt.size())
-    predict = Image<value_type>::create (opt[0][0], header);
+    predict = Image<value_type>::create (opt[0][0], header_out);
 
-  header.size(3) = 6;
-  auto dt = Image<value_type>::create (argument[1], header);
+  header_out.size(3) = 6;
+  auto dt = Image<value_type>::create (argument[1], header_out);
 
   Image<value_type> b0;
   opt = get_options ("b0");
   if (opt.size()) {
-    header.ndim() = 3;
-    b0 = Image<value_type>::create (opt[0][0], header);
+    header_out.ndim() = 3;
+    b0 = Image<value_type>::create (opt[0][0], header_out);
   }
 
   Image<value_type> dkt;
   opt = get_options ("dkt");
   if (opt.size()) {
-    header.ndim() = 4;
-    header.size(3) = 15;
-    dkt = Image<value_type>::create (opt[0][0], header);
+    header_out.ndim() = 4;
+    header_out.size(3) = 15;
+    dkt = Image<value_type>::create (opt[0][0], header_out);
   }
 
-  Eigen::MatrixXd b = -DWI::grad2bmatrix<double> (grad, opt.size()>0);
+  Eigen::MatrixXd b = -DWI::grad2bmatrix<double> (grad, dkt.valid());
 
+  auto dwi = Header::open (argument[0]).get_image<value_type>();
   ThreadedLoop("computing tensors", dwi, 0, 3).run (processor (b, ols, iter, mask, b0, dkt, predict), dwi, dt);
 }
 
diff --git a/cmd/dwiextract.cpp b/cmd/dwiextract.cpp
@@ -66,10 +66,10 @@ void usage ()
 
 void run()
 {
-  auto input_image = Image<float>::open (argument[0]);
-  if (input_image.ndim() < 4)
+  auto header_in = Header::open (argument[0]);
+  if (header_in.ndim() < 4)
     throw Exception ("Epected input image to contain more than three dimensions");
-  auto grad = DWI::get_DW_scheme (input_image);
+  auto grad = DWI::get_DW_scheme (header_in);
 
   // Want to support non-shell-like data if it's just a straight extraction
   //   of all dwis or all bzeros i.e. don't initialise the Shells class
@@ -101,7 +101,7 @@ void run()
   }
 
   auto opt = get_options ("pe");
-  const auto pe_scheme = Metadata::PhaseEncoding::get_scheme (input_image);
+  const auto pe_scheme = Metadata::PhaseEncoding::get_scheme (header_in);
   if (opt.size()) {
     if (!pe_scheme.rows())
       throw Exception ("Cannot filter volumes by phase-encoding: No such information present");
@@ -134,24 +134,25 @@ void run()
 
   std::sort (volumes.begin(), volumes.end());
 
-  Header header (input_image);
-  Stride::set_from_command_line (header);
-  header.size (3) = volumes.size();
+  Header header_out (header_in);
+  Stride::set_from_command_line (header_out);
+  header_out.size (3) = volumes.size();
 
   Eigen::MatrixXd new_grad (volumes.size(), grad.cols());
   for (size_t i = 0; i < volumes.size(); i++)
     new_grad.row (i) = grad.row (volumes[i]);
-  DWI::set_DW_scheme (header, new_grad);
+  DWI::set_DW_scheme (header_out, new_grad);
 
   if (pe_scheme.rows()) {
     Eigen::MatrixXd new_scheme (volumes.size(), pe_scheme.cols());
     for (size_t i = 0; i != volumes.size(); ++i)
       new_scheme.row(i) = pe_scheme.row (volumes[i]);
-    Metadata::PhaseEncoding::set_scheme (header.keyval(), new_scheme);
+    Metadata::PhaseEncoding::set_scheme (header_out.keyval(), new_scheme);
   }
 
-  auto output_image = Image<float>::create (argument[1], header);
-  DWI::export_grad_commandline (header);
+  auto input_image = Image<float>::open (argument[0]);
+  auto output_image = Image<float>::create (argument[1], header_out);
+  DWI::export_grad_commandline (header_out);
 
   auto input_volumes = Adapter::make<Adapter::Extract1D> (input_image, 3, volumes);
   threaded_copy_with_progress_message ("extracting volumes", input_volumes, output_image);

diff --git a/cmd/tckglobal.cpp b/cmd/tckglobal.cpp
@@ -218,8 +218,7 @@ void run ()
   using namespace DWI::Tractography::GT;
 
   // Inputs -----------------------------------------------------------------------------
-
-  auto dwi = Image<float>::open(argument[0]).with_direct_io(3);
+  Header header_in = Header::open (argument[0]);
 
   Properties properties;
   properties.resp_WM = load_matrix<float> (argument[1]);
@@ -237,7 +236,7 @@ void run ()
   opt = get_options("mask");
   if (opt.size()) {
     mask = Image<bool>::open(opt[0][0]);
-    check_dimensions(dwi, mask, 0, 3);
+    check_dimensions(header_in, mask, 0, 3);
   }
 
   // Parameters -------------------------------------------------------------------------
@@ -304,9 +303,10 @@ void run ()
   if (opt.size())
     stats.open_stream(opt[0][0]);
 
+  auto dwi = header_in.get_image<float>().with_direct_io(3);
   ParticleGrid pgrid (dwi);
 
-  ExternalEnergyComputer* Eext = new ExternalEnergyComputer(stats, dwi, properties);
+  ExternalEnergyComputer* Eext = new ExternalEnergyComputer(stats, header_in, properties);
   InternalEnergyComputer* Eint = new InternalEnergyComputer(stats, pgrid);
   Eint->setConnPot(cpot);
   EnergySumComputer* Esum = new EnergySumComputer(stats, Eint, properties.lam_int, Eext, properties.lam_ext / ( wmscale2 * properties.weight*properties.weight));
@@ -345,14 +345,14 @@ void run ()
 
 
   // Save fiso, tod and eext
-  Header header (dwi);
-  header.datatype() = DataType::Float32;
+  Header header_out (dwi);
+  header_out.datatype() = DataType::Float32;
 
   opt = get_options("fod");
   if (opt.size()) {
     INFO("Saving fODF image to file");
-    header.size(3) = Math::SH::NforL(properties.Lmax);
-    auto fODF = Image<float>::create (opt[0][0], header);
+    header_out.size(3) = Math::SH::NforL(properties.Lmax);
+    auto fODF = Image<float>::create (opt[0][0], header_out);
     auto f = __copy_fod<float>(properties.Lmax, properties.weight, !get_options("noapo").size());
     ThreadedLoop(Eext->getTOD(), 0, 3).run(f, Eext->getTOD(), fODF);
   }
@@ -361,8 +361,8 @@ void run ()
   if (opt.size()) {
     if (properties.resp_ISO.size() > 0) {
       INFO("Saving isotropic fractions to file");
-      header.size(3) = properties.resp_ISO.size();
-      auto Fiso = Image<float>::create (opt[0][0], header);
+      header_out.size(3) = properties.resp_ISO.size();
+      auto Fiso = Image<float>::create (opt[0][0], header_out);
       threaded_copy(Eext->getFiso(), Fiso);
     }
     else {
@@ -373,8 +373,8 @@ void run ()
   opt = get_options("eext");
   if (opt.size()) {
     INFO("Saving external energy to file");
-    header.ndim() = 3;
-    auto EextI = Image<float>::create (opt[0][0], header);
+    header_out.ndim() = 3;
+    auto EextI = Image<float>::create (opt[0][0], header_out);
     threaded_copy(Eext->getEext(), EextI);
   }
 

diff --git a/core/axes.h b/core/axes.h
@@ -19,6 +19,7 @@
 
 
 #include <array>
+#include <set>
 
 #include "types.h"
 
@@ -31,15 +32,21 @@ namespace MR
 
 
 
-    using permutations_type = std::array<size_t, 3>;
+    using permutations_type = std::array<ssize_t, 3>;
     using flips_type = std::array<bool, 3>;
-    class Shuffle {
+    class Shuffle { NOMEMALIGN
     public:
-      Shuffle() : permutations ({0, 1, 2}), flips ({false, false, false}) {}
-      operator bool() const {
+      Shuffle() : permutations ({-1, -1, -1}), flips ({false, false, false}) {}
+      bool is_identity() const {
         return (permutations[0] != 0 || permutations[1] != 1 || permutations[2] != 2 || //
                 flips[0] || flips[1] || flips[2]);
       }
+      bool is_set() const {
+        return permutations != permutations_type{-1, -1, -1};
+      }
+      bool valid() const {
+        return std::set<ssize_t>(permutations.begin(), permutations.end()) == std::set<ssize_t>({0, 1, 2});
+      }
       permutations_type permutations;
       flips_type flips;
     };

diff --git a/core/dwi/gradient.cpp b/core/dwi/gradient.cpp
@@ -115,6 +115,8 @@ namespace MR
 
     Eigen::MatrixXd load_bvecs_bvals (const Header& header, const std::string& bvecs_path, const std::string& bvals_path)
     {
+      assert (header.realignment().orig_transform().matrix().allFinite());
+
       Eigen::MatrixXd bvals, bvecs;
       try {
         bvals = load_matrix<> (bvals_path);
@@ -167,8 +169,8 @@ namespace MR
       // transform have negative determinant:
       if (adjusted_transform.linear().determinant() > 0.0)
         bvecs.row(0) = -bvecs.row(0);
-      save_matrix(bvecs, bvecs_path, KeyValues(), false);
-      save_matrix(grad.col(3), bvals_path, KeyValues(), false);
+      save_matrix (bvecs, bvecs_path, KeyValues(), false);
+      save_matrix (grad.col(3).transpose(), bvals_path, KeyValues(), false);
     }
 
 
@@ -257,6 +259,7 @@ namespace MR
                 "(maximum scaling factor = " + str(max_scaling_factor) + ")");
           }
         }
+        assert (grad.allFinite());
 
         // write the scheme as interpreted back into the header if:
         // - vector normalisation effect is large, regardless of whether or not b-value scaling was applied