diff --git a/cmd/amp2sh.cpp b/cmd/amp2sh.cpp
index e28d4995c1..cd137b7b19 100644
--- 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
index 49812e0950..f49d41b169 100644
--- 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
index 466f62c5ad..c95451cc9d 100644
--- 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
index cca2d52b87..8d77c0bfbd 100644
--- a/cmd/dwi2tensor.cpp
+++ b/cmd/dwi2tensor.cpp
@@ -202,14 +202,14 @@ 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();
@@ -217,37 +217,38 @@ void run ()
   // 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
index d60d947f1e..791eebd3b7 100644
--- 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
index 7c4f4e560e..03a7352809 100644
--- 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
index f7d2b88cfd..12de59cbcd 100644
--- 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
index 8d93dc2704..95e55fffee 100644
--- 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
diff --git a/core/file/nifti_utils.cpp b/core/file/nifti_utils.cpp
index 9ddf624d06..ae71d7c272 100644
--- a/core/file/nifti_utils.cpp
+++ b/core/file/nifti_utils.cpp
@@ -588,7 +588,7 @@ namespace MR
         strides.resize (3);
         auto order = Stride::order (strides);
         Axes::Shuffle result;
-        result.permutations = {order[0], order[1], order[2]};
+        result.permutations = {ssize_t(order[0]), ssize_t(order[1]), ssize_t(order[2])};
         result.flips = {strides[order[0]] < 0, strides[order[1]] < 0, strides[order[2]] < 0};
         return result;
       }
