DWI metadata handling fixes

cmd/amp2sh.cpp

@@ -21,7 +21,7 @@
 #include "file/matrix.h"
 #include "image.h"
 #include "math/SH.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 #include "progressbar.h"
 
 using namespace MR;
@@ -219,7 +219,7 @@ void run() {
       DWI::stash_DW_scheme(header, grad);
     }
   }
-  PhaseEncoding::clear_scheme(header);
+  Metadata::PhaseEncoding::clear_scheme(header.keyval());
 
   auto sh2amp = DWI::compute_SH2amp_mapping(dirs, true, 8);
 

cmd/dwi2adc.cpp

@@ -19,7 +19,7 @@
 #include "dwi/gradient.h"
 #include "image.h"
 #include "math/least_squares.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 #include "progressbar.h"
 
 using namespace MR;
@@ -44,7 +44,7 @@ void usage ()
     "adopts a 2-stage fitting strategy, in which the ADC is first estimated based on "
     "the DWI data with b > cutoff, and the other parameters are estimated subsequently. "
     "The output consists of 4 volumes, respectively S(0), D, f, and D'."
-    
+
   + "Note that this command ignores the gradient orientation entirely. This approach is "
     "therefore only suited for mean DWI (trace-weighted) images or for DWI data collected "
     "with isotropically-distributed gradient directions.";
@@ -60,12 +60,12 @@ void usage ()
     +   Argument ("bval").type_integer (0, 1000)
 
     + DWI::GradImportOptions();
-  
+
   REFERENCES
   + "Le Bihan, D.; Breton, E.; Lallemand, D.; Aubin, M.L.; Vignaud, J.; Laval-Jeantet, M. "
     "Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging. "
     "Radiology, 1988, 168, 497–505."
-    
+
   + "Jalnefjord, O.; Andersson, M.; Montelius; M.; Starck, G.; Elf, A.; Johanson, V.; Svensson, J.; Ljungberg, M. "
     "Comparison of methods for estimation of the intravoxel incoherent motion (IVIM) "
     "diffusion coefficient (D) and perfusion fraction (f). "
