DWI metadata fixes for master

cmd/amp2sh.cpp

@@ -16,12 +16,12 @@
 
 #include "command.h"
 #include "image.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
 #include "algo/threaded_loop.h"
 #include "dwi/gradient.h"
 #include "dwi/shells.h"
 #include "math/SH.h"
+#include "metadata/phase_encoding.h"
 
 
 using namespace MR;
@@ -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);
     }
   }
-  PhaseEncoding::clear_scheme (header);
+  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()) {

cmd/dwi2adc.cpp

@@ -16,11 +16,11 @@
 
 #include "command.h"
 #include "image.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
 #include "algo/threaded_copy.h"
-#include "math/least_squares.h"
 #include "dwi/gradient.h"
+#include "math/least_squares.h"
+#include "metadata/phase_encoding.h"
 
 
 using namespace MR;
@@ -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);
-  PhaseEncoding::clear_scheme (header);
-  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);

cmd/dwi2fod.cpp

@@ -17,11 +17,11 @@
 #include "command.h"
 #include "header.h"
 #include "image.h"
-#include "phase_encoding.h"
 #include "algo/threaded_loop.h"
 #include "dwi/gradient.h"
 #include "dwi/shells.h"
 #include "math/SH.h"
+#include "metadata/phase_encoding.h"
 
 #include "dwi/sdeconv/csd.h"
 #include "dwi/sdeconv/msmt_csd.h"
@@ -270,7 +270,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/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);

cmd/dwi2tensor.cpp

@@ -15,12 +15,12 @@
  */
 
 #include "command.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
 #include "image.h"
 #include "algo/threaded_copy.h"
 #include "dwi/gradient.h"
 #include "dwi/tensor.h"
+#include "metadata/phase_encoding.h"
 
 using namespace MR;
 using namespace App;
@@ -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);
-  PhaseEncoding::clear_scheme (header);
+  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);
 }
 

cmd/dwiextract.cpp

@@ -16,11 +16,11 @@
 
 #include "command.h"
 #include "image.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
-#include "dwi/gradient.h"
-#include "algo/loop.h"
 #include "adapter/extract.h"
+#include "algo/loop.h"
+#include "dwi/gradient.h"
+#include "metadata/phase_encoding.h"
 
 
 using namespace MR;
@@ -54,8 +54,8 @@ void usage ()
     + DWI::GradImportOptions()
     + DWI::ShellsOption
     + DWI::GradExportOptions()
-    + PhaseEncoding::ImportOptions
-    + PhaseEncoding::SelectOptions
+    + Metadata::PhaseEncoding::ImportOptions
+    + Metadata::PhaseEncoding::SelectOptions
     + Stride::Options;
 }
 
@@ -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 = 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]);
-    PhaseEncoding::set_scheme (header, 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);

cmd/mrcalc.cpp

@@ -182,7 +182,6 @@ UNARY_OP (atanh, "atanh (%1)", NORMAL, "inverse hyperbolic tangent", { return st
 #include "command.h"
 #include "image.h"
 #include "memory.h"
-#include "phase_encoding.h"
 #include "math/rng.h"
 #include "algo/threaded_copy.h"
 #include "dwi/gradient.h"
@@ -731,7 +730,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());
 }
 
 

cmd/mrcat.cpp

@@ -17,7 +17,6 @@
 #include "command.h"
 #include "image.h"
 #include "algo/loop.h"
-#include "phase_encoding.h"
 #include "progressbar.h"
 #include "dwi/gradient.h"