@@ -600,7 +600,7 @@ namespace MR
       {
         const Axes::Shuffle shuffle = axes_on_write(H);
         axes = shuffle.permutations;
-        if (!shuffle)
+        if (shuffle.is_identity())
           return H.transform();
 
         const auto& M_in = H.transform().matrix();
diff --git a/core/filter/dwi_brain_mask.h b/core/filter/dwi_brain_mask.h
index a6b2684b97..0b2d9f739e 100644
--- a/core/filter/dwi_brain_mask.h
+++ b/core/filter/dwi_brain_mask.h
@@ -46,8 +46,9 @@ namespace MR
      *
      * Typical usage:
      * \code
+     * auto header = Header::open (argument[0]);
+     * auto grad = DWI::get_DW_scheme (header);
      * auto input = Image<value_type>::open (argument[0]);
-     * auto grad = DWI::get_DW_scheme (input);
      * Filter::DWIBrainMask dwi_brain_mask_filter (input, grad);
      * auto output = Image<bool>::create (argument[1], dwi_brain_mask_filter);
      * dwi_brain_mask_filter (input, output);
diff --git a/core/header.cpp b/core/header.cpp
index 6151a154fe..24408db1cb 100644
--- a/core/header.cpp
+++ b/core/header.cpp
@@ -632,7 +632,7 @@ namespace MR
       return;
 
     realignment_.shuffle_ = Axes::get_shuffle_to_make_RAS(transform());
-    if (!realignment_)
+    if (realignment_.is_identity())
       return;
 
     auto M (transform());
@@ -811,9 +811,9 @@ namespace MR
 
 
 
-  Header::Realignment::Realignment() :
-      applied_transform_(applied_transform_type::Identity()) {
-    orig_transform_.matrix().fill(std::numeric_limits<default_type>::quiet_NaN());
+  Header::Realignment::Realignment() {
+    applied_transform_.matrix().fill(0);
+    orig_transform_.matrix().fill(std::numeric_limits<default_type>::signaling_NaN());
   }
 
 
diff --git a/core/header.h b/core/header.h
index 6abaa98435..0822dba4b3 100644
--- a/core/header.h
+++ b/core/header.h
@@ -111,9 +111,11 @@ namespace MR
           Header (static_cast<const Header&> (original)) { }
 
       //! copy constructor from type of class other than Header
-      /*! This copies all relevant parameters over from \a original. */
+      /*! This copies all relevant parameters over from \a original.
+       * Note that information about transform realignment on image load will not be available.
+       */
       template <class HeaderType, typename std::enable_if<!std::is_base_of<Header, HeaderType>::value, void*>::type = nullptr>
-        Header (const HeaderType& original) :
+        explicit Header (const HeaderType& original) :
           transform_ (original.transform()),
           name_ (original.name()),
           keyval_ (original.keyval()),
@@ -203,18 +205,20 @@ namespace MR
 
       // Class to store all information relating to internal transform realignment
       class Realignment {
+        MEMALIGN(Realignment)
         public:
           // From one image space to another image space;
           //   linear component is permutations & flips only,
           //   transformation is in voxel count,
           //   therefore can store as integer
+          // TODO Calculate translations; turn into affine transform; verify
           using applied_transform_type = Eigen::Matrix<int, 3, 3>;
           Realignment();
-          Realignment (Header&);
-          operator bool() const { return bool(shuffle_); }
-          const std::array<size_t, 3>& permutations() const { return shuffle_.permutations; }
+          bool is_identity() const { return shuffle_.is_identity(); }
+          bool valid() const { return shuffle_.valid(); }
+          const Axes::permutations_type& permutations() const { return shuffle_.permutations; }
           size_t permutation (const size_t axis) const { assert(axis < 3); return shuffle_.permutations[axis]; }
-          const std::array<bool, 3>& flips() const { return shuffle_.flips; }
+          const Axes::flips_type& flips() const { return shuffle_.flips; }
           bool flip (const size_t axis) const { assert(axis < 3); return shuffle_.flips[axis]; }
           const transform_type& orig_transform() const { return orig_transform_; }
           const applied_transform_type& applied_transform() const { return applied_transform_; }
diff --git a/core/image.h b/core/image.h
index 3dbb3aab6a..d628c0d786 100644
--- a/core/image.h
+++ b/core/image.h
@@ -61,6 +61,7 @@ namespace MR
 
         FORCE_INLINE const std::string& name() const { return buffer->name(); }
         FORCE_INLINE const transform_type& transform() const { return buffer->transform(); }
+        //FORCE_INLINE const Header::Realignment& realignment() const { return buffer->realignment(); }
 
         FORCE_INLINE size_t  ndim () const { return buffer->ndim(); }
         FORCE_INLINE ssize_t size (size_t axis) const { return buffer->size (axis); }
diff --git a/core/metadata/phase_encoding.cpp b/core/metadata/phase_encoding.cpp
index 4794502ba0..d130f8cb7f 100644
--- a/core/metadata/phase_encoding.cpp
+++ b/core/metadata/phase_encoding.cpp
@@ -212,7 +212,7 @@ namespace MR {
           DEBUG("No phase encoding information found for transformation with load of image \"" + H.name() + "\"");
           return;
         }
-        if (!H.realignment()) {
+        if (H.realignment().is_identity()) {
           INFO("No transformation of phase encoding data for load of image \"" + H.name() + "\" required");
           return;
         }
@@ -223,7 +223,7 @@ namespace MR {
 
 
       scheme_type transform_for_image_load(const scheme_type& pe_scheme, const Header& H) {
-        if (!H.realignment())
+        if (H.realignment().is_identity())
           return pe_scheme;
         scheme_type result(pe_scheme.rows(), pe_scheme.cols());
         for (ssize_t row = 0; row != pe_scheme.rows(); ++row) {
@@ -251,7 +251,7 @@ namespace MR {
         if (pe_scheme.rows() == 0)
           return pe_scheme;
         Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H);
-        if (!shuffle) {
+        if (shuffle.is_identity()) {
           INFO("No transformation of phase encoding data required for export to file:"
                " output image will be RAS");
           return pe_scheme;
@@ -358,6 +358,19 @@ namespace MR {
 
 
 
+      void save(const scheme_type& PE, const std::string& path) {
+        File::OFStream out(path);
+        for (ssize_t row = 0; row != PE.rows(); ++row) {
+          // Write phase-encode direction as integers; other information as floating-point
+          out << PE.template block<1, 3>(row, 0).template cast<int>();
+          if (PE.cols() > 3)
+            out << " " << PE.block(row, 3, 1, PE.cols() - 3);
+          out << "\n";
+        }
+      }
+
+
+
     }
   }
 }
diff --git a/core/metadata/phase_encoding.h b/core/metadata/phase_encoding.h
index 296984002e..e10f30ff4d 100644
--- a/core/metadata/phase_encoding.h
+++ b/core/metadata/phase_encoding.h
@@ -93,20 +93,17 @@ namespace MR {
       void transform_for_nifti_write(KeyValues& keyval, const Header& H);
       scheme_type transform_for_nifti_write(const scheme_type& pe_scheme, const Header& H);
 
-      namespace {
-        void _save(const scheme_type& PE, const std::string& path) {
-          File::OFStream out(path);
-          for (ssize_t row = 0; row != PE.rows(); ++row) {
-            // Write phase-encode direction as integers; other information as floating-point
-            out << PE.template block<1, 3>(row, 0).template cast<int>();
-            if (PE.cols() > 3)
-              out << " " << PE.block(row, 3, 1, PE.cols() - 3);
-            out << "\n";
-          }
-        }
-      } // namespace
+      void save(const scheme_type& PE, const std::string& path);
+
+      template <class HeaderType>
+      void save(const HeaderType &header, const std::string &path) {
+        const scheme_type scheme = get_scheme(header);
+        if (scheme.rows() == 0)
+          throw Exception ("No phase encoding scheme in header of image \"" + header.name() + "\" to save");
+        save(scheme, header, path);
+      }
 
-      //! Save a phase-encoding scheme to file
+      //! Save a phase-encoding scheme associated with an image to file
       // Note that because the output table requires permutation / sign flipping
       //   only if the output target image is a NIfTI, the output file name must have
       //   already been set
@@ -119,9 +116,9 @@ namespace MR {
         }
 
         if (Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) {
-          _save(transform_for_nifti_write(PE, header), path);
+          save(transform_for_nifti_write(PE, header), path);
         } else {
-          _save(PE, path);
+          save(PE, path);
         }
       }
 
diff --git a/core/metadata/slice_encoding.cpp b/core/metadata/slice_encoding.cpp
index e9109826d9..04235540ec 100644
--- a/core/metadata/slice_encoding.cpp
+++ b/core/metadata/slice_encoding.cpp
@@ -33,7 +33,7 @@ namespace MR {
         auto slice_encoding_it = keyval.find("SliceEncodingDirection");
         auto slice_timing_it = keyval.find("SliceTiming");
         if (!(slice_encoding_it == keyval.end() && slice_timing_it == keyval.end())) {
-          if (!header.realignment()) {
+          if (header.realignment().is_identity()) {
             INFO("No transformation of slice encoding direction for load of image \"" + header.name() + "\" required");
             return;
           }
@@ -71,7 +71,7 @@ namespace MR {
           return;
 
         const Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H);
-        if (!shuffle) {
+        if (shuffle.is_identity()) {
           INFO("No need to transform slice encoding information for NIfTI image write:"
                " image is already RAS");
           return;
@@ -106,8 +106,8 @@ namespace MR {
       std::string resolve_slice_timing (const std::string &one, const std::string &two) {
         if (one == "variable" || two == "variable")
           return "variable";
-        std::vector<std::string> one_split = split(one, ",");
-        std::vector<std::string> two_split = split(two, ",");
+        const vector<std::string> one_split = split(one, ",");
+        const vector<std::string> two_split = split(two, ",");
         if (one_split.size() != two_split.size()) {
           DEBUG("Slice timing vectors of inequal length");
           return "invalid";
diff --git a/src/dwi/tractography/GT/externalenergy.cpp b/src/dwi/tractography/GT/externalenergy.cpp
index 6014d56a27..012f92055a 100644
--- a/src/dwi/tractography/GT/externalenergy.cpp
+++ b/src/dwi/tractography/GT/externalenergy.cpp
@@ -28,17 +28,17 @@ namespace MR {
     namespace Tractography {
       namespace GT {
 
-        ExternalEnergyComputer::ExternalEnergyComputer(Stats& stat, const Image<float>& dwimage, const Properties& props)
+        ExternalEnergyComputer::ExternalEnergyComputer(Stats& stat, Header& dwiheader, const Properties& props)
           : EnergyComputer(stat),
-            dwi(dwimage),
-            T(Transform(dwimage).scanner2voxel),
+            dwi(dwiheader.get_image<float>().with_direct_io(3)),
+            T(Transform(dwiheader).scanner2voxel),
             lmax(props.Lmax), ncols(Math::SH::NforL(lmax)), nf(props.resp_ISO.size()),
             beta(props.beta), mu(props.ppot*M_sqrt4PI), dE(0.0)
         {
           DEBUG("Initialise computation of external energy.");
 
           // Create images --------------------------------------------------------------
-          Header header (dwimage);
+          Header header (dwiheader);
           header.datatype() = DataType::Float32;
           header.size(3) = ncols;
           tod = Image<float>::scratch(header, "TOD image");
@@ -54,7 +54,7 @@ namespace MR {
           eext = Image<float>::scratch(header, "external energy");
 
           // Set kernel matrices --------------------------------------------------------
-          auto grad = DWI::get_DW_scheme(dwimage);
+          auto grad = DWI::get_DW_scheme(dwiheader);
           nrows = grad.rows();
           DWI::Shells shells (grad);
 
diff --git a/src/dwi/tractography/GT/externalenergy.h b/src/dwi/tractography/GT/externalenergy.h
index 332c005d3a..9267c4a394 100644
--- a/src/dwi/tractography/GT/externalenergy.h
+++ b/src/dwi/tractography/GT/externalenergy.h
@@ -30,83 +30,83 @@ namespace MR {
   namespace DWI {
     namespace Tractography {
       namespace GT {
-        
+
         class ExternalEnergyComputer : public EnergyComputer
         { MEMALIGN(ExternalEnergyComputer)
         public:
-          
-          ExternalEnergyComputer(Stats& stat, const Image<float>& dwimage, const Properties& props);
-          
-          
+
+          ExternalEnergyComputer(Stats& stat, Header& dwiheader, const Properties& props);
+
+
           Image<float>& getTOD() { return tod; }
           Image<float>& getFiso() { return fiso; }
           Image<float>& getEext() { return eext; }
-          
+
           void resetEnergy();
-          
+
           double stageAdd(const Point_t& pos, const Point_t& dir)
           {
             add(pos, dir, 1.0);
             return eval();
           }
-          
+
           double stageShift(const Particle* par, const Point_t& pos, const Point_t& dir)
           {
             add(par->getPosition(), par->getDirection(), -1.0);
             add(pos, dir, 1.0);
             return eval();
           }
-          
+
           double stageRemove(const Particle* par)
           {
             add(par->getPosition(), par->getDirection(), -1.0);
             return eval();
           }
-          
+
           void acceptChanges();
-          
+
           void clearChanges();
-          
+
           EnergyComputer* clone() const { return new ExternalEnergyComputer(*this); }
-          
-          
-          
+
+
+
         protected:
-          
+
           Image<float> dwi;
           Image<float> tod;
           Image<float> fiso;
           Image<float> eext;
-          
+
           transform_type T;
-          
-          int lmax; 
+
+          int lmax;
           size_t nrows, ncols, nf;
           double beta, mu, dE;
           Eigen::MatrixXd K, Ak;
           Eigen::VectorXd y, t, d, fk;
-          
+
           Math::ICLS::Problem<double> nnls;
-          
+
           vector<Eigen::Vector3i > changes_vox;
           vector<Eigen::VectorXd > changes_tod;
           vector<Eigen::VectorXd > changes_fiso;
           vector<double> changes_eext;
-          
-          
+
+
           void add(const Point_t& pos, const Point_t& dir, const double factor = 1.);
-          
+
           void add2vox(const Eigen::Vector3i& vox, const double w);
-          
+
           double eval();
-          
+
           double calcEnergy();
-          
+
           inline double hanning(const double w) const
           {
             return (w <= (1.0-beta)/2) ? 0.0 : (w >= (1.0+beta)/2) ? 1.0 : (1 - std::cos(Math::pi * (w-(1.0-beta)/2)/beta )) / 2;
           }
-          
+
         };
 
 
diff --git a/src/dwi/tractography/GT/particlegrid.cpp b/src/dwi/tractography/GT/particlegrid.cpp
index fb2b8018eb..67019098e7 100644
--- a/src/dwi/tractography/GT/particlegrid.cpp
+++ b/src/dwi/tractography/GT/particlegrid.cpp
@@ -20,15 +20,33 @@ namespace MR {
   namespace DWI {
     namespace Tractography {
       namespace GT {
-        
-        
+
+
+        ParticleGrid::ParticleGrid(const Header& H)
+        {
+          DEBUG("Initialise particle grid.");
+          dims[0] = Math::ceil<size_t>( H.size(0) * H.spacing(0) / (2.0*Particle::L) );
+          dims[1] = Math::ceil<size_t>( H.size(1) * H.spacing(1) / (2.0*Particle::L) );
+          dims[2] = Math::ceil<size_t>( H.size(2) * H.spacing(2) / (2.0*Particle::L) );
+          grid.resize(dims[0]*dims[1]*dims[2]);
+
+          // Initialise scanner-to-grid transform
+          Eigen::DiagonalMatrix<default_type, 3> newspacing (2.0*Particle::L, 2.0*Particle::L, 2.0*Particle::L);
+          Eigen::Vector3d shift (H.spacing(0)/2.0 - Particle::L,
+                                 H.spacing(1)/2.0 - Particle::L,
+                                 H.spacing(2)/2.0 - Particle::L);
+          T_s2g = H.transform() * newspacing;
+          T_s2g = T_s2g.inverse().translate(shift);
+        }
+
+
         void ParticleGrid::add(const Point_t &pos, const Point_t &dir)
         {
           Particle* p = pool.create(pos, dir);
           size_t gidx = pos2idx(pos);
           grid[gidx].push_back(p);
         }
-        
+
         void ParticleGrid::shift(Particle *p, const Point_t& pos, const Point_t& dir)
         {
           size_t gidx0 = pos2idx(p->getPosition());
@@ -38,27 +56,27 @@ namespace MR {
           p->setDirection(dir);
           grid[gidx1].push_back(p);
         }
-        
+
         void ParticleGrid::remove(Particle* p)
         {
           size_t gidx0 = pos2idx(p->getPosition());
           grid[gidx0].erase(std::remove(grid[gidx0].begin(), grid[gidx0].end(), p), grid[gidx0].end());
           pool.destroy(p);
         }
-        
+
         void ParticleGrid::clear()
         {
           grid.clear();
           pool.clear();
         }
-        
+
         const ParticleGrid::ParticleVectorType* ParticleGrid::at(const ssize_t x, const ssize_t y, const ssize_t z) const
         {
           if ((x < 0) || (size_t(x) >= dims[0]) || (y < 0) || (size_t(y) >= dims[1]) || (z < 0) || (size_t(z) >= dims[2]))  // out of bounds
             return nullptr;
           return &grid[xyz2idx(x, y, z)];
         }
-        
+
         void ParticleGrid::exportTracks(Tractography::Writer<float> &writer)
         {
           std::lock_guard<std::mutex> lock (mutex);
@@ -70,7 +88,7 @@ namespace MR {
           // Loop through all unvisited particles
           for (ParticleVectorType& gridvox : grid)
           {
-            for (Particle* par0 : gridvox) 
+            for (Particle* par0 : gridvox)
             {
               par = par0;
               if (!par->isVisited())
@@ -114,7 +132,7 @@ namespace MR {
             }
           }
         }
-        
+
 
       }
     }
diff --git a/src/dwi/tractography/GT/particlegrid.h b/src/dwi/tractography/GT/particlegrid.h
index d0e86bb8e5..fadd63bbaf 100644
--- a/src/dwi/tractography/GT/particlegrid.h
+++ b/src/dwi/tractography/GT/particlegrid.h
@@ -32,62 +32,45 @@ namespace MR {
   namespace DWI {
     namespace Tractography {
       namespace GT {
-        
+
         /**
          * @brief The ParticleGrid class
          */
         class ParticleGrid
         { MEMALIGN(ParticleGrid)
         public:
-          
+
           using ParticleVectorType = vector<Particle*>;
-          
-          template <class HeaderType>
-          ParticleGrid(const HeaderType& image)
-          {
-            DEBUG("Initialise particle grid.");
-            dims[0] = Math::ceil<size_t>( image.size(0) * image.spacing(0) / (2.0*Particle::L) );
-            dims[1] = Math::ceil<size_t>( image.size(1) * image.spacing(1) / (2.0*Particle::L) );
-            dims[2] = Math::ceil<size_t>( image.size(2) * image.spacing(2) / (2.0*Particle::L) );
-            grid.resize(dims[0]*dims[1]*dims[2]);
-            
-            // Initialise scanner-to-grid transform
-            Eigen::DiagonalMatrix<default_type, 3> newspacing (2.0*Particle::L, 2.0*Particle::L, 2.0*Particle::L);
-            Eigen::Vector3d shift (image.spacing(0)/2.0 - Particle::L,
-                                  image.spacing(1)/2.0 - Particle::L,
-                                  image.spacing(2)/2.0 - Particle::L);
-            T_s2g = image.transform() * newspacing;
-            T_s2g = T_s2g.inverse().translate(shift);
-          }
-          
+
+          ParticleGrid(const Header&);
           ParticleGrid(const ParticleGrid&) = delete;
           ParticleGrid& operator=(const ParticleGrid&) = delete;
-          
+
           ~ParticleGrid() {
             clear();
           }
-          
+
           inline unsigned int getTotalCount() const {
             return pool.size();
           }
-          
+
           void add(const Point_t& pos, const Point_t& dir);
-          
+
           void shift(Particle* p, const Point_t& pos, const Point_t& dir);
-          
+
           void remove(Particle* p);
-          
+
           void clear();
-          
+
           const ParticleVectorType* at(const ssize_t x, const ssize_t y, const ssize_t z) const;
-          
+
           inline Particle* getRandom() {
             return pool.random();
           }
-          
+
           void exportTracks(Tractography::Writer<float>& writer);
-          
-          
+
+
         protected:
           std::mutex mutex;
           ParticlePool pool;
@@ -95,15 +78,15 @@ namespace MR {
           Math::RNG rng;
           transform_type T_s2g;
           size_t dims[3];
-          
-          
+
+
           inline size_t pos2idx(const Point_t& pos) const
           {
             size_t x, y, z;
             pos2xyz(pos, x, y, z);
             return xyz2idx(x, y, z);
           }
-          
+
         public:
           inline void pos2xyz(const Point_t& pos, size_t& x, size_t& y, size_t& z) const
           {
@@ -113,13 +96,13 @@ namespace MR {
             y = Math::round<size_t>(gpos[1]);
             z = Math::round<size_t>(gpos[2]);
           }
-          
+
         protected:
           inline size_t xyz2idx(const size_t x, const size_t y, const size_t z) const
           {
             return z + dims[2] * (y + dims[1] * x);
           }
-          
+
         };
 
       }
diff --git a/src/dwi/tractography/algorithms/tensor_det.h b/src/dwi/tractography/algorithms/tensor_det.h
index 44fc1781a0..870401d118 100644
--- a/src/dwi/tractography/algorithms/tensor_det.h
+++ b/src/dwi/tractography/algorithms/tensor_det.h
@@ -65,7 +65,7 @@ namespace MR
           properties["method"] = "TensorDet";
 
           try {
-            auto grad = DWI::get_DW_scheme (source);
+            auto grad = DWI::get_DW_scheme (source_header);
             auto bmat_double = grad2bmatrix<double> (grad);
             binv = Math::pinv (bmat_double).cast<float>();
             bmat = bmat_double.cast<float>();
diff --git a/src/dwi/tractography/tracking/shared.cpp b/src/dwi/tractography/tracking/shared.cpp
index 7109a9d2bd..01b8034a52 100644
--- a/src/dwi/tractography/tracking/shared.cpp
+++ b/src/dwi/tractography/tracking/shared.cpp
@@ -29,7 +29,8 @@ namespace MR
 
 
         SharedBase::SharedBase (const std::string& diff_path, Properties& property_set) :
-            source (Image<float>::open (diff_path).with_direct_io (3)),
+            source_header (Header::open (diff_path)),
+            source (source_header.get_image<float>().with_direct_io (3)),
             properties (property_set),
             init_dir ({ NaN, NaN, NaN }),
             min_num_points_preds (0),
@@ -63,7 +64,7 @@ namespace MR
           properties.set (rk4, "rk4");
           properties.set (stop_on_all_include, "stop_on_all_include");
 
-          properties["source"] = source.name();
+          properties["source"] = source_header.name();
 
           max_num_seeds = Defaults::seed_to_select_ratio * max_num_tracks;
           properties.set (max_num_seeds, "max_num_seeds");
diff --git a/src/dwi/tractography/tracking/shared.h b/src/dwi/tractography/tracking/shared.h
index e2bb586e76..911e4a44f4 100644
--- a/src/dwi/tractography/tracking/shared.h
+++ b/src/dwi/tractography/tracking/shared.h
@@ -56,6 +56,7 @@ namespace MR
             virtual ~SharedBase();
 
 
+            Header source_header;
             Image<float> source;
             Properties& properties;
             Eigen::Vector3f init_dir;
@@ -76,7 +77,7 @@ namespace MR
 
             float vox () const
             {
-              return std::pow (source.spacing(0)*source.spacing(1)*source.spacing(2), float (1.0/3.0));
+              return std::pow (source_header.spacing(0)*source_header.spacing(1)*source_header.spacing(2), float (1.0/3.0));
             }
 
             void set_step_and_angle (const float stepsize, const float angle, bool is_higher_order);