@@ -159,7 +159,7 @@ void run() {
   Header header(dwi);
   header.datatype() = DataType::Float32;
   DWI::stash_DW_scheme(header, grad);
-  PhaseEncoding::clear_scheme(header);
+  Metadata::PhaseEncoding::clear_scheme(header.keyval());
   header.ndim() = 4;
   header.size(3) = ivim ? 4 : 2;
 

cmd/dwi2fod.cpp

@@ -21,7 +21,7 @@
 #include "header.h"
 #include "image.h"
 #include "math/SH.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 
 #include "dwi/sdeconv/csd.h"
 #include "dwi/sdeconv/msmt_csd.h"
@@ -258,7 +258,7 @@ void run() {
     shared.init();
 
     DWI::stash_DW_scheme(header_out, shared.grad);
-    PhaseEncoding::clear_scheme(header_out);
+    Metadata::PhaseEncoding::clear_scheme(header_out.keyval());
 
     header_out.size(3) = shared.nSH();
     auto fod = Image<float>::create(argument[3], header_out);

cmd/dwi2tensor.cpp

@@ -22,7 +22,7 @@
 #include "file/matrix.h"
 #include "image.h"
 #include "math/constrained_least_squares.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 #include "progressbar.h"
 
 using namespace MR;
@@ -296,7 +296,7 @@ void run() {
   header.datatype() = DataType::Float32;
   header.ndim() = 4;
   DWI::stash_DW_scheme(header, grad);
-  PhaseEncoding::clear_scheme(header);
+  Metadata::PhaseEncoding::clear_scheme(header.keyval());
 
   Image<value_type> predict;
   opt = get_options("predicted_signal");

cmd/dwiextract.cpp

@@ -19,7 +19,7 @@
 #include "command.h"
 #include "dwi/gradient.h"
 #include "image.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 #include "progressbar.h"
 
 using namespace MR;
@@ -65,8 +65,8 @@ void usage() {
     + DWI::GradImportOptions()
     + DWI::ShellsOption
     + DWI::GradExportOptions()
-    + PhaseEncoding::ImportOptions
-    + PhaseEncoding::SelectOptions
+    + Metadata::PhaseEncoding::ImportOptions
+    + Metadata::PhaseEncoding::SelectOptions
     + Stride::Options;
 }
 // clang-format off
@@ -108,7 +108,7 @@ void run() {
   }
 
   auto opt = get_options("pe");
-  const auto pe_scheme = PhaseEncoding::get_scheme(input_image);
+  const auto pe_scheme = Metadata::PhaseEncoding::get_scheme(input_image);
   if (!opt.empty()) {
     if (!pe_scheme.rows())
       throw Exception("Cannot filter volumes by phase-encoding: No such information present");
@@ -154,7 +154,7 @@ void run() {
     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]);
-    PhaseEncoding::set_scheme(header, new_scheme);
+    Metadata::PhaseEncoding::set_scheme(header.keyval(), new_scheme);
   }
 
   auto output_image = Image<float>::create(argument[1], header);

cmd/mrcalc.cpp

@@ -320,7 +320,6 @@ UNARY_OP(
 #include "image.h"
 #include "math/rng.h"
 #include "memory.h"
-#include "phase_encoding.h"
 
 using namespace MR;
 using namespace App;
@@ -824,7 +823,7 @@ void get_header(const StackEntry &entry, Header &header) {
       header.spacing(n) = entry.image->spacing(n);
   }
 
-  header.merge_keyval(*entry.image);
+  header.merge_keyval(entry.image->keyval());
 }
 
 class ThreadFunctor {

cmd/mrcat.cpp

@@ -18,7 +18,6 @@
 #include "command.h"
 #include "dwi/gradient.h"
 #include "image.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
 
 using namespace MR;

cmd/mrconvert.cpp

@@ -24,7 +24,7 @@
 #include "file/ofstream.h"
 #include "header.h"
 #include "image.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 #include "transform.h"
 #include "types.h"
 
@@ -244,8 +244,8 @@ void usage() {
   + DWI::bvalue_scaling_option
   + DWI::GradExportOptions()
 
-  + PhaseEncoding::ImportOptions
-  + PhaseEncoding::ExportOptions;
+  + Metadata::PhaseEncoding::ImportOptions
+  + Metadata::PhaseEncoding::ExportOptions;
 
 }
 // clang-format on
@@ -271,7 +271,7 @@ void permute_DW_scheme(Header &H, const std::vector<int> &axes) {
 }
 
 void permute_PE_scheme(Header &H, const std::vector<int> &axes) {
-  auto in = PhaseEncoding::parse_scheme(H);
+  auto in = Metadata::PhaseEncoding::parse_scheme(H.keyval(), H);
   if (!in.rows())
     return;
 
@@ -285,16 +285,27 @@ void permute_PE_scheme(Header &H, const std::vector<int> &axes) {
   for (int row = 0; row != in.rows(); ++row)
     out.block<1, 3>(row, 0) = in.block<1, 3>(row, 0) * permute;
 
-  PhaseEncoding::set_scheme(H, out);
+  Metadata::PhaseEncoding::set_scheme(H.keyval(), out);
 }
 
 void permute_slice_direction(Header &H, const std::vector<int> &axes) {
-  auto it = H.keyval().find("SliceEncodingDirection");
-  if (it == H.keyval().end())
+  using Metadata::BIDS::axis_vector_type;
+  auto slice_encoding_it = H.keyval().find("SliceEncodingDirection");
+  auto slice_timing_it = H.keyval().find("SliceTiming");
+  if (slice_encoding_it == H.keyval().end() && slice_timing_it == H.keyval().end())
     return;
-  const Eigen::Vector3d orig_dir = Axes::id2dir(it->second);
-  const Eigen::Vector3d new_dir(orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]);
-  it->second = Axes::dir2id(new_dir);
+  if (slice_encoding_it == H.keyval().end()) {
+    const axis_vector_type orig_dir({0, 0, 1});
+    const axis_vector_type new_dir(orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]);
+    slice_encoding_it->second = Metadata::BIDS::vector2axisid(new_dir);
+    WARN("Header field \"SliceEncodingDirection\" inferred to be \"k\" in input image"
+         " and then transformed according to axis permutation"
+         " in order to preserve validity of existing header field \"SliceTiming\"");
+    return;
+  }
+  const axis_vector_type orig_dir = Metadata::BIDS::axisid2vector(slice_encoding_it->second);
+  const axis_vector_type new_dir(orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]);
+  slice_encoding_it->second = Metadata::BIDS::vector2axisid(new_dir);
 }
 
 template <class ImageType> inline std::vector<int> set_header(Header &header, const ImageType &input) {
@@ -352,7 +363,7 @@ void copy_permute(const InputType &in, Header &header_out, const std::string &ou
   const auto axes = set_header(header_out, in);
   auto out = Image<T>::create(output_filename, header_out, add_to_command_history);
   DWI::export_grad_commandline(out);
-  PhaseEncoding::export_commandline(out);
+  Metadata::PhaseEncoding::export_commandline(out);
   auto perm = Adapter::make<Adapter::PermuteAxes>(in, axes);
   threaded_copy_with_progress(perm, out, 0, std::numeric_limits<size_t>::max(), 2);
 }
@@ -381,20 +392,17 @@ void run() {
       throw;
     e.display(2);
   }
+  auto opt = get_options("json_import");
+  if (!opt.empty())
+    File::JSON::load(header_in, opt[0][0]);
+  if (!get_options("import_pe_table").empty() || !get_options("import_pe_eddy").empty())
+    Metadata::PhaseEncoding::set_scheme(header_in.keyval(), Metadata::PhaseEncoding::get_scheme(header_in));
 
   Header header_out(header_in);
   header_out.datatype() = DataType::from_command_line(header_out.datatype());
-
   if (header_in.datatype().is_complex() && !header_out.datatype().is_complex())
     WARN("requested datatype is real but input datatype is complex - imaginary component will be ignored");
 
-  if (!get_options("import_pe_table").empty() || !get_options("import_pe_eddy").empty())
-    PhaseEncoding::set_scheme(header_out, PhaseEncoding::get_scheme(header_in));
-
-  auto opt = get_options("json_import");
-  if (!opt.empty())
-    File::JSON::load(header_out, opt[0][0]);
-
   opt = get_options("copy_properties");
   if (!opt.empty()) {
     header_out.keyval().clear();
@@ -478,17 +486,17 @@ void run() {
         }
         Eigen::MatrixXd pe_scheme;
         try {
-          pe_scheme = PhaseEncoding::get_scheme(header_in);
+          pe_scheme = Metadata::PhaseEncoding::get_scheme(header_in);
           if (pe_scheme.rows()) {
             Eigen::MatrixXd extract_scheme(pos[3].size(), pe_scheme.cols());
             for (size_t vol = 0; vol != pos[3].size(); ++vol)
               extract_scheme.row(vol) = pe_scheme.row(pos[3][vol]);
-            PhaseEncoding::set_scheme(header_out, extract_scheme);
+            Metadata::PhaseEncoding::set_scheme(header_out.keyval(), extract_scheme);
           }
         } catch (...) {
-          WARN("Phase encoding scheme of input file does not match number of image volumes; omitting information from "
-               "output image");
-          PhaseEncoding::set_scheme(header_out, Eigen::MatrixXd());
+          WARN("Phase encoding scheme of input file does not match number of image volumes;"
+               " omitting information from output image");
+          Metadata::PhaseEncoding::clear_scheme(header_out.keyval());
         }
       }
     }

cmd/mrdegibbs.cpp

@@ -19,6 +19,7 @@
 #include "command.h"
 #include "degibbs/unring2d.h"
 #include "degibbs/unring3d.h"
+#include "metadata/bids.h"
 
 using namespace MR;
 using namespace App;
@@ -147,7 +148,8 @@ void run() {
       WARN("If data were acquired using multi-slice encoding, run in default 2D mode.");
     } else {
       try {
-        const Eigen::Vector3d slice_encoding_axis_onehot = Axes::id2dir(slice_encoding_it->second);
+        const Metadata::BIDS::axis_vector_type
+            slice_encoding_axis_onehot = Metadata::BIDS::axisid2vector(slice_encoding_it->second);
         std::vector<size_t> auto_slice_axes = {0, 0};
         if (slice_encoding_axis_onehot[0])
           auto_slice_axes = {1, 2};

cmd/mrinfo.cpp

@@ -23,7 +23,7 @@
 #include "file/json_utils.h"
 #include "header.h"
 #include "image_io/pipe.h"
-#include "phase_encoding.h"
+#include "metadata/phase_encoding.h"
 #include "types.h"
 
 using namespace MR;
@@ -119,7 +119,7 @@ void usage() {
       + Option ("shell_sizes", "list the number of volumes in each shell")
       + Option ("shell_indices", "list the image volumes attributed to each b-value shell")
 
-    + PhaseEncoding::ExportOptions
+    + Metadata::PhaseEncoding::ExportOptions
       + Option ("petable", "print the phase encoding table")
 
     + OptionGroup ("Handling of piped images")
@@ -252,7 +252,7 @@ void run() {
     ImageIO::Pipe::delete_piped_images = false;
 
   const bool export_grad = check_option_group(GradExportOptions);
-  const bool export_pe = check_option_group(PhaseEncoding::ExportOptions);
+  const bool export_pe = check_option_group(Metadata::PhaseEncoding::ExportOptions);
 
   if (export_grad && argument.size() > 1)
     throw Exception("can only export DW gradient table to file if a single input image is provided");
@@ -310,7 +310,7 @@ void run() {
     if (transform)
       print_transform(header);
     if (petable)
-      std::cout << PhaseEncoding::get_scheme(header) << "\n";
+      std::cout << Metadata::PhaseEncoding::get_scheme(header) << "\n";
 
     for (size_t n = 0; n < properties.size(); ++n)
       print_properties(header, properties[n][0]);
@@ -329,7 +329,7 @@ void run() {
     }
 
     DWI::export_grad_commandline(header);
-    PhaseEncoding::export_commandline(header);
+    Metadata::PhaseEncoding::export_commandline(header);
 
     if (json_keyval)
       File::JSON::write(header, *json_keyval, (argument.size() > 1 ? std::string("") : std::string(argument[0])));

cmd/mrmath.cpp

@@ -24,8 +24,8 @@
 #include "math/math.h"
 #include "math/median.h"
 #include "memory.h"
+#include "metadata/phase_encoding.h"
 #include "misc/voxel2vector.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
 
 #include <limits>
@@ -359,7 +359,7 @@ void run() {
       } catch (...) {
       }
       DWI::clear_DW_scheme(header_out);
-      PhaseEncoding::clear_scheme(header_out);
+      Metadata::PhaseEncoding::clear_scheme(header_out.keyval());
     }
 
     header_out.datatype() = DataType::from_command_line(DataType::Float32);
@@ -448,7 +448,7 @@ void run() {
           throw Exception("Image " + path + " has axis with non-unary dimension beyond first input image " +
                           header.name());
       }
-      header.merge_keyval(temp);
+      header.merge_keyval(temp.keyval());
     }
 
     // Instantiate a kernel depending on the operation requested

cmd/transformconvert.cpp

@@ -14,6 +14,7 @@
  * For more details, see http://www.mrtrix.org/.
  */
 
+#include "axes.h"
 #include "command.h"
 #include "file/key_value.h"
 #include "file/matrix.h"
@@ -64,15 +65,14 @@ void usage() {
 // clang-format on
 
 transform_type get_flirt_transform(const Header &header) {
-  std::vector<size_t> axes;
-  transform_type nifti_transform = File::NIfTI::adjust_transform(header, axes);
-  if (nifti_transform.matrix().topLeftCorner<3, 3>().determinant() < 0.0)
-    return nifti_transform;
+  const transform_type ondisk_transform = header.realignment().orig_transform();
+  if (ondisk_transform.matrix().topLeftCorner<3, 3>().determinant() < 0.0)
+    return ondisk_transform;
   transform_type coord_switch;
   coord_switch.setIdentity();
   coord_switch(0, 0) = -1.0f;
-  coord_switch(0, 3) = (header.size(axes[0]) - 1) * header.spacing(axes[0]);
-  return nifti_transform * coord_switch;
+  coord_switch(0, 3) = (header.size(header.realignment().permutation(0)) - 1) * header.spacing(header.realignment().permutation(0));
+  return ondisk_transform * coord_switch;
 }
 
 // transform_type parse_surfer_transform (const Header& header) {