diff --git a/Doxyfile b/Doxyfile
index eebb7a8881..3a291637ee 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -31,7 +31,7 @@ PROJECT_NAME           = MRtrix
 # This could be handy for archiving the generated documentation or 
 # if some version control system is used.
 
-PROJECT_NUMBER         = 0.3.2
+PROJECT_NUMBER         = 0.3.3
 
 # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) 
 # base path where the generated documentation will be put. 
diff --git a/cmd/bogus.cpp b/cmd/bogus.cpp
index 602f7716f6..5bcd937c86 100644
--- a/cmd/bogus.cpp
+++ b/cmd/bogus.cpp
@@ -44,31 +44,31 @@ typedef float T;
 
 EXECUTE {
   Image::Header header;
-  header.axes.resize (3);
-  header.axes[0].dim = 1024;
-  header.axes[1].dim = 1024;
-  header.axes[2].dim = 1024;
+  header.axes.ndim() = 3;
+  header.axes.dim(0) = 1024;
+  header.axes.dim(1) = 1024;
+  header.axes.dim(2) = 1024;
 
-  header.axes[0].vox = 1.0;
-  header.axes[1].vox = 1.0;
-  header.axes[2].vox = 1.0;
+  header.axes.vox(0) = 1.0;
+  header.axes.vox(1) = 1.0;
+  header.axes.vox(2) = 1.0;
 
-  header.axes[0].order = 0;
-  header.axes[1].order = 1;
-  header.axes[2].order = 2;
+  header.axes.order(0) = 0;
+  header.axes.order(1) = 1;
+  header.axes.order(2) = 2;
 
-  VAR (header.data_type.description());
+  VAR (header.datatype().description());
   //header.data_type = DataType::UInt8;
 
-  Image::Object obj;
+  const Image::Header obj;
   obj.create ("poo.mif", header);
   
-  Image::Voxel vox (obj);
-  vox[0] = 1023;
-  vox[1] = 1023;
-  vox[2] = 1023;
+  Image::Voxel<float> vox (obj);
+  vox.pos(0, 1023);
+  vox.pos(1, 1023);
+  vox.pos(2, 1023);
 
-  vox.value() = 0.0;
+  vox.value (0.0);
 
 
   /*
diff --git a/cmd/csdeconv.cpp b/cmd/csdeconv.cpp
index 4113a4d7d6..9e31e45a5e 100644
--- a/cmd/csdeconv.cpp
+++ b/cmd/csdeconv.cpp
@@ -110,7 +110,7 @@ class Item {
 class Allocator {
   public:
     Allocator (size_t data_size) : N (data_size) { }
-    Item* alloc () { Item* item = new Item; item->data.resize (N); return (item); }
+    Item* alloc () { Item* item = new Item; item->data.allocate (N); return (item); }
     void reset (Item* item) { }
     void dealloc (Item* item) { delete item; }
   private:
@@ -163,23 +163,23 @@ class DataLoader {
           value_type norm = 0.0;
           if (P.normalise) {
             for (size_t n = 0; n < P.bzeros.size(); n++) {
-              D.pos(3, P.bzeros[n]);
+              D[3] = P.bzeros[n];
               norm += D.value ();
             }
             norm /= P.bzeros.size();
           }
 
           for (size_t n = 0; n < P.dwis.size(); n++) {
-            D.pos(3, P.dwis[n]);
+            D[3] = P.dwis[n];
             item->data[n] = D.value(); 
             if (!finite (item->data[n])) return;
             if (item->data[n] < 0.0) item->data[n] = 0.0;
             if (P.normalise) item->data[n] /= norm;
           }
 
-          item->pos[0] = D.pos(0);
-          item->pos[1] = D.pos(1);
-          item->pos[2] = D.pos(2);
+          item->pos[0] = D[0];
+          item->pos[1] = D[1];
+          item->pos[2] = D[2];
 
           if (!item.write()) throw Exception ("error writing to work queue");
         }
@@ -213,12 +213,12 @@ class Processor {
             str(item->pos[0]) + " " + str(item->pos[1]) + " " + str(item->pos[2]) +
             " ] failed to converge"); 
 
-        SH.pos(0, item->pos[0]);
-        SH.pos(1, item->pos[1]);
-        SH.pos(2, item->pos[2]);
+        SH[0] = item->pos[0];
+        SH[1] = item->pos[1];
+        SH[2] = item->pos[2];
 
-        for (SH.pos(3,0); SH.pos(3) < SH.dim(3); SH.move(3,1))
-          SH.value (sdeconv.FOD()[SH.pos(3)]);
+        for (SH[3] = 0; SH[3] < SH.dim(3); ++SH[3])
+          SH.value() = sdeconv.FOD()[SH[3]];
       }
     }
 
@@ -251,7 +251,7 @@ EXECUTE {
   else {
     if (!header.DW_scheme.is_set()) 
       throw Exception ("no diffusion encoding found in image \"" + header.name() + "\"");
-    grad = header.DW_scheme;
+    grad.copy (header.DW_scheme);
   }
 
   if (grad.rows() < 7 || grad.columns() != 4) 
@@ -301,15 +301,15 @@ EXECUTE {
   if (opt.size()) HR_dirs.load (opt[0][0].get_string());
   else {
     HR_dirs.allocate (300,2);
-    HR_dirs.view() = Math::MatrixView<float> (default_directions, 300, 2);
+    HR_dirs = Math::Matrix<float> (default_directions, 300, 2);
   }
 
   header.axes.dim(3) = Math::SH::NforL (lmax);
   header.datatype() = DataType::Float32;
-  header.axes.order(0) = 1; header.axes.forward(0) = true;
-  header.axes.order(1) = 2; header.axes.forward(1) = true;
-  header.axes.order(2) = 3; header.axes.forward(2) = true;
-  header.axes.order(3) = 0; header.axes.forward(3) = true;
+  header.axes.stride(0) = 2;
+  header.axes.stride(1) = 3;
+  header.axes.stride(2) = 4;
+  header.axes.stride(3) = 1;
 
   const Image::Header* mask_header = NULL;
   opt = get_options (2); // mask
diff --git a/cmd/gendir.cpp b/cmd/gendir.cpp
index 3360a20730..a014b70dd1 100644
--- a/cmd/gendir.cpp
+++ b/cmd/gendir.cpp
@@ -131,7 +131,7 @@ EXECUTE {
   ProgressBar::init (0, "Optimising directions");
   for (power = -1.0; power >= -target_power/2.0; power *= 2.0) {
     info ("setting power = " + str (-power*2.0));
-    gsl_multimin_fdfminimizer_set (minimizer, &fdf, &v, 0.01, 1e-4);
+    gsl_multimin_fdfminimizer_set (minimizer, &fdf, v.gsl(), 0.01, 1e-4);
 
     for (uint iter = 0; iter < niter; iter++) {
 
@@ -147,7 +147,7 @@ EXECUTE {
 
       ProgressBar::inc();
     }
-    gsl_vector_memcpy (&v, minimizer->x);
+    gsl_vector_memcpy (v.gsl(), minimizer->x);
   }
   ProgressBar::done();
 
diff --git a/cmd/mrconvert.cpp b/cmd/mrconvert.cpp
index d1a0b1d382..97b8a5edc2 100644
--- a/cmd/mrconvert.cpp
+++ b/cmd/mrconvert.cpp
@@ -85,7 +85,7 @@ OPTIONS = {
 template <class Set, class Set2> void copy_replace_NaN_kernel (Set& destination, Set2& source) { 
   typedef typename Set::value_type T;
   T val = source.value();
-  destination.value (isnan(val) ? 0.0 : val);
+  destination.value() = ( isnan(val) ? 0.0 : val );
 }
 
 
@@ -93,20 +93,18 @@ template <class Set, class Set2> void copy (Set& destination, Set2& source, bool
 { 
   std::string progress_message ("copying from \"" + source.name() + "\" to \"" + destination.name() + "\"...");
 
-  typedef DataSet::Reorder<Set> S1;
-  typedef DataSet::Reorder<Set2> S2;
+  typedef DataSet::Reorder<Set2> SetR;
 
-  S1 dest (destination, NULL, destination.name());
-  S2 src (source, dest.layout(), source.name());
+  SetR src (source, destination, source.name());
 
-  if (replace_NaN) DataSet::loop2 (progress_message, DataSet::copy_kernel<S1,S2>, dest, src);
-  else DataSet::loop2 (progress_message, copy_replace_NaN_kernel<S1,S2>, dest, src);
+  if (replace_NaN) DataSet::loop2 (progress_message, DataSet::copy_kernel<Set,SetR>, destination, src);
+  else DataSet::loop2 (progress_message, copy_replace_NaN_kernel<Set,SetR>, destination, src);
 }
 
 
 
 EXECUTE {
-  std::vector<OptBase> opt = get_options (1); // vox
+  OptionList opt = get_options (1); // vox
   std::vector<float> vox;
   if (opt.size()) 
     vox = parse_floats (opt[0][0].get_string());
@@ -158,14 +156,12 @@ EXECUTE {
 
   opt = get_options (7); // layout
   if (opt.size()) {
-    std::vector<Image::Axes::Order> ax = parse_axes_specifier (header.axes, opt[0][0].get_string());
+    std::vector<ssize_t> ax = Image::Axes::parse (header.ndim(), opt[0][0].get_string());
     if (ax.size() != header.axes.ndim()) 
       throw Exception (std::string("specified layout \"") + opt[0][0].get_string() + "\" does not match image dimensions");
 
-    for (size_t i = 0; i < ax.size(); i++) {
-      header.axes.order(i) = ax[i].order;
-      header.axes.forward(i) = ax[i].forward;
-    }
+    for (size_t i = 0; i < ax.size(); i++)
+      header.axes.stride(i) = ax[i];
   }
 
 
@@ -219,6 +215,3 @@ EXECUTE {
 }
 
 
-
-
-
diff --git a/cmd/mrmath.cpp b/cmd/mrmath.cpp
new file mode 100644
index 0000000000..d73af05f7f
--- /dev/null
+++ b/cmd/mrmath.cpp
@@ -0,0 +1,282 @@
+/*
+    Copyright 2009 Brain Research Institute, Melbourne, Australia
+
+    Written by J-Donald Tournier, 03/12/09.
+
+    This file is part of MRtrix.
+
+    MRtrix is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    MRtrix is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with MRtrix.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+#include "app.h"
+#include "image/voxel.h"
+#include "math/rng.h"
+#include "dataset/loop.h"
+
+using namespace MR; 
+
+SET_VERSION_DEFAULT;
+
+DESCRIPTION = {
+  "apply generic mathematical operations to images.",
+  NULL
+};
+
+
+ARGUMENTS = { 
+  Argument ("image", "input image", "the input image.").type_image_in (),
+  Argument ("output", "output image", "the output image.").type_image_out (),
+  Argument::End 
+}; 
+
+OPTIONS = { 
+  Option ("datatype", "data type", "specify output image data type.") // 0
+    .append (Argument ("spec", "specifier", "the data type specifier.").type_choice (DataType::identifiers)),
+
+  Option ("abs", "absolute value", "take absolute value of current voxel value"), // 1
+  Option ("neg", "negative value", "take negative value of current voxel value"), // 2
+  Option ("sqrt", "square root", "take square root of current voxel value"), // 3
+  Option ("exp", "exp", "take e raised to the power of current voxel value"), // 4
+  Option ("log", "log", "take natural logarithm of current voxel value"), // 5
+  Option ("cos", "cos", "take cosine of current voxel value"), // 6
+  Option ("sin", "sin", "take sine of current voxel value"), // 7
+  Option ("tan", "tanh", "take tangent of current voxel value"), // 8
+  Option ("cosh", "cosh", "take hyperbolic cosine of current voxel value"), // 9
+  Option ("sinh", "sinh", "take hyperbolic sine of current voxel value"), // 10
+  Option ("tanh", "tanh", "take hyperbolic tangent of current voxel value"), // 11
+  Option ("acos", "acos", "take inverse cosine of current voxel value"), // 12
+  Option ("asin", "asin", "take inverse sine of current voxel value"), // 13
+  Option ("atan", "atan", "take inverse tangent of current voxel value"), // 14
+  Option ("acosh", "acosh", "take inverse hyperbolic cosine of current voxel value"), // 15
+  Option ("asinh", "asinh", "take inverse hyperbolic sine of current voxel value"), // 16
+  Option ("atanh", "atanh", "take inverse hyperbolic tangent of current voxel value"), // 17
+
+  Option::End 
+};
+
+typedef double value_type;
+
+std::vector<ssize_t> dimension;
+
+void check_dim (const Image::Header& header) {
+  try {
+    if (header.ndim() != dimension.size()) throw 0;
+    for (size_t n = 0; n < dimension.size(); ++n) 
+      if (header.dim(n) != dimension[n]) throw 0;
+    return;
+  }
+  catch (...) {
+    throw Exception ("mismatch between dimensions of input images");
+  }
+}
+
+class Functor {
+  public:
+    Functor (Functor* input) : in (input) { }
+    virtual ~Functor () { }
+    virtual void get (const std::vector<ssize_t>& pos, std::vector<value_type>& values) { }
+    virtual size_t ndim () const { return (0); }
+    virtual ssize_t dim (size_t axis) const { return (0); }
+  protected:
+    Ptr<Functor> in;
+};
+
+
+class Source : public Functor {
+  public:
+    Source (const Image::Header& header) :
+      Functor (NULL),
+      V (header) { }
+    void get (const std::vector<ssize_t>& pos, std::vector<value_type>& values) {
+      ssize_t axis = -1;
+      for (size_t n = 0; n < V.ndim(); ++n) {
+        if (pos[n] < 0) axis = n;
+        else V.pos(n,pos[n]);
+      }
+      assert (axis >= 0);
+      assert (size_t (V.dim(axis)) == values.size());
+      for (V.pos(axis,0); V.pos(axis) < V.dim(axis); V.move(axis,1)) 
+        values[V.pos(axis)] = V.value();
+    }
+    virtual size_t ndim () const { return (V.ndim()); }
+    virtual ssize_t dim (size_t axis) const { return (V.dim (axis)); }
+  private:
+    Image::Voxel<value_type> V;
+};
+
+
+class Constant : public Functor {
+  public:
+    Constant (value_type value) :
+      Functor (NULL),
+      val (value) { }
+    virtual void get (const std::vector<ssize_t>& pos, std::vector<value_type>& values) {
+      for (size_t n = 0; n < values.size(); ++n) values[n] = val;
+    }
+    virtual size_t ndim () const { return (0); }
+    virtual ssize_t dim (size_t axis) const { return (0); }
+  private:
+    value_type val;
+};
+
+
+class Unary : public Functor {
+  public:
+    Unary (Functor* input) : Functor (input) { }
+    virtual size_t ndim () const { return (in->ndim()); }
+    virtual ssize_t dim (size_t axis) const { return (in->dim (axis)); }
+};
+
+class Binary : public Functor {
+  public:
+    Binary (Functor* input1, Functor* input2) : Functor (input1), in2 (input2) { }
+    virtual size_t ndim () const { assert (in->ndim() == in2->ndim()); return (in->ndim()); }
+    virtual ssize_t dim (size_t axis) const { return (in->dim (axis) > 1 ? in->dim(axis) : in2->dim(axis)); }
+  protected:
+    Ptr<Functor> in2;
+    std::vector<value_type> values2;
+};
+
+
+class Kernel {
+  public:
+    Kernel (Image::Voxel<value_type>& output, Functor* input) :
+      out (output), in (input), N (in->ndim()), x (in->ndim()) { 
+        for (size_t n = 0; n < N.size(); ++n) {
+          x[n] = 0;
+          N[n] = in->dim(n);
+        }
+        x[0] = -1;
+        values.resize (N[0]);
+      }
+    ~Kernel () { ProgressBar::done(); }
+    std::string name () { return ("math kernel"); }
+    size_t ndim () const { return (N.size()); }
+    ssize_t dim (size_t axis) const { return (N[axis]); }
+    ssize_t pos (size_t axis) const { return (x[axis]); }
+    void pos (size_t axis, ssize_t position) { x[axis] = position; out.pos (axis, position); }
+    void move (size_t axis, ssize_t inc) { x[axis] += inc; out.move (axis, inc); }
+    void check (size_t from_axis, size_t to_axis) const { 
+      if (!DataSet::dimensions_match (out, *this))
+        throw Exception ("dimensions mismatch in math kernel");
+      ProgressBar::init (DataSet::voxel_count (*this, 1, SIZE_MAX), "performing mathematical operations...");
+    }
+    void operator() () { 
+      in->get (x, values);
+      for (out.pos(0,0); out.pos(0) < N[0]; out.move(0,1)) 
+        out.value (values[out.pos(0)]);
+      ProgressBar::inc();
+    }
+  protected:
+    Image::Voxel<value_type>& out;
+    Functor* in;
+    std::vector<ssize_t> N, x;
+    std::vector<value_type> values;
+};
+
+#define DEF_UNARY(name, operation) class name : public Unary { \
+  public: \
+    name (Functor* input) : Unary (input) { } \
+    virtual void get (const std::vector<ssize_t>& pos, std::vector<value_type>& values) { \
+      in->get (pos, values); \
+      for (ssize_t n = 0; n < dimension[0]; ++n) { \
+        value_type& val = values[n]; \
+        values[n] = operation; \
+      } \
+    } \
+}
+
+#define DEF_BINARY(name, operation) class name : public Binary { \
+  public: \
+    name (Functor* input1, Functor* input2) : Binary (input1, input2) { } \
+    virtual void get (const std::vector<ssize_t>& pos, std::vector<value_type>& values) { \
+      values2.resize (values.size()); \
+      in->get (pos, values); \
+      in2->get (pos, values2); \
+      for (ssize_t n = 0; n < dimension[0]; ++n) { \
+        value_type& val1 = values[n]; \
+        value_type& val2 = values2[n]; \
+        values[n] = operation; \
+      } \
+    } \
+}
+
+DEF_UNARY (Abs, Math::abs(val));
+DEF_UNARY (Neg, -val);
+DEF_UNARY (Sqrt, Math::sqrt (val));
+DEF_UNARY (Exp, Math::exp (val));
+DEF_UNARY (Log, Math::log (val));
+DEF_UNARY (Cos, Math::cos (val));
+DEF_UNARY (Sin, Math::sin (val));
+DEF_UNARY (Tan, Math::tan (val));
+DEF_UNARY (Cosh, Math::cosh (val));
+DEF_UNARY (Sinh, Math::sinh (val));
+DEF_UNARY (Tanh, Math::tanh (val));
+DEF_UNARY (Acos, Math::acos (val));
+DEF_UNARY (Asin, Math::asin (val));
+DEF_UNARY (Atan, Math::atan (val));
+DEF_UNARY (Acosh, Math::acosh (val));
+DEF_UNARY (Asinh, Math::asinh (val));
+DEF_UNARY (Atanh, Math::atanh (val));
+
+DEF_BINARY (Mult, val1*val2);
+
+EXECUTE {
+  const Image::Header header_in (argument[0].get_image());
+  Image::Header header (header_in);
+
+  header.datatype() = DataType::Float32;
+  OptionList opt = get_options (0); // datatype
+  if (opt.size()) header.datatype().parse (DataType::identifiers[opt[0][0].get_int()]);
+
+  Image::Header header_out (argument[1].get_image (header));
+  dimension.resize (header.ndim());
+  for (size_t n = 0; n < header.ndim(); ++n) 
+    dimension[n] = header.dim(n);
+
+  Functor* last = new Source (header_in);
+
+  for (size_t n = 0; n < option.size(); ++n) {
+    switch (option[n].index) {
+      case 0: continue;
+      case 1: last = new Abs (last); break;
+      case 2: last = new Neg (last); break;
+      case 3: last = new Sqrt (last); break;
+      case 4: last = new Exp (last); break;
+      case 5: last = new Log (last); break;
+      case 6: last = new Cos (last); break;
+      case 7: last = new Sin (last); break;
+      case 8: last = new Tan (last); break;
+      case 9: last = new Cosh (last); break;
+      case 10: last = new Sinh (last); break;
+      case 11: last = new Tanh (last); break;
+      case 12: last = new Acos (last); break;
+      case 13: last = new Asin (last); break;
+      case 14: last = new Atan (last); break;
+      case 15: last = new Acosh (last); break;
+      case 16: last = new Asinh (last); break;
+      case 17: last = new Atanh (last); break;
+      default: assert (0);
+    }
+  }
+
+  Image::Voxel<value_type> vox (header_out);
+  Kernel kernel (vox, last);
+  DataSet::loop (kernel,1);
+
+  delete last;
+}
+
+
diff --git a/cmd/threshold.cpp b/cmd/threshold.cpp
index 7c275185a8..004e52c379 100644
--- a/cmd/threshold.cpp
+++ b/cmd/threshold.cpp
@@ -76,7 +76,7 @@ class ByValue {
     ByValue (float threshold, float value_if_less, float value_if_more) :
       val (threshold), lt (value_if_less), gt (value_if_more) { }
     void operator() (Image::Voxel<float>& dest, Image::Voxel<float>& src) {
-      dest.value (src.value() < val ? lt : gt);
+      dest.value() = src.value() < val ? lt : gt;
     }
   private:
     float val, lt, gt;
@@ -88,7 +88,8 @@ class ByValue {
 class ByGreatestN {
   public:
     ByGreatestN (size_t num) : N (num) { }
-    void operator() (const Image::Voxel<float>& D) {
+
+    void operator() (Image::Voxel<float>& D) {
       assert (list.size() <= N);
       float val = D.value();
       if (list.size() == N) {
@@ -96,7 +97,7 @@ class ByGreatestN {
         list.erase (list.begin());
       }
       std::vector<ssize_t> pos (D.ndim());
-      for (size_t n = 0; n < D.ndim(); ++n) pos[n] = D.pos(n);
+      for (size_t n = 0; n < D.ndim(); ++n) pos[n] = D[n];
       list.insert (std::pair<float,std::vector<ssize_t> > (val, pos));
     }
 
@@ -106,8 +107,8 @@ class ByGreatestN {
 
       for (std::multimap<float,std::vector<ssize_t> >::const_iterator i = list.begin(); i != list.end(); ++i) {
         for (size_t n = 0; n < dest.ndim(); ++n)
-          dest.pos(n, i->second[n]);
-        dest.value (val_if_included);
+          dest[n] = i->second[n];
+        dest.value() = val_if_included;
       }
     }
 
@@ -115,7 +116,7 @@ class ByGreatestN {
     class ZeroKernel {
       public:
         ZeroKernel (float value) : val (value) { }
-        void operator() (Image::Voxel<float>& D) { D.value (val); }
+        void operator() (Image::Voxel<float>& D) { D.value() = val; }
       private:
         float val;
     };
@@ -128,7 +129,7 @@ class ByGreatestN {
 class BySmallestN : public ByGreatestN {
   public:
     BySmallestN (size_t num) : ByGreatestN (num) { }
-    void operator() (const Image::Voxel<float>& D) {
+    void operator() (Image::Voxel<float>& D) {
       assert (list.size() <= N);
       float val = D.value();
       if (list.size() == N) {
@@ -138,7 +139,7 @@ class BySmallestN : public ByGreatestN {
         list.erase (i);
       }
       std::vector<ssize_t> pos (D.ndim());
-      for (size_t n = 0; n < D.ndim(); ++n) pos[n] = D.pos(n);
+      for (size_t n = 0; n < D.ndim(); ++n) pos[n] = D[n];
       list.insert (std::pair<float,std::vector<ssize_t> > (val, pos));
     }
 };
@@ -148,7 +149,7 @@ class BySmallestN : public ByGreatestN {
 
 EXECUTE {
   float val (NAN), percentile (NAN);
-  size_t topN (0), bottomN (0), nopt (0);
+  off64_t topN (0), bottomN (0), nopt (0);
 
   OptionList opt = get_options (0); // abs
   if (opt.size()) {
diff --git a/configure b/configure
index 0676233398..0c13662ed1 100755
--- a/configure
+++ b/configure
@@ -268,13 +268,11 @@ if 'LDLIB_FLAGS' in os.environ.keys(): ld_lib_flags = os.environ['LDLIB_FLAGS'].
 obj_suffix = '.o'
 exe_suffix = ''
 lib_prefix = 'lib'
-lib_suffix = '.so'
-gl_cflags = []
 
 system = platform.system().lower()
 report ('Detecting OS: ' + system + '\n')
 if system == 'linux':
-  gl_cflags = [ '-DGL_GLEXT_PROTOTYPES' ]
+  lib_suffix = '.so'
 elif system == 'windows':
   cpp_flags += [ '-DWINDOWS', '-mno-cygwin', '-mms-bitfields' ]
   exe_suffix = '.exe'
@@ -725,7 +723,6 @@ int main() { Foo f; }
 
 
   report ('Checking for OpenGL: ')
-  qt_cflags += gl_cflags
   try:
     gl_test = compile ('''
 #include <QtGui>
@@ -797,6 +794,45 @@ WARNING: not all required OpenGL extensions are supported (see above).
 
 
 
+  report ('Checking whether OpenGL extensions compile with GL_GLEXT_PROTOTYPES: ')
+  try:
+    gl_test = compile ('''
+#include <QtGui>
+#include <QtOpenGL>
+#include <QApplication>
+#include <QGLWidget>
+#include <iostream>
+#include <GL/glext.h>
+
+class GLWidget : public QGLWidget {
+  public:
+    GLWidget (QWidget *parent = 0) : QGLWidget (parent) { }
+    ~GLWidget() { }
+  protected:
+    void initializeGL () { exit (0); }
+    void paintGL () { 
+      char* log;
+      int chars;
+      GLhandleARB obj;
+      glGetInfoLogARB (obj, 1, &chars, log); 
+    }
+    void resizeGL (int width, int height) { }
+};
+
+int main (int argc, char *argv[]) {
+  QApplication app (argc, argv);
+  GLWidget window;
+  window.show();
+  return app.exec();
+}
+''', qt_cflags + [ '-DGL_GLEXT_PROTOTYPES' ], qt_ldflags)
+    qt_cflags += [ '-DGL_GLEXT_PROTOTYPES' ]
+    report ('yes\n')
+  except:
+    report ('no\n')
+
+
+
 
 
 
diff --git a/lib/app.cpp b/lib/app.cpp
index 38021c49e9..ca5afb15e8 100644
--- a/lib/app.cpp
+++ b/lib/app.cpp
@@ -168,6 +168,9 @@ namespace MR {
     if (candidates.size() == 0) return (SIZE_MAX);
     if (candidates.size() == 1) return (candidates[0]);
 
+    for (size_t n = 0; n < candidates.size(); ++n) 
+      if (option_name (candidates[n]) == s) return (candidates[n]);
+
     s = "several matches possible for option \"-" + s + "\": \"-" + option_name (candidates[0]) + "\", \"-" + option_name (candidates[1]) + "\"";
     for (size_t n = 2; n < candidates.size(); n++) { s += ", \"-"; s += option_name (candidates[n]); s += "\""; }
     throw Exception (s);
diff --git a/lib/app.h b/lib/app.h
index cbcf95dfa6..ad57e0669c 100644
--- a/lib/app.h
+++ b/lib/app.h
@@ -100,7 +100,7 @@ namespace MR {
       std::vector<OptBase> option;
 
       virtual void execute () = 0;
-      std::vector<OptBase> get_options (size_t index);
+      OptionList get_options (size_t index);
 
       void parse_arguments ();
 
@@ -139,9 +139,9 @@ namespace MR {
 
 
 
-  inline std::vector<OptBase> App::get_options (size_t index)
+  inline OptionList App::get_options (size_t index)
   {
-    std::vector<OptBase> a;
+    OptionList a;
     for (size_t n = 0; n < option.size(); n++) 
       if (option[n].index == index)
         a.push_back (option[n]);
diff --git a/lib/dataset/copy.h b/lib/dataset/copy.h
index 9683a9451f..96216a38aa 100644
--- a/lib/dataset/copy.h
+++ b/lib/dataset/copy.h
@@ -32,7 +32,7 @@ namespace MR {
     // @{
 
     template <class Set, class Set2> 
-      void copy_kernel (Set& destination, Set2& source) { destination.value (source.value()); }
+      void copy_kernel (Set& destination, Set2& source) { destination.value() = source.value(); }
 
     template <class Set, class Set2> 
       void copy (Set& destination, Set2& source, size_t from_axis = 0, size_t to_axis = SIZE_MAX) { 
diff --git a/lib/dataset/extract.h b/lib/dataset/extract.h
index afa3306f29..0ef795470a 100644
--- a/lib/dataset/extract.h
+++ b/lib/dataset/extract.h
@@ -24,6 +24,7 @@
 #define __dataset_extract_h__
 
 #include "math/matrix.h"
+#include "dataset/value.h"
 
 namespace MR {
   namespace DataSet {
@@ -42,21 +43,26 @@ namespace MR {
         int     dim (size_t axis) const { return (P[axis].size()); }
         float   vox (size_t axis) const { return (D.vox (axis)); }
         const Math::Matrix<float>& transform () const { return (D.transform()); }
-        const size_t* layout () const { return (D.layout()); }
+        const std::vector<ssize_t>& stride () const { return (D.stride()); }
 
         void reset () { memset (x, 0, sizeof(size_t)*ndim()); for (size_t a = 0; a < ndim(); ++a) D.pos (a, P[a][0]); }
 
-        ssize_t pos (size_t axis) const { return (x[axis]); } 
-        void    pos (size_t axis, ssize_t position) const { x[axis] = position; D.pos (axis, P[axis][position]); }
-        void    move (size_t axis, ssize_t increment) const { x[axis] += increment; D.pos (axis, P[axis][x[axis]]); }
-
-        value_type   value () const { return (D.value()); }
-        void         value (value_type val) { D.value (val); }
+        Position<Extract<Set> > operator[] (size_t axis) { return (Position<Extract<Set> > (*this, axis)); }
+        Value<Extract<Set> > value () { return (Value<Extract<Set> > (*this)); }
 
       private:
         Set& D;
         size_t* x;
         const std::vector<std::vector<int> > P;
+
+        value_type get_value () const { return (D.value()); }
+        void set_value (value_type val) { D.value (val); }
+        ssize_t get_pos (size_t axis) const { return (x[axis]); } 
+        void set_pos (size_t axis, ssize_t position) const { x[axis] = position; D[axis] = P[axis][position]; }
+        void move_pos (size_t axis, ssize_t increment) const { x[axis] += increment; D[axis] += P[axis][x[axis]]; }
+
+        friend class Position<Extract<Set> >;
+        friend class Value<Extract<Set> >;
     };
     
     //! @}
diff --git a/lib/dataset/interp.h b/lib/dataset/interp.h
index 8f2a1c4953..a0b77b2088 100644
--- a/lib/dataset/interp.h
+++ b/lib/dataset/interp.h
@@ -96,14 +96,14 @@ namespace MR {
         value_type value () const {
           if (out_of_bounds) return (NAN);
           value_type val = 0.0;
-          if (faaa) val  = faaa * data.value (); data.move(2,1);
-          if (faab) val += faab * data.value (); data.move(1,1);
-          if (fabb) val += fabb * data.value (); data.move(2,-1);
-          if (faba) val += faba * data.value (); data.move(0,1);
-          if (fbba) val += fbba * data.value (); data.move(1,-1);
-          if (fbaa) val += fbaa * data.value (); data.move(2,1);
-          if (fbab) val += fbab * data.value (); data.move(1,1);
-          if (fbbb) val += fbbb * data.value (); data.move(0,-1); data.move(1,-1); data.move(2,-1);
+          if (faaa) val  = faaa * data.value(); data[2]++;
+          if (faab) val += faab * data.value(); data[1]++;
+          if (fabb) val += fabb * data.value(); data[2]--;
+          if (faba) val += faba * data.value(); data[0]++;
+          if (fbba) val += fbba * data.value(); data[1]--;
+          if (fbaa) val += fbaa * data.value(); data[2]++;
+          if (fbab) val += fbab * data.value(); data[1]++;
+          if (fbbb) val += fbbb * data.value(); data[0]--;-1 data[1]--; data[2]--;
           return (val);
         }
 
diff --git a/lib/dataset/loop.h b/lib/dataset/loop.h
index bbdf577b7a..21bcf6400b 100644
--- a/lib/dataset/loop.h
+++ b/lib/dataset/loop.h
@@ -24,6 +24,7 @@
 #define __dataset_loop_h__
 
 #include "progressbar.h"
+#include "dataset/position.h"
 
 namespace MR {
   namespace DataSet {
@@ -31,24 +32,20 @@ namespace MR {
     //! \cond skip
     namespace {
 
-      template <class Set> bool increment (Set& D, size_t axis, size_t last_axis) { 
-        if (axis >= last_axis) return (false);
-        if (D.pos(axis)+1 < ssize_t(D.dim(axis))) { D.move(axis,1); return (true); }
-        D.pos (axis,0); 
-        return (increment (D, axis+1, last_axis));
-      }
-
       template <class Functor, class Set> class Kernel1 {
         public:
           Kernel1 (Functor& func, Set& set) : F (func), D (set) { assert (voxel_count (D)); }
           const std::string& name () { return (D.name()); }
           size_t ndim () const { return (D.ndim()); }
           ssize_t dim (size_t axis) const { return (D.dim (axis)); }
-          ssize_t pos (size_t axis) const { return (D.pos (axis)); }
-          void pos (size_t axis, ssize_t position) { D.pos (axis, position); }
-          void move (size_t axis, ssize_t inc) { D.move (axis, inc); }
+          Position<Kernel1<Functor,Set> > operator[] (size_t axis) {
+            return (Position<Kernel1<Functor,Set> > (*this, axis)); }
           void check (size_t from_axis, size_t to_axis) { }
           void operator() () { F (D); }
+
+          ssize_t get_pos (size_t axis) const { return (D[axis]); }
+          void set_pos (size_t axis, ssize_t position) { D[axis] = position; }
+          void move_pos (size_t axis, ssize_t inc) { D[axis] += inc; }
         private:
           Functor& F;
           Set& D;
@@ -60,14 +57,17 @@ namespace MR {
           const std::string& name () { return (D.name()); }
           size_t ndim () const { return (D.ndim()); }
           ssize_t dim (size_t axis) const { return (D.dim (axis)); }
-          ssize_t pos (size_t axis) const { return (D.pos (axis)); }
-          void pos (size_t axis, ssize_t position) { D.pos (axis, position); D2.pos (axis, position); }
-          void move (size_t axis, ssize_t inc) { D.move (axis, inc); D2.move (axis, inc); }
+          Position<Kernel2<Functor,Set,Set2> > operator[] (size_t axis) {
+            return (Position<Kernel2<Functor,Set,Set2> > (*this, axis)); }
           void check (size_t from_axis, size_t to_axis) const { 
             if (!dimensions_match (D, D2, from_axis, to_axis))
               throw Exception ("dimensions mismatch between \"" + D.name() + "\" and \"" + D2.name() + "\"");
           }
           void operator() () { F (D, D2); }
+
+          ssize_t get_pos (size_t axis) const { return (D[axis]); }
+          void set_pos (size_t axis, ssize_t position) { D[axis] = position; D2[axis] = position; }
+          void move_pos (size_t axis, ssize_t inc) { D[axis] += inc; D2[axis] += inc; }
         private:
           Functor& F;
           Set& D;
@@ -81,19 +81,22 @@ namespace MR {
           const std::string& name () { return (M.name()); }
           size_t ndim () const { return (M.ndim()); }
           ssize_t dim (size_t axis) const { return (M.dim (axis)); }
-          ssize_t pos (size_t axis) const { return (M.pos (axis)); }
-          void pos (size_t axis, ssize_t position) { M.pos (axis, position); }
-          void move (size_t axis, ssize_t inc) { M.move (axis, inc); }
+          Position<KernelMask<FunctorKernel,Mask> > operator[] (size_t axis) {
+            return (Position<KernelMask<FunctorKernel,Mask> > (*this, axis)); }
           void check (size_t from_axis, size_t to_axis) const { 
             if (!dimensions_match (M, F, from_axis, to_axis))
               throw Exception ("dimensions mismatch between \"" + M.name() + "\" and \"" + F.name() + "\"");
           }
           void operator() () { 
             if (M.value()) { 
-              for (size_t i = 0; i < M.ndim(); ++i) F.pos (i, M.pos(i));
+              for (size_t i = 0; i < M.ndim(); ++i) F[i] = M[i];
               F ();
             }
           }
+
+          ssize_t get_pos (size_t axis) const { return (M[axis]); }
+          void set_pos (size_t axis, ssize_t position) { M[axis] = position; }
+          void move_pos (size_t axis, ssize_t inc) { M[axis] += inc; }
         private:
           FunctorKernel& F;
           Mask& M;
@@ -107,9 +110,8 @@ namespace MR {
           const std::string& name () { return (F.name()); }
           size_t ndim () const { return (F.ndim()); }
           ssize_t dim (size_t axis) const { return (F.dim (axis)); }
-          ssize_t pos (size_t axis) const { return (F.pos (axis)); }
-          void pos (size_t axis, ssize_t position) { F.pos (axis, position); }
-          void move (size_t axis, ssize_t inc) { F.move (axis, inc); }
+          Position<ProgressKernel<Functor> > operator[] (size_t axis) {
+            return (Position<ProgressKernel<Functor> > (*this, axis)); }
           void check (size_t from_axis, size_t to_axis) const { 
             F.check (from_axis, to_axis);
             size_t count = 1;
@@ -117,6 +119,10 @@ namespace MR {
             ProgressBar::init (count, m);
           }
           void operator() () { F(); ProgressBar::inc(); }
+
+          ssize_t get_pos (size_t axis) const { return (F[axis]); }
+          void set_pos (size_t axis, ssize_t position) { F[axis] = position; }
+          void move_pos (size_t axis, ssize_t inc) { F[axis] += inc; }
         private:
           Functor& F;
           const std::string& m;
@@ -127,11 +133,11 @@ namespace MR {
         inline void loop_exec (Kernel& K, size_t from_axis, size_t to_axis) {
           --to_axis;
           assert (K.dim(to_axis));
-          K.pos (to_axis,0);
+          K[to_axis] = 0;
 loop_start:
           if (to_axis == from_axis) K();
           else loop_exec (K, from_axis, to_axis);
-          if (K.pos(to_axis) < K.dim(to_axis)-1) { K.move (to_axis,1); goto loop_start; }
+          if (K[to_axis] < K.dim(to_axis)-1) { ++K[to_axis]; goto loop_start; }
         }
 
     }
@@ -140,6 +146,17 @@ namespace MR {
     //! \addtogroup DataSet
     // @{
 
+
+    template <class Set> bool increment_position (Set& D, size_t from_axis, size_t to_axis) { 
+      if (from_axis >= to_axis) return (false);
+      if (D[from_axis]+1 < D.dim(from_axis)) { ++D[from_axis]; return (true); }
+      D[from_axis] = 0; 
+      return (increment_position (D, from_axis+1, to_axis));
+    }
+
+    template <class Set> inline bool increment_position (Set& D) { return (increment_position (D, 0, D.ndim())); }
+
+
     template <class Kernel> 
       void loop (Kernel& K, size_t from_axis = 0, size_t to_axis = SIZE_MAX) 
       {
diff --git a/lib/dataset/misc.h b/lib/dataset/misc.h
index d0fad72c3f..3d529f8fd3 100644
--- a/lib/dataset/misc.h
+++ b/lib/dataset/misc.h
@@ -33,11 +33,12 @@ namespace MR {
     // @{
 
     //! returns the number of voxel in the data set, or a relevant subvolume
-    template <class Set> inline off64_t voxel_count (const Set& ds, size_t up_to_dim = SIZE_MAX) 
+    template <class Set> inline off64_t voxel_count (const Set& ds, size_t from_axis = 0, size_t to_axis = SIZE_MAX) 
     { 
-      if (up_to_dim > ds.ndim()) up_to_dim = ds.ndim();
+      if (to_axis > ds.ndim()) to_axis = ds.ndim();
+      assert (from_axis < to_axis);
       off64_t fp = 1;
-      for (size_t n = 0; n < up_to_dim; n++) 
+      for (size_t n = from_axis; n < to_axis; n++) 
         fp *= ds.dim(n);
       return (fp); 
     }
@@ -57,9 +58,29 @@ namespace MR {
       template <typename T> inline bool is_complex__ () { return (false); }
       template <> inline bool is_complex__<cfloat> () { return (true); }
       template <> inline bool is_complex__<cdouble> () { return (true); }
+      template <class Set> class CompareStride {
+        public:
+          CompareStride (const Set& set) : S (set) { }
+          bool operator() (const size_t a, const size_t b) const { return (abs(S.stride(a)) < abs(S.stride(b))); }
+        private:
+          const Set& S;
+      };
     }
     //! \endcond
 
+
+
+    //! sort strides in increasing absolute order and return vector of indices
+    template <class Set> std::vector<size_t> stride_order (const Set& set) {
+      std::vector<size_t> ret (set.ndim());
+      for (size_t i = 0; i < ret.size(); ++i) ret[i] = i;
+      CompareStride<Set> compare (set);
+      std::sort (ret.begin(), ret.end(), compare);
+      return (ret);
+    }
+
+
+
     //! return whether the Set contains complex data
     template <class Set> inline bool is_complex (const Set& ds) { 
       typedef typename Set::value_type T;
diff --git a/lib/dataset/position.h b/lib/dataset/position.h
new file mode 100644
index 0000000000..56fe6839b4
--- /dev/null
+++ b/lib/dataset/position.h
@@ -0,0 +1,54 @@
+/*
+    Copyright 2009 Brain Research Institute, Melbourne, Australia
+
+    Written by J-Donald Tournier, 07/01/10.
+
+    This file is part of MRtrix.
+
+    MRtrix is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    MRtrix is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with MRtrix.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+#ifndef __dataset_position_h__
+#define __dataset_position_h__
+
+namespace MR {
+  namespace DataSet {
+
+    //! \addtogroup DataSet
+    // @{
+
+    template <class Set> class Position 
+    {
+      public:
+        Position (Set& parent, size_t corresponding_axis) : S (parent), axis (corresponding_axis) { }
+        operator ssize_t () const          { return (S.get_pos (axis)); }
+        ssize_t operator++ ()              { ssize_t p = S.get_pos(axis); S.move_pos (axis,1); return (p); }
+        ssize_t operator-- ()              { ssize_t p = S.get_pos(axis); S.move_pos (axis,-1); return (p); }
+        ssize_t operator++ (int notused)   { S.move_pos (axis,1); return (S.get_pos (axis)); }
+        ssize_t operator-- (int notused)   { S.move_pos (axis,-1); return (S.get_pos (axis)); }
+        ssize_t operator+= (ssize_t increment) { S.move_pos (axis, increment); return (S.get_pos(axis)); }
+        ssize_t operator-= (ssize_t increment) { S.move_pos (axis, -increment); return (S.get_pos(axis)); }
+        ssize_t operator= (ssize_t position)   { S.set_pos (axis, position); return (position); }
+      private:
+        Set& S;
+        size_t axis;
+    };
+
+    //! @}
+  }
+}
+
+#endif
+
diff --git a/lib/dataset/reorder.h b/lib/dataset/reorder.h
index 340d6d8463..ab34734d93 100644
--- a/lib/dataset/reorder.h
+++ b/lib/dataset/reorder.h
@@ -24,6 +24,9 @@
 #define __dataset_reorder_h__
 
 #include "math/matrix.h"
+#include "dataset/misc.h"
+#include "dataset/value.h"
+#include "dataset/position.h"
 
 namespace MR {
   namespace DataSet {
@@ -35,35 +38,49 @@ namespace MR {
       public:
         typedef typename Set::value_type value_type;
 
-        Reorder (Set& original, const size_t* ordering = NULL, const std::string& description = "") : 
+        Reorder (Set& original, const std::string& description = "") : 
           D (original),
-          order (new size_t [ndim()]),
-          descriptor (description.empty() ? D.name() + " [reordered]" : description) {
-            memcpy (order, ordering ? ordering : D.layout(), ndim()*sizeof(size_t));
-        }
+          S (stride_order (original)),
+          descriptor (description.empty() ? D.name() + " [reordered]" : description) { }
+
+        template <class Set2> Reorder (Set& original, const Set2& reference, const std::string& description = "") : 
+          D (original),
+          S (stride_order (reference)),
+          descriptor (description.empty() ? D.name() + " [reordered]" : description) { }
+
+        Reorder (Set& original, const std::vector<size_t>& ordering, const std::string& description = "") : 
+          D (original),
+          S (ordering),
+          descriptor (description.empty() ? D.name() + " [reordered]" : description) { }
 
         const std::string& name () const { return (descriptor); }
         size_t  ndim () const { return (D.ndim()); }
-        int     dim (size_t axis) const { return (D.dim(order[axis])); }
-        const size_t* layout () const { return (order); }
+        int     dim (size_t axis) const { return (D.dim(S[axis])); }
+        ssize_t stride (size_t axis) const { return (D.stride (S[axis])); }
+        const std::vector<size_t>& order () const { return (S); }
 
-        float   vox (size_t axis) const { return (D.vox(order[axis])); }
+        float   vox (size_t axis) const { return (D.vox(S[axis])); }
 
         const Math::Matrix<float>& transform () const { return (D.transform()); }
 
         void    reset () { D.reset(); }
 
-        ssize_t pos (size_t axis) const { return (D.pos(order[axis])); }
-        void    pos (size_t axis, ssize_t position) const { D.pos(order[axis], position); }
-        void    move (size_t axis, ssize_t increment) const { D.move (order[axis], increment); }
-
-        value_type   value () const { return (D.value()); }
-        void         value (value_type val) { D.value (val); }
+        Position<Reorder<Set> > operator[] (size_t axis) { return (Position<Reorder<Set> > (*this, axis)); }
+        Value<Reorder<Set> > value () { return (Value<Reorder<Set> > (*this)); }
 
       private:
         Set& D;
-        size_t* order;
+        std::vector<size_t> S;
         std::string descriptor;
+
+        ssize_t get_pos (size_t axis) const { return (D[S[axis]]); }
+        void    set_pos (size_t axis, ssize_t position) const { D[S[axis]] = position; }
+        void    move_pos (size_t axis, ssize_t increment) const { D[S[axis]] += increment; }
+        value_type get_value () const { return (D.value()); }
+        void set_value (value_type val) { D.value() = val; }
+
+        friend class Position<Reorder<Set> >;
+        friend class Value<Reorder<Set> >;
     };
 
     //! @}
diff --git a/lib/dataset/subset.h b/lib/dataset/subset.h
index d55217f18d..d6c55b1c58 100644
--- a/lib/dataset/subset.h
+++ b/lib/dataset/subset.h
@@ -24,6 +24,8 @@
 #define __dataset_subset_h__
 
 #include "math/matrix.h"
+#include "dataset/value.h"
+#include "dataset/position.h"
 
 namespace MR {
   namespace DataSet {
@@ -60,12 +62,8 @@ namespace MR {
 
         void    reset () { for (size_t n = 0; n < NDIM; ++n) pos(n,0); }
 
-        ssize_t pos (size_t axis) const { return (D.pos(axis)-C[axis]); }
-        void    pos (size_t axis, ssize_t position) const { D.pos(axis, position + C[axis]); }
-        void    move (size_t axis, ssize_t increment) const { D.move (axis, increment); }
-
-        value_type   value () const { return (D.value()); }
-        void         value (value_type val) { D.value (val); }
+        Value<Subset<Set> > value () { return (Value<Subset<Set> > (*this)); }
+        Position<Subset<Set> > operator[] (size_t axis) { return (Position<Subset<Set> > (*this, axis)); }
 
       private:
         Set& D;
@@ -73,6 +71,15 @@ namespace MR {
         size_t N [NDIM];
         std::string descriptor;
         Math::Matrix<float> transform_matrix;
+
+        value_type get_value () const { return (D.value()); }
+        void set_value (value_type val) { D.value() = val; }
+        ssize_t get_pos (size_t axis) const { return (D[axis]-C[axis]); }
+        void set_pos (size_t axis, ssize_t position) const { D[axis] = position + C[axis]; }
+        void move_pos (size_t axis, ssize_t increment) const { D[axis] += increment; }
+
+        friend class Value<Subset<Set> >;
+        friend class Position<Subset<Set> >;
     };
 
     //! @}
diff --git a/lib/dataset/transform.h b/lib/dataset/transform.h
index 4753f4f716..adfd7b0d80 100644
--- a/lib/dataset/transform.h
+++ b/lib/dataset/transform.h
@@ -34,9 +34,9 @@ namespace MR {
 
     namespace Transform {
 
-      template <class Set, typename T> inline Math::MatrixView<T>& P2I (Math::MatrixView<T>& M, const Set& ds) 
+      template <class Set, typename T> inline Math::Matrix<T>& P2I (Math::Matrix<T>& M, const Set& ds) 
       {
-        assert (M.rows() == 4 && M.columns() == 4);
+        M.allocate(4,4);
         M.zero();
         M(0,0) = ds.vox(0);
         M(1,1) = ds.vox(1);
@@ -44,12 +44,10 @@ namespace MR {
         M(3,3) = 1.0;
         return (M);
       }
-      template <class Set, typename T> inline Math::Matrix<T>& P2I (Math::Matrix<T>& M, const Set& ds) 
-      { M.allocate(4,4); P2I(M.view(), ds); return (M); }
 
 
-      template <class Set, typename T> inline Math::MatrixView<T>& I2P (Math::MatrixView<T>& M, const Set& ds) {
-        assert (M.rows() == 4 && M.columns() == 4);
+      template <class Set, typename T> inline Math::Matrix<T>& I2P (Math::Matrix<T>& M, const Set& ds) {
+        M.allocate(4,4);
         M.zero();
         M(0,0) = 1.0/ds.vox(0);
         M(1,1) = 1.0/ds.vox(1);
@@ -57,23 +55,17 @@ namespace MR {
         M(3,3) = 1.0;
         return (M);
       }
-      template <class Set, typename T> inline Math::Matrix<T>& I2P (Math::Matrix<T>& M, const Set& ds) 
-      { M.allocate(4,4); I2P(M.view(), ds); return (M); }
 
 
-      template <class Set, typename T> inline Math::MatrixView<T>& R2I (Math::MatrixView<T>& M, const Set& ds) {
-        assert (M.rows() == 4 && M.columns() == 4);
-        assert (ds.transform().rows() == 4 && ds.transform().columns() == 4);
+      template <class Set, typename T> inline Math::Matrix<T>& R2I (Math::Matrix<T>& M, const Set& ds) {
+        M.allocate(4,4);
         M = ds.transform();
         return (M);
       }
-      template <class Set, typename T> inline Math::Matrix<T>& R2I (Math::Matrix<T>& M, const Set& ds) 
-      { M.allocate(4,4); R2I(M.view(), ds); return (M); }
 
 
-      template <class Set, typename T> inline Math::MatrixView<T>& I2R (Math::MatrixView<T>& M, const Set& ds) {
-        assert (M.rows() == 4 && M.columns() == 4);
-        assert (ds.transform().rows() == 4 && ds.transform().columns() == 4);
+      template <class Set, typename T> inline Math::Matrix<T>& I2R (Math::Matrix<T>& M, const Set& ds) {
+        M.allocate(4,4);
         int signum;
         Math::Permutation p (4);
         Math::Matrix<T> D (ds.transform());
@@ -82,26 +74,20 @@ namespace MR {
         M(3,0) = M(3,1) = M(3,2) = 0.0; M(3,3) = 1.0;
         return (M);
       }
-      template <class Set, typename T> inline Math::Matrix<T>& I2R (Math::Matrix<T>& M, const Set& ds) 
-      { M.allocate(4,4); I2R(M.view(), ds); return (M); }
 
 
-      template <class Set, typename T> inline Math::MatrixView<T>& P2R (Math::MatrixView<T>& M, const Set& ds) {
-        assert (M.rows() == 4 && M.columns() == 4);
-        assert (ds.transform().rows() == 4 && ds.transform().columns() == 4);
+      template <class Set, typename T> inline Math::Matrix<T>& P2R (Math::Matrix<T>& M, const Set& ds) {
+        M.allocate(4,4);
         M = ds.transform();
         for (size_t i = 0; i < 3; i++) 
           for (size_t j = 0; j < 3; j++) 
             M(i,j) *= ds.vox(i);
         return (M);
       }
-      template <class Set, typename T> inline Math::Matrix<T>& P2R (Math::Matrix<T>& M, const Set& ds) 
-      { M.allocate(4,4); P2R(M.view(), ds); return (M); }
 
 
-      template <class Set, typename T> inline Math::MatrixView<T>& R2P (Math::MatrixView<T>& M, const Set& ds) {
-        assert (M.rows() == 4 && M.columns() == 4);
-        assert (ds.transform().rows() == 4 && ds.transform().columns() == 4);
+      template <class Set, typename T> inline Math::Matrix<T>& R2P (Math::Matrix<T>& M, const Set& ds) {
+        M.allocate(4,4);
         int signum;
         Math::Permutation p (4);
         Math::Matrix<T> D (4,4);
@@ -111,8 +97,6 @@ namespace MR {
         M(3,0) = M(3,1) = M(3,2) = 0.0; M(3,3) = 1.0;
         return (M);
       }
-      template <class Set, typename T> inline Math::Matrix<T>& R2P (Math::Matrix<T>& M, const Set& ds) 
-      { M.allocate(4,4); R2P(M.view(), ds); return (M); }
 
     }
     //! @}
diff --git a/lib/dataset/value.h b/lib/dataset/value.h
new file mode 100644
index 0000000000..194b24489f
--- /dev/null
+++ b/lib/dataset/value.h
@@ -0,0 +1,53 @@
+/*
+    Copyright 2009 Brain Research Institute, Melbourne, Australia
+
+    Written by J-Donald Tournier, 07/01/10.
+
+    This file is part of MRtrix.
+
+    MRtrix is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    MRtrix is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with MRtrix.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+#ifndef __dataset_value_h__
+#define __dataset_value_h__
+
+namespace MR {
+  namespace DataSet {
+
+    //! \addtogroup DataSet
+    // @{
+
+    template <class Set> class Value {
+      public:
+        typedef typename Set::value_type value_type;
+
+        Value (Set& parent) : S (parent) { }
+        operator value_type () const { return (S.get_value()); }
+        value_type operator= (value_type value) { S.set_value (value); return (value); }
+        value_type operator= (const Value& V) { value_type value = V.S.get_value(); S.set_value (value); return (value); }
+        value_type operator+= (value_type value) { value += S.get_value(); S.set_value (value); return (value); }
+        value_type operator-= (value_type value) { value = S.get_value() - value; S.set_value (value); return (value); }
+        value_type operator*= (value_type value) { value *= S.get_value(); S.set_value (value); return (value); }
+        value_type operator/= (value_type value) { value = S.get_value() / value; S.set_value (value); return (value); }
+      private:
+        Set& S;
+    };
+
+    //! @}
+  }
+}
+
+#endif
+
diff --git a/lib/documentation.h b/lib/documentation.h
index d240fc324a..1230971afc 100644
--- a/lib/documentation.h
+++ b/lib/documentation.h
@@ -438,18 +438,13 @@ namespace MR {
   *
   *     int     ndim () const           { return (3); }
   *     int     dim (int axis) const    { return (N[axis]); }
-  *     int     pos (int axis)          { return (X[axis]); }
-  *     void    pos (int axis, int i)   { X[axis] = i; }
-  *     void    move (int axis, int i)  { X[axis] += i; }
-  *     float   value ()                { return (data[offset()]); }
-  *     void    value (float val)       { data[offset()] = val; }
+  *     int&    operator[] (int axis)   { return (X[axis]); }
+  *     float&  value ()                { return (data[X[0]+N[0]*(X[1]+N[1]*X[2])]); }
   *
   *   private:
   *     float*  data
   *     int     N[3];
   *     int     X[3];
-  *
-  *     int     offset () const         { return (X[0]+N[0]*(X[1]+N[1]*X[2])); }
   * };
   * \endcode
   *
@@ -462,10 +457,10 @@ namespace MR {
   * \code
   * template <class Set> void scale (Set& data, float factor)
   * {
-  *   for (data.pos(2,0); data.pos(2) < data.dim(2); data.move(2,1))
-  *     for (data.pos(1,0); data.pos(1) < data.dim(1); data.move(1,1))
-  *       for (data.pos(0,0); data.pos(0) < data.dim(0); data.move(0,1))
-  *         data.value (factor * data.value());
+  *   for (data[2] = 0; data[2] < data.dim(2); ++data[2])
+  *     for (data[1] = 0; data[1] < data.dim(1); ++data[1])
+  *       for (data[0] = 0; data[0] < data.dim(0); ++data[0])
+  *         data.value() *= factor;
   * }
   * \endcode
   *
@@ -483,7 +478,7 @@ namespace MR {
   * required for the \a scale() function to compile and run. There is also
   * plenty of scope for the compiler to optimise this particular function,
   * since all member functions of \c my_image are declared inline. Note that
-  * this does not mean that this class can be used with any of the other
+  * this does not mean that this class can be used with all of the other
   * template functions, some of which might rely on some of the other member
   * functions having been defined.
   *
@@ -514,6 +509,9 @@ namespace MR {
   * interface were required, it would be trivial to define such an abstract class and
   * use it with the template functions provided by MRtrix.
   *
+  * \sa the DataSet::Position and DataSet::Value template classes are designed
+  * to simplify the process of returning a modifiable object for non-trivial
+  * DataSet classes.
   */
 
 
diff --git a/lib/file/nifti1_utils.cpp b/lib/file/nifti1_utils.cpp
index 3a606c470f..b7ae5e1928 100644
--- a/lib/file/nifti1_utils.cpp
+++ b/lib/file/nifti1_utils.cpp
@@ -305,11 +305,11 @@ namespace MR {
         if (permutation[0] != 0 || permutation[1] != 1 || permutation[2] != 2 || 
             !H.axes.forward(0)  || !H.axes.forward(1)  || !H.axes.forward(2)) {
 
-          Math::VectorView<float> translation = M.column(3).view(0,3);
+          Math::Vector<float> translation = M.column(3).sub(0,3);
           for (size_t i = 0; i < 3; i++) {
             if (!H.axes.forward(i)) {
               float length = float(H.axes.dim(i)-1) * H.axes.vox(i);
-              Math::VectorView<float> axis = M.column(i).view(0,3);
+              Math::Vector<float> axis = M.column(i).sub(0,3);
               for (size_t n = 0; n < 3; n++) {
                 axis[n] = -axis[n];
                 translation[n] -= length*axis[n];
@@ -318,7 +318,7 @@ namespace MR {
           }
 
           for (size_t i = 0; i < 3; i++) {
-            Math::VectorView<float> row = M.row(i).view(0,3);
+            Math::Vector<float> row = M.row(i).sub(0,3);
             permutation.apply (row); 
           }
         }
diff --git a/lib/image/axis.cpp b/lib/image/axis.cpp
index 5770308102..dddbad25a5 100644
--- a/lib/image/axis.cpp
+++ b/lib/image/axis.cpp
@@ -93,7 +93,7 @@ namespace MR {
 
 
 
-    std::vector<ssize_t> parse_axes_specifier (size_t ndim, const std::string& specifier)
+    std::vector<ssize_t> Axes::parse (size_t ndim, const std::string& specifier)
     {
       std::vector<ssize_t> parsed (ndim);
 
@@ -125,7 +125,7 @@ namespace MR {
       if (current != ndim) 
         throw Exception ("incorrect number of axes in axes specification \"" + specifier + "\"");
 
-      check_axes_specifier (parsed, ndim);
+      check (parsed, ndim);
 
       return (parsed);
     }
@@ -136,7 +136,7 @@ namespace MR {
 
 
 
-    void check_axes_specifier (const std::vector<ssize_t>& parsed, size_t ndim)
+    void Axes::check (const std::vector<ssize_t>& parsed, size_t ndim)
     {
       if (parsed.size() != ndim) 
         throw Exception ("incorrect number of dimensions for axes specifier");
@@ -158,6 +158,7 @@ namespace MR {
 
 
 
+
     std::ostream& operator<< (std::ostream& stream, const Axes& axes)
     {
       stream << "dim [ ";
diff --git a/lib/image/axis.h b/lib/image/axis.h
index 6de1d8c0ae..629cd3092e 100644
--- a/lib/image/axis.h
+++ b/lib/image/axis.h
@@ -29,16 +29,6 @@
 namespace MR {
   namespace Image {
 
-    /*
-        class CompareStride {
-          public:
-            CompareStride (const Axes& parent) : A (parent) { }
-            bool operator() (const size_t a, const size_t b) const { return (abs(A.stride(a)) < abs (A.stride(b))); }
-          private:
-            const Axes& A;
-        };
-*/
-
     //! \addtogroup Image 
     // @{
     
@@ -113,6 +103,10 @@ namespace MR {
         bool  forward (size_t index) const { return (axes[index].stride > 0); }
         ssize_t direction (size_t index) const { return (axes[index].stride > 0 ? 1 : -1); }
 
+        friend std::ostream& operator<< (std::ostream& stream, const Axes& axes);
+        static std::vector<ssize_t> parse (size_t ndim, const std::string& specifier);
+        static void check (const std::vector<ssize_t>& parsed, size_t ndim);
+
       protected:
         std::vector<Axis> axes;
         static const std::string axes_name;
@@ -143,10 +137,6 @@ namespace MR {
     
     //! @}
 
-    std::ostream& operator<< (std::ostream& stream, const Axes& axes);
-    std::vector<ssize_t> parse_axes_specifier (size_t ndim, const std::string& specifier);
-    void check_axes_specifier (const std::vector<ssize_t>& parsed, size_t ndim);
-
   }
 }
 
diff --git a/lib/image/format/mrtrix.cpp b/lib/image/format/mrtrix.cpp
index 23908fb750..6119beb1e9 100644
--- a/lib/image/format/mrtrix.cpp
+++ b/lib/image/format/mrtrix.cpp
@@ -93,7 +93,7 @@ namespace MR {
 
 
         if (layout.empty()) throw Exception ("missing \"layout\" specification for MRtrix image \"" + H.name() + "\"");
-        std::vector<ssize_t> ax = parse_axes_specifier (H.ndim(), layout);
+        std::vector<ssize_t> ax = Axes::parse (H.ndim(), layout);
         for (size_t i = 0; i < ax.size(); i++) H.axes.stride(i) = ax[i];
 
         for (size_t n = 0; n < MIN (H.axes.ndim(), labels.size()); n++) H.axes.description(n) = labels[n];
@@ -107,7 +107,7 @@ namespace MR {
             for (int col = 0; col < 4; col++) 
               T(row,col) = transform[count++];
           T(3,0) = T(3,1) = T(3,2) = 0.0; T(3,3) = 1.0;
-          H.transform() = T;
+          H.transform().copy (T);
         }
 
 
@@ -119,7 +119,7 @@ namespace MR {
             for (size_t row = 0; row < M.rows(); row++) 
               for (size_t col = 0; col < 4; col++) 
                 M(row,col) = dw_scheme[count++];
-            H.DW_scheme = M;
+            H.DW_scheme.copy (M);
           }
         }
 
diff --git a/lib/image/header.cpp b/lib/image/header.cpp
index a8592bd192..f34575a11a 100644
--- a/lib/image/header.cpp
+++ b/lib/image/header.cpp
@@ -77,9 +77,9 @@ namespace MR {
       transform_matrix(3,0) = transform_matrix(3,1) = transform_matrix(3,2) = 0.0; transform_matrix(3,3) = 1.0;
 
       Math::Permutation permutation (3);
-      Math::absmax (transform_matrix.row(0).view(0,3), permutation[0]);
-      Math::absmax (transform_matrix.row(1).view(0,3), permutation[1]);
-      Math::absmax (transform_matrix.row(2).view(0,3), permutation[2]);
+      Math::absmax (transform_matrix.row(0).sub(0,3), permutation[0]);
+      Math::absmax (transform_matrix.row(1).sub(0,3), permutation[1]);
+      Math::absmax (transform_matrix.row(2).sub(0,3), permutation[2]);
 
       assert (permutation[0] != permutation[1] && permutation[1] != permutation[2] && permutation[2] != permutation[0]);
 
@@ -97,16 +97,16 @@ namespace MR {
         axes[2] = a[2];
 
         for (size_t i = 0; i < 3; i++) {
-          Math::VectorView<float> row = transform_matrix.row(i).view(0,3);
+          Math::Vector<float> row = transform_matrix.row(i).sub(0,3);
           permutation.apply (row); 
         }
 
-        Math::VectorView<float> translation = transform_matrix.column(3).view(0,3);
+        Math::Vector<float> translation = transform_matrix.column(3).sub(0,3);
         for (size_t i = 0; i < 3; i++) {
           if (flip[i]) {
             axes.stride(i) = -axes.stride(i);
             float length = float(axes[i].dim-1) * axes[i].vox;
-            Math::VectorView<float> axis = transform_matrix.column(i).view(0,3);
+            Math::Vector<float> axis = transform_matrix.column(i).sub(0,3);
             for (size_t n = 0; n < 3; n++) {
               axis[n] = -axis[n];
               translation[n] -= length*axis[n];
diff --git a/lib/image/header.h b/lib/image/header.h
index 629ff15daf..6c246166b6 100644
--- a/lib/image/header.h
+++ b/lib/image/header.h
@@ -118,7 +118,8 @@ namespace MR {
         template <typename T> float scale_to_storage (T val) const   { return (scale_to_storage (val, scale, offset)); }
 
         //! returns the memory footprint of the Image
-        off64_t footprint (size_t up_to_dim = SIZE_MAX) { return (footprint_for_count (DataSet::voxel_count (*this, up_to_dim))); }
+        off64_t footprint (size_t from_dim = 0, size_t up_to_dim = SIZE_MAX) { 
+          return (footprint_for_count (DataSet::voxel_count (*this, from_dim, up_to_dim))); }
 
         //! returns the memory footprint of a DataSet
         off64_t footprint (const char* specifier) { return (footprint_for_count (DataSet::voxel_count (*this, specifier))); }
@@ -132,7 +133,6 @@ namespace MR {
         std::string          identifier;
         DataType             dtype;
         Math::Matrix<float>  transform_matrix;
-        mutable Array<size_t>::Ptr axes_layout; 
 
         void merge (const Header& H);
 
diff --git a/lib/image/voxel.h b/lib/image/voxel.h
index 2fd31593bb..138b2fb963 100644
--- a/lib/image/voxel.h
+++ b/lib/image/voxel.h
@@ -26,6 +26,8 @@
 #include "get_set.h"
 #include "image/header.h"
 #include "math/complex.h"
+#include "dataset/value.h"
+#include "dataset/position.h"
 
 #define MAX_FILES_PER_IMAGE 256U
 #define MAX_NDIM 16
@@ -47,6 +49,8 @@ namespace MR {
       template <typename value_type, typename S> void __putBE (value_type val, void* data, size_t i) { return (MR::putBE<S> (val, data, i)); }
 
       // specialisation for conversion to bool
+      template <> bool __getLE<bool,float> (const void* data, size_t i) { return (Math::round(MR::getLE<float> (data, i))); }
+      template <> bool __getBE<bool,float> (const void* data, size_t i) { return (Math::round(MR::getBE<float> (data, i))); }
       template <> bool __getLE<bool,double> (const void* data, size_t i) { return (Math::round(MR::getLE<double> (data, i))); }
       template <> bool __getBE<bool,double> (const void* data, size_t i) { return (Math::round(MR::getBE<double> (data, i))); }
     }
@@ -84,21 +88,66 @@ namespace MR {
           offset = handler.start();
 
           switch (H.datatype()()) {
-            case DataType::Bit:        get_func = &__get<value_type,bool>;       put_func = &__put<value_type,bool>;       return;
-            case DataType::Int8:       get_func = &__get<value_type,int8_t>;     put_func = &__put<value_type,int8_t>;     return;
-            case DataType::UInt8:      get_func = &__get<value_type,uint8_t>;    put_func = &__put<value_type,uint8_t>;    return;
-            case DataType::Int16LE:    get_func = &__getLE<value_type,int16_t>;  put_func = &__putLE<value_type,int16_t>;  return;
-            case DataType::UInt16LE:   get_func = &__getLE<value_type,uint16_t>; put_func = &__putLE<value_type,uint16_t>; return;
-            case DataType::Int16BE:    get_func = &__getBE<value_type,int16_t>;  put_func = &__putBE<value_type,int16_t>;  return;
-            case DataType::UInt16BE:   get_func = &__getBE<value_type,uint16_t>; put_func = &__putBE<value_type,uint16_t>; return;
-            case DataType::Int32LE:    get_func = &__getLE<value_type,int32_t>;  put_func = &__putLE<value_type,int32_t>;  return;
-            case DataType::UInt32LE:   get_func = &__getLE<value_type,uint32_t>; put_func = &__putLE<value_type,uint32_t>; return;
-            case DataType::Int32BE:    get_func = &__getBE<value_type,int32_t>;  put_func = &__putBE<value_type,int32_t>;  return;
-            case DataType::UInt32BE:   get_func = &__getBE<value_type,uint32_t>; put_func = &__putBE<value_type,uint32_t>; return;
-            case DataType::Float32LE:  get_func = &__getLE<value_type,float>;    put_func = &__putLE<value_type,float>;    return;
-            case DataType::Float32BE:  get_func = &__getBE<value_type,float>;    put_func = &__putBE<value_type,float>;    return;
-            case DataType::Float64LE:  get_func = &__getLE<value_type,double>;   put_func = &__putLE<value_type,double>;   return;
-            case DataType::Float64BE:  get_func = &__getBE<value_type,double>;   put_func = &__putBE<value_type,double>;   return;
+            case DataType::Bit:
+              get_func = &__get<value_type,bool>;   
+              put_func = &__put<value_type,bool>;
+              return;
+            case DataType::Int8:
+              get_func = &__get<value_type,int8_t>;
+              put_func = &__put<value_type,int8_t>;
+              return;
+            case DataType::UInt8:
+              get_func = &__get<value_type,uint8_t>;
+              put_func = &__put<value_type,uint8_t>;
+              return;
+            case DataType::Int16LE:
+              get_func = &__getLE<value_type,int16_t>;
+              put_func = &__putLE<value_type,int16_t>;
+              return;
+            case DataType::UInt16LE:
+              get_func = &__getLE<value_type,uint16_t>;
+              put_func = &__putLE<value_type,uint16_t>;
+              return;
+            case DataType::Int16BE:
+              get_func = &__getBE<value_type,int16_t>;
+              put_func = &__putBE<value_type,int16_t>;
+              return;
+            case DataType::UInt16BE:
+              get_func = &__getBE<value_type,uint16_t>;
+              put_func = &__putBE<value_type,uint16_t>;
+              return;
+            case DataType::Int32LE:
+              get_func = &__getLE<value_type,int32_t>;
+              put_func = &__putLE<value_type,int32_t>;
+              return;
+            case DataType::UInt32LE:
+              get_func = &__getLE<value_type,uint32_t>;
+              put_func = &__putLE<value_type,uint32_t>;
+              return;
+            case DataType::Int32BE:
+              get_func = &__getBE<value_type,int32_t>;
+              put_func = &__putBE<value_type,int32_t>;
+              return;
+            case DataType::UInt32BE:
+              get_func = &__getBE<value_type,uint32_t>;
+              put_func = &__putBE<value_type,uint32_t>;
+              return;
+            case DataType::Float32LE:
+              get_func = &__getLE<value_type,float>;
+              put_func = &__putLE<value_type,float>;
+              return;
+            case DataType::Float32BE:
+              get_func = &__getBE<value_type,float>;
+              put_func = &__putBE<value_type,float>;
+              return;
+            case DataType::Float64LE:
+              get_func = &__getLE<value_type,double>;
+              put_func = &__putLE<value_type,double>;
+              return;
+            case DataType::Float64BE:
+              get_func = &__getBE<value_type,double>;
+              put_func = &__putBE<value_type,double>;
+              return;
             default: throw Exception ("invalid data type in image header");
           }
         }
@@ -117,7 +166,7 @@ namespace MR {
           return (false);
         }
 
-        const size_t* layout () const { return (H.layout()); }
+        ssize_t stride (size_t axis) const { return (handler.stride (axis)); }
         size_t  ndim () const { return (H.ndim()); }
         ssize_t dim (size_t axis) const { return (H.dim(axis)); }
         float   vox (size_t axis) const { return (H.vox(axis)); }
@@ -127,7 +176,7 @@ namespace MR {
           ssize_t shift = 0;
           for (size_t n = 0; n < ndim(); n++) {
             x[n] = V.pos(n);
-            shift += handler.stride(n) * x[n];
+            shift += stride(n) * x[n];
           }
           offset = handler.start() + shift;
           return (*this);
@@ -135,23 +184,13 @@ namespace MR {
 
         //! reset all coordinates to zero. 
         void reset () { offset = handler.start(); for (size_t i = 0; i < ndim(); i++) x[i] = 0; }
-
-        ssize_t pos (size_t axis) const { return (x[axis]); }
-        void    pos (size_t axis, ssize_t position) { offset += handler.stride(axis) * (position - x[axis]); x[axis] = position; }
-        void    move (size_t axis, ssize_t increment) { offset += handler.stride(axis) * increment; x[axis] += increment; }
-
-        value_type value () const { 
-          ssize_t nseg (offset / handler.segment_size());
-          return (H.scale_from_storage (get_func (handler.segment(nseg), offset - nseg*handler.segment_size()))); 
-        }
-        void value (value_type val) {
-          ssize_t nseg (offset / handler.segment_size());
-          put_func (H.scale_to_storage (val), handler.segment(nseg), offset - nseg*handler.segment_size()); 
-        }
+       
+        DataSet::Position<Voxel<T> > operator[] (size_t axis) { return (DataSet::Position<Voxel<T> > (*this, axis)); }
+        DataSet::Value<Voxel<T> > value () { return (DataSet::Value<Voxel<T> > (*this)); }
 
         friend std::ostream& operator<< (std::ostream& stream, const Voxel& V) {
           stream << "position for image \"" << V.name() << "\" = [ ";
-          for (size_t n = 0; n < V.ndim(); n++) stream << V[n] << " ";
+          for (size_t n = 0; n < V.ndim(); ++n) stream << const_cast<Voxel&>(V)[n] << " ";
           stream << "]\n  current offset = " << V.offset;
           return (stream);
         }
@@ -164,6 +203,27 @@ namespace MR {
 
         value_type (*get_func) (const void* data, size_t i);
         void       (*put_func) (value_type val, void* data, size_t i);
+
+        value_type get_value () const { 
+          ssize_t nseg (offset / handler.segment_size());
+          return (H.scale_from_storage (get_func (handler.segment(nseg), offset - nseg*handler.segment_size()))); 
+        }
+        void set_value (value_type val) {
+          ssize_t nseg (offset / handler.segment_size());
+          put_func (H.scale_to_storage (val), handler.segment(nseg), offset - nseg*handler.segment_size()); 
+        }
+        ssize_t get_pos (size_t axis) const { return (x[axis]); }
+        void set_pos (size_t axis, ssize_t position) { 
+          offset += stride(axis) * (position - x[axis]); 
+          x[axis] = position;
+        }
+        void move_pos (size_t axis, ssize_t increment) {
+          offset += stride(axis) * increment;
+          x[axis] += increment;
+        }
+        
+        friend class DataSet::Position<Voxel<T> >;
+        friend class DataSet::Value<Voxel<T> >;
     };
 
 
diff --git a/lib/math/LU.h b/lib/math/LU.h
index 98e190c96c..dcd0fcc2d8 100644
--- a/lib/math/LU.h
+++ b/lib/math/LU.h
@@ -55,45 +55,35 @@ namespace MR {
 
       //! %LU decomposition of A
       /** \note the contents of \a A will be overwritten with its %LU decomposition */
-      template <typename T> inline MatrixView<T>& decomp (MatrixView<T>& A, Permutation& p, int& signum) { 
+      template <typename T> inline Matrix<T>& decomp (Matrix<T>& A, Permutation& p, int& signum) { 
         gsl_linalg_LU_decomp (&A, &p, &signum);
         return (A); 
       }
 
       //! inverse of A given its %LU decomposition D,p
-      template <typename T> inline MatrixView<T>& inv (MatrixView<T>& I, const MatrixView<T>D, const Permutation& p) { 
+      template <typename T> inline Matrix<T>& inv (Matrix<T>& I, const Matrix<T>D, const Permutation& p) { 
+	I.allocate (D); 
         gsl_linalg_LU_invert (&D, &p, &I);
         return (I); 
       }
 
       //! solve A*x = b given %LU decomposition D,p of A
-      template <typename T> inline VectorView<T>& solve (VectorView<T>& x, const MatrixView<T>& D, const Permutation& p, const VectorView<T>& b) {
+      template <typename T> inline Vector<T>& solve (Vector<T>& x, const Matrix<T>& D, const Permutation& p, const Vector<T>& b) {
+        x.allocate (D.rows());
         gsl_linalg_LU_solve (&D, &p, &b, &x);
         return (x);
       }
 
       //! solve A*x = b given %LU decomposition D,p of A, in place (b passed in as x).
-      template <typename T> inline VectorView<T>& solve (VectorView<T>& x, const MatrixView<T>& D, const Permutation& p)
+      template <typename T> inline Vector<T>& solve (Vector<T>& x, const Matrix<T>& D, const Permutation& p)
       {
         gsl_linalg_LU_svx (&D, &p, &x);
         return (x);
       }
 
-      //! solve A*x = b given %LU decomposition D,p of A
-      template <typename T> inline Vector<T>& solve (Vector<T>& x, const MatrixView<T>& D, const Permutation& p, const VectorView<T>& b) {
-        x.allocate (D.rows());
-        solve (x.view(), D, p, b);
-        return (x);
-      }
-
-      //! inverse of A given its %LU decomposition D,p
-      template <typename T> inline MatrixView<T>& inv (Matrix<T>& I, const MatrixView<T>D, const Permutation& p) {
-	I.allocate (D); 
-	return (inv (I.view(), D, p)); 
-      }
-
       //! inverse of A by %LU decomposition
-      template <typename T> inline MatrixView<T>& inv (MatrixView<T>& I, const MatrixView<T>& A) { 
+      template <typename T> inline Matrix<T>& inv (Matrix<T>& I, const Matrix<T>& A) { 
+        I.allocate (A);
 	Permutation p (A.rows());
         int signum;
         Matrix<T> D (A);
@@ -101,12 +91,6 @@ namespace MR {
         return (inv (I, D, p));
       }
 
-      //! inverse of A by %LU decomposition
-      template <typename T> inline Matrix<T>& inv (Matrix<T>& I, const MatrixView<T>& A) { 
-        I.allocate (A);
-        inv (I.view(), A);
-        return (I); 
-      }
 
       /** @} */
       /** @} */
diff --git a/lib/math/SH.h b/lib/math/SH.h
index 0743211261..ed3d857602 100644
--- a/lib/math/SH.h
+++ b/lib/math/SH.h
@@ -112,8 +112,9 @@ namespace MR {
       }
 
 
-      template <typename T> inline Math::VectorView<T>& delta (Math::VectorView<T>& D, const Point& unit_dir, int lmax)
+      template <typename T> inline Math::Vector<T>& delta (Math::Vector<T>& D, const Point& unit_dir, int lmax)
       {
+        D.allocate (NforL (lmax));
         T az = Math::atan2 (unit_dir[1], unit_dir[0]);
         T AL [lmax+1];
         Legendre::Plm_sph (AL, lmax, 0, T(unit_dir[2]));
@@ -131,15 +132,8 @@ namespace MR {
       }
 
 
-      template <typename T> inline Math::Vector<T>& delta (Math::Vector<T>& D, const Point& unit_dir, int lmax)
-      {
-        D.allocate (NforL (lmax));
-        delta (D.view(), unit_dir, lmax);
-        return (D);
-      }
 
-
-      template <typename T> inline Math::Vector<T>& SH2RH (Math::Vector<T>& RH, const Math::VectorView<T>& SH)
+      template <typename T> inline Math::Vector<T>& SH2RH (Math::Vector<T>& RH, const Math::Vector<T>& SH)
       {
         RH.allocate (SH.size());
         int lmax = 2*SH.size()+1;
@@ -151,10 +145,10 @@ namespace MR {
 
 
 
-      template <typename T> inline Math::VectorView<T>& sconv (Math::VectorView<T>& C, const Math::VectorView<T>& RH, const Math::VectorView<T>& SH)
+      template <typename T> inline Math::Vector<T>& sconv (Math::Vector<T>& C, const Math::Vector<T>& RH, const Math::Vector<T>& SH)
       {
         assert (SH.size() >= NforL (2*(RH.size()-1)));
-        assert (C.size() >= NforL (2*(RH.size()-1)));
+        C.allocate (NforL (2*(RH.size()-1)));
         for (int i = 0; i < int (RH.size()); ++i) {
           int l = 2*i;
           for (int m = -l; m <= l; ++m) 
@@ -163,12 +157,6 @@ namespace MR {
         return (C);
       }
 
-      template <typename T> inline Math::Vector<T>& sconv (Math::Vector<T>& C, const Math::VectorView<T>& RH, const Math::VectorView<T>& SH)
-      {
-        C.allocate (NforL (2*(RH.size()-1)));
-        sconv (C.view(), RH, SH);
-        return (C);
-      }
 
 
 
diff --git a/lib/math/cholesky.h b/lib/math/cholesky.h
index 01e5fc6380..73d0d94d15 100644
--- a/lib/math/cholesky.h
+++ b/lib/math/cholesky.h
@@ -49,32 +49,32 @@ namespace MR {
 
       //! %Cholesky decomposition of \a A into lower triangular matrix
       /** \note the contents of \a A will be overwritten with its %Cholesky decomposition */
-      template <typename T> inline MatrixView<T>& decomp (MatrixView<T>& A) {
-       	gsl_linalg_cholesky_decomp (&A); 
+      template <typename T> inline Matrix<T>& decomp (Matrix<T>& A) {
+       	gsl_linalg_cholesky_decomp (A.gsl()); 
 	return (A);
       }
 
       //! solve A*x = b given its %Cholesky decomposition \a D 
-      template <typename T> inline VectorView<T>& solve (VectorView<T>& x, const MatrixView<T>& D, const VectorView<T>& b) {
-        gsl_linalg_cholesky_solve (&D, &b, &x);
+      template <typename T> inline Vector<T>& solve (Vector<T>& x, const Matrix<T>& D, const Vector<T>& b) {
+        gsl_linalg_cholesky_solve (D.gsl(), b.gsl(), x.gsl());
         return (x);
       }
 
       //! solve A*x = b given its %Cholesky decomposition \a D, in place
       /** \note b is passed in as \a x */
-      template <typename T> inline VectorView<T>& solve (VectorView<T>& x, const MatrixView<T>& D) {
-        gsl_linalg_cholesky_svx (&D, &x);
+      template <typename T> inline Vector<T>& solve (Vector<T>& x, const Matrix<T>& D) {
+        gsl_linalg_cholesky_svx (D.gsl(), x.gsl());
         return (x);
       }
 
       //! invert A given its %Cholesky decomposition \a D, in place.
-      template <typename T> inline MatrixView<T>& inv_from_decomp (MatrixView<T>& D) {
-	gsl_linalg_cholesky_invert (&D);
+      template <typename T> inline Matrix<T>& inv_from_decomp (Matrix<T>& D) {
+	gsl_linalg_cholesky_invert (D.gsl());
 	return (D);
       }
 
       //! invert \a A using %Cholesky decomposition, in place.
-      template <typename T> inline MatrixView<T>& inv (MatrixView<T>& A) { return (inv_from_decomp (decomp (A))); }
+      template <typename T> inline Matrix<T>& inv (Matrix<T>& A) { return (inv_from_decomp (decomp (A))); }
 
       /** @} */
       /** @} */
diff --git a/lib/math/eigen.h b/lib/math/eigen.h
index bdc004d040..a2c743dfb3 100644
--- a/lib/math/eigen.h
+++ b/lib/math/eigen.h
@@ -46,11 +46,11 @@ namespace MR {
           Symm (size_t n) : work (gsl_eigen_symm_alloc(n)) { }
           ~Symm () { gsl_eigen_symm_free (work); }
 
-          VectorView<T>& operator () (VectorView<T>& eval, MatrixView<T>& A) {
+          Vector<T>& operator () (Vector<T>& eval, Matrix<T>& A) {
             assert (A.rows() == A.columns());
             assert (A.rows() == eval.size());
             assert (A.rows() == work->size);
-            int status = gsl_eigen_symm (&A, &eval, work);
+            int status = gsl_eigen_symm (A.gsl(), eval.gsl(), work);
             if (status) throw Exception (std::string ("eigenvalue decomposition failed: ") + gsl_strerror (status)); 
             return (eval);
           }
@@ -68,13 +68,13 @@ namespace MR {
           SymmV (size_t n) : work (gsl_eigen_symmv_alloc(n)) { }
           ~SymmV () { gsl_eigen_symmv_free (work); }
 
-          VectorView<T>& operator () (VectorView<T>& eval, MatrixView<T>& A, MatrixView<T>& evec) {
+          Vector<T>& operator () (Vector<T>& eval, Matrix<T>& A, Matrix<T>& evec) {
             assert (A.rows() == A.columns());
             assert (A.rows() == evec.rows());
             assert (evec.rows() == evec.columns());
             assert (A.rows() == eval.size());
             assert (A.rows() == work->size);
-            int status = gsl_eigen_symmv (&A, &eval, &evec, work);
+            int status = gsl_eigen_symmv (A.gsl(), eval.gsl(), evec.gsl(), work);
             if (status) throw Exception (std::string ("eigenvalue decomposition failed: ") + gsl_strerror (status)); 
             return (eval);
           }
@@ -84,10 +84,10 @@ namespace MR {
 
 
       //! Eigenvalue sorting
-      inline VectorView<double>& sort (VectorView<double>& eval) { gsl_sort_vector (&eval); return (eval); } 
+      inline Vector<double>& sort (Vector<double>& eval) { gsl_sort_vector (eval.gsl()); return (eval); } 
       //! Eigenvalue sorting
-      template <typename T> inline VectorView<T>& sort (VectorView<T>& eval, MatrixView<T>& evec) { 
-        gsl_eigen_symmv_sort (&eval, &evec, GSL_EIGEN_SORT_VAL_ASC);
+      template <typename T> inline Vector<T>& sort (Vector<T>& eval, Matrix<T>& evec) { 
+        gsl_eigen_symmv_sort (eval.gsl(), evec.gsl(), GSL_EIGEN_SORT_VAL_ASC);
         return (eval);
       }
 
diff --git a/lib/math/least_squares.h b/lib/math/least_squares.h
index 27bea3d542..99d26e8512 100644
--- a/lib/math/least_squares.h
+++ b/lib/math/least_squares.h
@@ -37,7 +37,7 @@ namespace MR {
 
 
     //! solve least-squares problem Mx = b
-    template <typename T> inline Vector<T>& solve_LS (Vector<T>& x, const MatrixView<T>& M, const VectorView<T>& b, MatrixView<T>& work) {
+    template <typename T> inline Vector<T>& solve_LS (Vector<T>& x, const Matrix<T>& M, const Vector<T>& b, Matrix<T>& work) {
       rankN_update (work, M, CblasTrans);
       Cholesky::decomp (work);
       mult (x, T(1.0), CblasTrans, M, b);
@@ -45,26 +45,20 @@ namespace MR {
     }
 
     //! compute Moore-Penrose pseudo-inverse of M given its transpose Mt
-    template <typename T> inline MatrixView<T>& pinv (MatrixView<T>& I, const MatrixView<T>& Mt, MatrixView<T>& work) {
+    template <typename T> inline Matrix<T>& pinv (Matrix<T>& I, const Matrix<T>& Mt, Matrix<T>& work) {
       mult (work, T(0.0), T(1.0), CblasNoTrans, Mt, CblasTrans, Mt);
       Cholesky::inv (work);
       return (mult (I, CblasLeft, T(0.0), T(1.0), CblasUpper, work, Mt));
     }
 
     //! compute Moore-Penrose pseudo-inverse of M 
-    template <typename T> inline MatrixView<T>& pinv (MatrixView<T>& I, const MatrixView<T>& M) {
+    template <typename T> inline Matrix<T>& pinv (Matrix<T>& I, const Matrix<T>& M) {
+      I.allocate (M.columns(), M.rows());
       Matrix<T> work (M.columns(), M.columns());
       Matrix<T> Mt (M.columns(), M.rows());
       return (pinv (I, transpose(Mt, M), work));
     }
 
-    //! compute Moore-Penrose pseudo-inverse of M 
-    template <typename T> inline Matrix<T>& pinv (Matrix<T>& I, const MatrixView<T>& M) {
-      I.allocate (M.columns(), M.rows());
-      pinv (I.view(), M);
-      return (I);
-    }
-
     /** @} */
     /** @} */
 
diff --git a/lib/math/math.h b/lib/math/math.h
index cc2cccb974..476385872a 100644
--- a/lib/math/math.h
+++ b/lib/math/math.h
@@ -62,6 +62,12 @@ namespace MR {
     DEFINE_ELEMENTARY_FUNCTION (acos);
     DEFINE_ELEMENTARY_FUNCTION (asin);
     DEFINE_ELEMENTARY_FUNCTION (atan);
+    DEFINE_ELEMENTARY_FUNCTION (cosh);
+    DEFINE_ELEMENTARY_FUNCTION (sinh);
+    DEFINE_ELEMENTARY_FUNCTION (tanh);
+    DEFINE_ELEMENTARY_FUNCTION (acosh);
+    DEFINE_ELEMENTARY_FUNCTION (asinh);
+    DEFINE_ELEMENTARY_FUNCTION (atanh);
     DEFINE_ELEMENTARY_FUNCTION (round);
     DEFINE_ELEMENTARY_FUNCTION (floor);
     DEFINE_ELEMENTARY_FUNCTION (ceil);
diff --git a/lib/math/matrix.h b/lib/math/matrix.h
index 5abb0f77e8..df3fa6af12 100644
--- a/lib/math/matrix.h
+++ b/lib/math/matrix.h
@@ -54,26 +54,133 @@ namespace MR {
     /** @addtogroup linalg 
       @{ */
 
-
-    //! provides access to data as a matrix
-    /** \note this class is not capable of managing its own data allocation. 
-      See the Matrix class for a more general interface. */
-    template <typename T> class MatrixView : protected GSLMatrix<T> 
+    //! A matrix class
+    template <typename T> class Matrix : protected GSLMatrix<T> 
     {
       public:
-        template <typename U> friend class MatrixView;
+        template <typename U> friend class Matrix;
+
+	//! construct empty matrix
+        Matrix () throw () {
+	  GSLMatrix<T>::size1 = GSLMatrix<T>::size2 = GSLMatrix<T>::tda = 0;
+	  GSLMatrix<T>::data = NULL; GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 1; 
+        }
+
+	//! copy constructor
+        /*! \note the new instance will reference the same data as the
+         * original, but will not be responsible for (or even capable of) any
+         * operation requiring data allocation or deallocation. */
+        Matrix (const Matrix<T>& M) throw () { 
+	  GSLMatrix<T>::size1 = M.GSLMatrix<T>::size1;
+          GSLMatrix<T>::size2 = M.GSLMatrix<T>::size2; 
+          GSLMatrix<T>::tda = M.GSLMatrix<T>::tda;
+	  GSLMatrix<T>::data = M.GSLMatrix<T>::data; 
+          GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 0; 
+        }
+
+	//! construct matrix of size \a nrows by \a ncolumns
+	/** \note the elements of the matrix are left uninitialised. */
+        Matrix (size_t nrows, size_t ncolumns) { 
+          if (nrows && ncolumns) {
+            GSLMatrix<T>::block = GSLBlock<T>::alloc (nrows * ncolumns);
+            if (!GSLMatrix<T>::block) throw Exception ("Failed to allocate memory for Matrix data");  
+          }
+          else GSLMatrix<T>::block = NULL;
+	  GSLMatrix<T>::size1 = nrows; GSLMatrix<T>::size2 = GSLMatrix<T>::tda = ncolumns; GSLMatrix<T>::owner = 1;
+          GSLMatrix<T>::data = GSLMatrix<T>::block ? GSLMatrix<T>::block->data : NULL;
+        }
 
 	//! construct from existing data array
-        explicit MatrixView (T* data, size_t nrows, size_t ncolumns) throw () { 
+        /*! \note with this constructor, the matrix is not responsible for
+         * managing data allocation or deallocation. The matrix will simply
+         * allow access to the data provided using the Matrix interface. The
+         * underlying data array must remain accessible during the lifetime of
+         * the Matrix class. */
+        Matrix (T* data, size_t nrows, size_t ncolumns) throw () { 
 	  GSLMatrix<T>::size1 = nrows; GSLMatrix<T>::size2 = GSLMatrix<T>::tda = ncolumns;
 	  GSLMatrix<T>::set(data); GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 0; 
         }
 	//! construct from existing data array with non-standard row stride
-        explicit MatrixView (T* data, size_t nrows, size_t ncolumns, size_t row_skip) throw () {  
+        /*! \note with this constructor, the matrix is not responsible for
+         * managing data allocation or deallocation. The matrix will simply
+         * allow access to the data provided using the Matrix interface. The
+         * underlying data array must remain accessible during the lifetime of
+         * the Matrix class. */
+        Matrix (T* data, size_t nrows, size_t ncolumns, size_t row_skip) throw () {  
 	  GSLMatrix<T>::size1 = nrows; GSLMatrix<T>::size2 = ncolumns; GSLMatrix<T>::tda = row_skip;
 	  GSLMatrix<T>::set (data); GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 0; 
         }
 
+	//! construct a matrix by reading from the text file \a filename
+        Matrix (const std::string& filename) { load (filename); } 
+
+	//! destructor
+        ~Matrix () { if (GSLMatrix<T>::block) { assert (GSLMatrix<T>::owner); GSLBlock<T>::free (GSLMatrix<T>::block); } }
+
+	//! deallocate the matrix data
+        Matrix<T>& clear () { 
+          if (GSLMatrix<T>::block) {
+            assert (GSLMatrix<T>::owner);
+	    GSLBlock<T>::free (GSLMatrix<T>::block);
+          }
+	  GSLMatrix<T>::size1 = GSLMatrix<T>::size2 = GSLMatrix<T>::tda = 0;
+	  GSLMatrix<T>::data = NULL; GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 1; 
+          return (*this);
+        }
+	//! allocate the matrix to have the same size as \a M
+        Matrix<T>& allocate (const Matrix<T>& M) { allocate (M.rows(), M.columns()); return (*this); }
+	//! allocate the matrix to have the same size as \a M
+        template <typename U> Matrix& allocate (const Matrix<U>& M) { allocate (M.rows(), M.columns()); return (*this); }
+
+	//! allocate the matrix to have size \a nrows by \a ncolumns
+        Matrix<T>& allocate (size_t nrows, size_t ncolumns) {
+          if (rows() == nrows && columns() == ncolumns) return (*this);
+          if (GSLMatrix<T>::block) {
+            assert (GSLMatrix<T>::owner);
+            if (GSLMatrix<T>::block->size < nrows * ncolumns) { 
+	      GSLBlock<T>::free (GSLMatrix<T>::block);
+	      GSLMatrix<T>::block = NULL;
+            }
+          }
+
+          if (!GSLMatrix<T>::block && nrows*ncolumns) {
+            GSLMatrix<T>::block = GSLBlock<T>::alloc (nrows * ncolumns);
+            if (!GSLMatrix<T>::block) throw Exception ("Failed to allocate memory for Matrix data");  
+          }
+	  GSLMatrix<T>::size1 = nrows; GSLMatrix<T>::size2 = GSLMatrix<T>::tda = ncolumns; GSLMatrix<T>::owner = 1;
+          GSLMatrix<T>::data = GSLMatrix<T>::block ? GSLMatrix<T>::block->data : NULL;
+          return (*this);
+        }
+
+        //! copy contents of \a M into current matrix
+        /*! \note this will also perform any reallocation required */
+        Matrix<T>& copy (const Matrix<T>& M) { allocate (M); return (operator= (M)); }
+
+        //! copy contents of \a M into current matrix
+        /*! \note this will also perform any reallocation required */
+        template <typename U> Matrix<T>& copy (const Matrix<U>& M) { allocate (M); return (operator= (M)); }
+
+	//! resize matrix to have size \a nrows by \a ncolumns
+        /*! \note this function will not perform any data reallocation. If the
+         * size requested exceeds the original size of the data, an exception
+         * will be thrown. */
+        Matrix<T>& resize (size_t nrows, size_t ncolumns) {
+          if (ncolumns > GSLMatrix<T>::tda || GSLMatrix<T>::block->size < nrows*GSLMatrix<T>::tda) 
+            throw Exception ("cannot resize matrix: size requested exceeds dimensions originally allocated");
+          GSLMatrix<T>::size1 = nrows;
+          GSLMatrix<T>::size2 = ncolumns;
+          return (*this); 
+        }
+
+	//! read matrix data from the text file \a filename
+        Matrix<T>& load (const std::string& filename) {
+          std::ifstream in (filename.c_str());
+          if (!in) throw Exception ("cannot open matrix file \"" + filename + "\": " + strerror (errno));
+          try { in >> *this; }
+          catch (Exception& E) { throw Exception (E, "error loading matrix file \"" + filename + "\""); }
+          return (*this);
+        }
+
 	//! write to text file \a filename
         void save (const std::string& filename) const {
           std::ofstream out (filename.c_str());
@@ -82,151 +189,175 @@ namespace MR {
         }
 
 	//! used to obtain a pointer to the underlying GSL structure
-	GSLMatrix<T>* operator& () { return (this); }
+	GSLMatrix<T>* gsl () { return (this); }
+
 	//! used to obtain a pointer to the underlying GSL structure
-        const GSLMatrix<T>* operator& () const { return (this); }
+        const GSLMatrix<T>* gsl () const { return (this); }
 
 	//! assign the specified \a value to all elements of the matrix
-        MatrixView<T>&  operator= (T value) throw () { LOOP (operator()(i,j) = value); return (*this); }
+        Matrix<T>&  operator= (T value) throw () { LOOP (operator()(i,j) = value); return (*this); }
+
 	//! assign the values in \a M to the corresponding elements of the matrix
-        MatrixView<T>&  operator= (const MatrixView& M) throw () { 
+        Matrix<T>&  operator= (const Matrix& M) throw () { 
           assert (rows() == M.rows() && columns() == M.columns());
           LOOP (operator()(i,j) = M(i,j)); 
           return (*this); 
         }
+
 	//! assign the values in \a M to the corresponding elements of the matrix
-        template <typename U> MatrixView<T>& operator= (const MatrixView<U>& M) {
+        template <typename U> Matrix<T>& operator= (const Matrix<U>& M) {
           assert (rows() == M.rows() && columns() == M.columns());
           LOOP (operator()(i,j) = M(i,j)); 
           return (*this); 
         }
 
 	//! comparison operator
-        bool          operator== (const MatrixView& M) throw () { return (! (*this != M)); }
+        bool operator== (const Matrix& M) throw () { return (! (*this != M)); }
+
 	//! comparison operator
-        bool          operator!= (const MatrixView& M) throw ()  { 
+        bool operator!= (const Matrix& M) throw ()  { 
           if (rows() != M.rows() || columns() != M.columns()) return (true);
           LOOP (if (operator()(i,j) != M(i,j)) return (true)); 
           return (false);
         }
+
 	//! swap contents with \a M without copying
-        void          swap (MatrixView& M) throw () { 
-          char c [sizeof (MatrixView)];
-          memcpy (&c, this, sizeof (MatrixView));
-          memcpy (this, &M, sizeof (MatrixView));
-          memcpy (&M, &c, sizeof (MatrixView));
+        void swap (Matrix& M) throw () { 
+          char c [sizeof (Matrix)];
+          memcpy (&c, this, sizeof (Matrix));
+          memcpy (this, &M, sizeof (Matrix));
+          memcpy (&M, &c, sizeof (Matrix));
         }
 
 	//! return reference to element at \a i, \a j
-        T&       operator() (size_t i, size_t j) throw () { 
+        T& operator() (size_t i, size_t j) throw () { 
 	  assert (i < rows()); assert (j < columns()); 
 	  return (ptr()[i * GSLMatrix<T>::tda + j]); 
 	}
-	//! return reference to element at \a i, \a j
+
+	//! return const reference to element at \a i, \a j
         const T& operator() (size_t i, size_t j) const throw () { 
 	  assert (i < rows()); assert (j < columns()); 
 	  return (ptr()[i * GSLMatrix<T>::tda + j]); 
 	}
+
 	//! return number of rows of matrix
-        size_t        rows () const throw () { return (GSLMatrix<T>::size1); }
+        size_t rows () const throw () { return (GSLMatrix<T>::size1); }
+
 	//! return number of columns of matrix
-        size_t        columns () const throw () { return (GSLMatrix<T>::size2); }
+        size_t columns () const throw () { return (GSLMatrix<T>::size2); }
 
 	//! true if matrix points to existing data
-        bool          is_set () const throw () { return (ptr()); }
+        bool is_set () const throw () { return (ptr()); }
+
+	//! return a pointer to the underlying data
+        T* ptr () throw () { return ((T*) (GSLMatrix<T>::data)); }
+
+	//! return a pointer to the underlying data
+        const T* ptr () const throw () { return ((const T*) (GSLMatrix<T>::data)); }
+
+	//! return the row stride
+        size_t row_stride () const throw () { return (GSLMatrix<T>::tda); }
+
 	//! set all elements of matrix to zero
-        MatrixView<T>&   zero () { LOOP (operator()(i,j) = 0.0); return (*this); }
+        Matrix<T>& zero () { LOOP (operator()(i,j) = 0.0); return (*this); }
+
 	//! set all diagonal elements of matrix to one, and all others to zero
-        MatrixView<T>&   identity () throw () { LOOP (operator()(i,j) = (i==j?1.0:0.0)); return (*this); }
+        Matrix<T>& identity () throw () { LOOP (operator()(i,j) = (i==j?1.0:0.0)); return (*this); }
 
 	//! add \a value to all elements of the matrix
-        MatrixView<T>& operator+= (T value) throw () { LOOP (operator()(i,j) += value); return (*this); }
+        Matrix<T>& operator+= (T value) throw () { LOOP (operator()(i,j) += value); return (*this); }
 	//! subtract \a value from all elements of the matrix
-        MatrixView<T>& operator-= (T value) throw () { LOOP (operator()(i,j) -= value); return (*this); }
+        Matrix<T>& operator-= (T value) throw () { LOOP (operator()(i,j) -= value); return (*this); }
 	//! multiply all elements of the matrix by \a value 
-        MatrixView<T>& operator*= (T value) throw () { LOOP (operator()(i,j) *= value); return (*this); }
+        Matrix<T>& operator*= (T value) throw () { LOOP (operator()(i,j) *= value); return (*this); }
 	//! divide all elements of the matrix by \a value 
-        MatrixView<T>& operator/= (T value) throw () { LOOP (operator()(i,j) /= value); return (*this); }
+        Matrix<T>& operator/= (T value) throw () { LOOP (operator()(i,j) /= value); return (*this); }
 
 	//! add each element of \a M to the corresponding element of the matrix
-        MatrixView<T>& operator+= (const MatrixView& M) throw () { 
+        Matrix<T>& operator+= (const Matrix& M) throw () { 
           assert (rows() == M.rows() && columns() == M.columns());
           LOOP (operator()(i,j) += M(i,j));
           return (*this); 
         }
+
 	//! subtract each element of \a M from the corresponding element of the matrix
-        MatrixView<T>& operator-= (const MatrixView& M) throw () { 
+        Matrix<T>& operator-= (const Matrix& M) throw () { 
           assert (rows() == M.rows() && columns() == M.columns());
           LOOP (operator()(i,j) -= M(i,j)); 
           return (*this);
         }
+
 	//! multiply each element of the matrix by the corresponding element of \a M
-        MatrixView<T>& operator*= (const MatrixView& M) throw () {
+        Matrix<T>& operator*= (const Matrix& M) throw () {
           assert (rows() == M.rows() && columns() == M.columns());
           LOOP (operator()(i,j) *= M(i,j)); 
           return (*this);
         }
+
 	//! divide each element of the matrix by the corresponding element of \a M
-        MatrixView<T>& operator/= (const MatrixView& M) throw () { 
+        Matrix<T>& operator/= (const Matrix& M) throw () { 
           assert (rows() == M.rows() && columns() == M.columns());
           LOOP (operator()(i,j) /= M(i,j));
           return (*this); 
         }
 
-	//! return a MatrixView corresponding to a submatrix of the matrix
-        MatrixView<T> view (size_t from_row, size_t to_row, size_t from_column, size_t to_column) throw () {
+	//! return a Matrix corresponding to a submatrix of the matrix
+        Matrix<T> sub (size_t from_row, size_t to_row, size_t from_column, size_t to_column) throw () {
           assert (from_row <= to_row && to_row <= rows());
           assert (from_column <= to_column && to_column <= columns());
-          MatrixView<T> M (MatrixView (ptr() + from_row*GSLMatrix<T>::tda + from_column, to_row-from_row, to_column-from_column, GSLMatrix<T>::tda));
-          return (M);
+          return (Matrix<T> (Matrix (ptr() + from_row*GSLMatrix<T>::tda + from_column, 
+                  to_row-from_row, to_column-from_column, GSLMatrix<T>::tda)));
         }
 
-	//! return a MatrixView corresponding to a submatrix of the matrix
-        const MatrixView<T> view (size_t from_row, size_t to_row, size_t from_column, size_t to_column) const throw () {
+	//! return a Matrix corresponding to a submatrix of the matrix
+        const Matrix<T> sub (size_t from_row, size_t to_row, size_t from_column, size_t to_column) const throw () {
           assert (from_row <= to_row && to_row <= rows());
           assert (from_column <= to_column && to_column <= columns());
-          MatrixView<T> M (MatrixView (ptr() + from_row*GSLMatrix<T>::tda + from_column, to_row-from_row, to_column-from_column, GSLMatrix<T>::tda));
-          return (M);
+          return (Matrix<T> (Matrix (ptr() + from_row*GSLMatrix<T>::tda + from_column, 
+                  to_row-from_row, to_column-from_column, GSLMatrix<T>::tda)));
         }
 
-	//! return a MatrixView corresponding to a row of the matrix
-        VectorView<T> row (size_t index = 0) throw () { 
+	//! return a Matrix corresponding to a row of the matrix
+        Vector<T> row (size_t index = 0) throw () { 
 	  assert (index < rows()); 
-	  return (VectorView<T> (ptr()+index*GSLMatrix<T>::tda, GSLMatrix<T>::size2, 1)); 
+	  return (Vector<T> (ptr()+index*GSLMatrix<T>::tda, GSLMatrix<T>::size2, 1)); 
 	}
-	//! return a MatrixView corresponding to a row of the matrix
-        const VectorView<T> row (size_t index = 0) const throw () { 
+	//! return a Matrix corresponding to a row of the matrix
+        const Vector<T> row (size_t index = 0) const throw () { 
 	  assert (index < rows()); 
-	  return (VectorView<T> (const_cast<T*> (ptr())+index*GSLMatrix<T>::tda, GSLMatrix<T>::size2, 1)); 
+	  return (Vector<T> (const_cast<T*> (ptr())+index*GSLMatrix<T>::tda, GSLMatrix<T>::size2, 1)); 
 	}
 
-	//! return a MatrixView corresponding to a column of the matrix
-        VectorView<T> column (size_t index = 0)  throw () { 
+	//! return a Matrix corresponding to a column of the matrix
+        Vector<T> column (size_t index = 0)  throw () { 
 	  assert (index < columns());
-	  return (VectorView<T> (ptr()+index, GSLMatrix<T>::size1, GSLMatrix<T>::tda)); 
+	  return (Vector<T> (ptr()+index, GSLMatrix<T>::size1, GSLMatrix<T>::tda)); 
 	}
-	//! return a MatrixView corresponding to a column of the matrix
-        const VectorView<T> column (size_t index = 0) const throw () { 
+
+	//! return a Matrix corresponding to a column of the matrix
+        const Vector<T> column (size_t index = 0) const throw () { 
 	  assert (index < columns());
-	  return (VectorView<T> (const_cast<T*> (ptr())+index, GSLMatrix<T>::size1, GSLMatrix<T>::tda)); 
+	  return (Vector<T> (const_cast<T*> (ptr())+index, GSLMatrix<T>::size1, GSLMatrix<T>::tda)); 
 	}
 
-	//! return a MatrixView corresponding to the diagonal of the matrix
-        VectorView<T> diagonal () const throw () {
+	//! return a Matrix corresponding to the diagonal of the matrix
+        Vector<T> diagonal () const throw () {
 	  assert (rows() > 0 && columns() > 0);
-	  return (VectorView<T> (ptr(), MIN(GSLMatrix<T>::size1,GSLMatrix<T>::size2), GSLMatrix<T>::tda+1)); 
+	  return (Vector<T> (ptr(), MIN(GSLMatrix<T>::size1,GSLMatrix<T>::size2), GSLMatrix<T>::tda+1)); 
 	}
-	//! return a MatrixView corresponding to a diagonal of the matrix
+
+	//! return a Matrix corresponding to a diagonal of the matrix
 	/** \param offset the diagonal to obtain. If \a offset > 0, return the corresponding upper diagonal. 
 	  If \a offset < 0, return the corresponding lower diagonal. */ 
-        VectorView<T> diagonal (int offset) const throw () {
+        Vector<T> diagonal (int offset) const throw () {
           assert (rows() > 0 && columns() > 0);
           if (offset == 0) return (diagonal());
-          if (offset < 0) return (VectorView<T> (
+          if (offset < 0) return (Vector<T> (
 		ptr()-GSLMatrix<T>::tda*offset, 
 		MIN(GSLMatrix<T>::size1+offset,GSLMatrix<T>::size2+offset), 
 		GSLMatrix<T>::tda+1));
-          return (VectorView<T> (
+          return (Vector<T> (
 		ptr()+offset,
 	       	MIN(GSLMatrix<T>::size1-offset,GSLMatrix<T>::size2-offset),
 	       	GSLMatrix<T>::tda+1));
@@ -235,37 +366,18 @@ namespace MR {
 	//! swap two rows of the matrix
         void swap_rows (size_t n, size_t m) {
 	  for (size_t i = 0; i < GSLMatrix<T>::size2; i++)
-	    Math::swap (operator()(n,i), operator()(m,i)); 
+	    std::swap (operator()(n,i), operator()(m,i)); 
 	}
+
 	//! swap two columns of the matrix
         void swap_columns (size_t n, size_t m) { 
 	  for (size_t i = 0; i < GSLMatrix<T>::size1; i++) 
-	    Math::swap (operator()(i,n), operator()(i,m)); 
+	    std::swap (operator()(i,n), operator()(i,m)); 
 	}
 
-	//! resize the matrix to have size \a nrows by \a ncolumns
-	/** \note no bounds checking and no data allocation/deallocation/copying is performed by this function. */
-        MatrixView<T>& resize (size_t nrows, size_t ncolumns) { GSLMatrix<T>::size1 = nrows; GSLMatrix<T>::size2 = ncolumns; return (*this); }
-	//! resize the matrix to refer the submatrix specified
-	/** \note no bounds checking and no data allocation/deallocation/copying is performed by this function. */
-        MatrixView<T>& resize (size_t from_row, size_t to_row, size_t from_column, size_t to_column) throw () {
-	  ptr() += from_row*GSLMatrix<T>::tda + from_column;
-	  GSLMatrix<T>::size1 = to_row-from_row;
-	  GSLMatrix<T>::size2 = to_column-from_column;
-          return (*this);
-        }
-
-	//! return a pointer to the underlying data
-        T* ptr () throw () { return ((T*) (GSLMatrix<T>::data)); }
-	//! return a pointer to the underlying data
-        const T* ptr () const throw () { return ((const T*) (GSLMatrix<T>::data)); }
-	//! return the row stride
-        size_t row_stride () const throw () { return (GSLMatrix<T>::tda); }
-
 
 	//! write the matrix \a M to \a stream as text
-        friend std::ostream& operator<< (std::ostream& stream, const MatrixView& M)
-        {
+        friend std::ostream& operator<< (std::ostream& stream, const Matrix& M) {
           for (size_t i = 0; i < M.rows(); i++) {
             for (size_t j = 0; j < M.columns(); j++) 
               stream << M(i,j) << " ";
@@ -274,133 +386,8 @@ namespace MR {
           return (stream);
         }
 
-
-      protected:
-        MatrixView () throw () {
-	  GSLMatrix<T>::size1 = GSLMatrix<T>::size2 = GSLMatrix<T>::tda = 0;
-	  GSLMatrix<T>::data = NULL; GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 0; }
-    };
-
-
-
-
-
-    //! The main matrix class 
-    template <typename T> class Matrix : public MatrixView<T> 
-    {
-      public:
-        template <typename U> friend class Matrix;
-
-	//! construct empty matrix
-        explicit Matrix () throw () { }
-	//! construct matrix of size \a nrows by \a ncolumns
-	/** \note the elements of the matrix are left uninitialised. */
-        explicit Matrix (size_t nrows, size_t ncolumns) { allocate (nrows, ncolumns); }
-	//! construct a matrix by copying the elements of \a M
-        explicit Matrix (const MatrixView<T>& M) { allocate (M); view() = M; }
-	//! construct a matrix by copying the elements of \a M
-        template <typename U> Matrix (const MatrixView<U>& M) { allocate (M); view() = M; }
-        explicit Matrix (const Matrix<T>& M) : MatrixView<T>() { allocate (M); view() = M.view(); }
-        template <typename U> Matrix (const Matrix<U>& M) : MatrixView<T>() { allocate (M); view() = M.view(); }
-
-
-	//! construct a matrix by reading from the text file \a filename
-        explicit Matrix (const std::string& filename) { load (filename); } 
-	//! destructor
-        ~Matrix () { if (GSLMatrix<T>::block) { assert (GSLMatrix<T>::owner); GSLBlock<T>::free (GSLMatrix<T>::block); } }
-
-	//! assignment operator: allocate matrix to have the same size as \a M, and copy contents from \a M
-        Matrix<T>&       operator= (const Matrix<T>& M) { return (operator= (M.view())); }
-	//! assignment operator: allocate matrix to have the same size as \a M, and copy contents from \a M
-        Matrix<T>&       operator= (const MatrixView<T>& M) { allocate (M); view() = M; return (*this); }
-	//! assignment operator: allocate matrix to have the same size as \a M, and copy contents from \a M
-        template <typename U> Matrix& operator= (const Matrix<U>& M) { return (operator=(M.view())); }
-	//! assignment operator: allocate matrix to have the same size as \a M, and copy contents from \a M
-        template <typename U> Matrix& operator= (const MatrixView<U>& M) { allocate (M); view() = M; return (*this); }
-
-	//! deallocate the matrix data
-        Matrix<T>&       clear () { 
-          if (GSLMatrix<T>::block) {
-            assert (GSLMatrix<T>::owner);
-	    GSLBlock<T>::free (GSLMatrix<T>::block);
-          }
-	  GSLMatrix<T>::size1 = GSLMatrix<T>::size2 = GSLMatrix<T>::tda = 0;
-	  GSLMatrix<T>::data = NULL; GSLMatrix<T>::block = NULL; GSLMatrix<T>::owner = 0; 
-          return (*this);
-        }
-	//! allocate the matrix to have the same size as \a M
-        Matrix<T>&       allocate (const MatrixView<T>& M) { allocate (M.rows(), M.columns()); return (*this); }
-	//! allocate the matrix to have the same size as \a M
-        template <typename U> Matrix& allocate (const MatrixView<U>& M) { allocate (M.rows(), M.columns()); return (*this); }
-	//! allocate the matrix to have size \a nrows by \a ncolumns
-        Matrix<T>&       allocate (size_t nrows, size_t ncolumns) {
-          if (GSLMatrix<T>::block) {
-            assert (GSLMatrix<T>::owner);
-            if (GSLMatrix<T>::block->size < nrows * ncolumns) { 
-	      GSLBlock<T>::free (GSLMatrix<T>::block);
-	      GSLMatrix<T>::block = NULL;
-            }
-          }
-
-          if (!GSLMatrix<T>::block && nrows*ncolumns) {
-            GSLMatrix<T>::block = GSLBlock<T>::alloc (nrows * ncolumns);
-            if (!GSLMatrix<T>::block) throw Exception ("Failed to allocate memory for Matrix data");  
-          }
-	  GSLMatrix<T>::size1 = nrows; GSLMatrix<T>::size2 = GSLMatrix<T>::tda = ncolumns; GSLMatrix<T>::owner = 1;
-          GSLMatrix<T>::data = GSLMatrix<T>::block ? GSLMatrix<T>::block->data : NULL;
-          return (*this);
-        }
-
-	//! resize matrix to have size \a nrows by \a ncolumns
-        Matrix<T>&        resize (size_t nrows, size_t ncolumns) {
-          if (nrows <= GSLMatrix<T>::size1 || ncolumns <= GSLMatrix<T>::size2) { 
-            GSLMatrix<T>::size1 = nrows;
-            GSLMatrix<T>::size2 = ncolumns;
-            return (*this); 
-          }
-          if (GSLMatrix<T>::block) {
-            assert (GSLMatrix<T>::owner);
-            if (GSLMatrix<T>::block->size >= (nrows-1)*GSLMatrix<T>::tda + ncolumns) {
-              GSLMatrix<T>::size1 = nrows; 
-              GSLMatrix<T>::size2 = ncolumns; 
-              return (*this);
-            }
-          }
-          Matrix M (nrows, ncolumns);
-          size_t nr = MIN(GSLMatrix<T>::size1, M.GSLMatrix<T>::size1), nc = MIN(GSLMatrix<T>::size2, M.GSLMatrix<T>::size2);
-          M.MatrixView<T>::view (0, nr, 0, nc) = MatrixView<T>::view (0, nr, 0, nc);
-          clear();
-          swap (M);
-          return (*this);
-        }
-
-	//! return as MatrixView class
-        MatrixView<T>& view () throw () { return (*this); }
-	//! return as MatrixView class
-        const MatrixView<T>& view () const throw () { return (*this); }
-
-	//! return a MatrixView corresponding to a submatrix of the matrix
-        MatrixView<T> view (size_t from_row, size_t to_row, size_t from_column, size_t to_column) throw () {
-	  return (MatrixView<T>::view (from_row, to_row, from_column, to_column));
-       	}
-	//! return a MatrixView corresponding to a submatrix of the matrix
-        const MatrixView<T> view (size_t from_row, size_t to_row, size_t from_column, size_t to_column) const throw () {
-	  return (MatrixView<T>::view (from_row, to_row, from_column, to_column));
-       	}
-
-	//! read matrix data from the text file \a filename
-        Matrix<T>& load (const std::string& filename) {
-          std::ifstream in (filename.c_str());
-          if (!in) throw Exception ("cannot open matrix file \"" + filename + "\": " + strerror (errno));
-          try { in >> *this; }
-          catch (Exception& E) { throw Exception (E, "error loading matrix file \"" + filename + "\""); }
-          return (*this);
-        }
-
-
 	//! read the matrix data from \a stream and assign to the matrix \a M
-        friend std::istream& operator>> (std::istream& stream, Matrix<T>& M)  
-        {
+        friend std::istream& operator>> (std::istream& stream, Matrix<T>& M)  {
           std::vector< RefPtr< std::vector<T> > > V;
           do {
             std::string sbuf;
@@ -443,6 +430,9 @@ namespace MR {
 
 
 
+
+
+
     /** @defgroup gemv Matrix-vector multiplication
       @{ */
 
@@ -454,11 +444,11 @@ namespace MR {
      *    - CblasNoTrans: use A
      *    - CblasTrans: use transpose of A
      *    - CblasConjTrans: use conjugate transpose of A
-     * \param A a MatrixView
-     * \param x a VectorView
+     * \param A a Matrix
+     * \param x a Vector
      * \return a reference to the target vector \a y
      */
-    template <typename T> inline VectorView<T>& mult (VectorView<T>& y, T beta, T alpha, CBLAS_TRANSPOSE op_A, const MatrixView<T>& A, const VectorView<T>& x) 
+    template <typename T> inline Vector<T>& mult (Vector<T>& y, T beta, T alpha, CBLAS_TRANSPOSE op_A, const Matrix<T>& A, const Vector<T>& x) 
     { 
       gemv (op_A, alpha, A, x, beta, y);
       return (y); 
@@ -471,39 +461,25 @@ namespace MR {
      *    - CblasNoTrans: use A
      *    - CblasTrans: use transpose of A
      *    - CblasConjTrans: use conjugate transpose of A
-     * \param A a MatrixView
+     * \param A a Matrix
      * \param x a Vector
      * \note this version will perform the appropriate allocation for \a y. 
      * \return a reference to the target vector \a y
      */
-    template <typename T> inline Vector<T>& mult (Vector<T>& y, T alpha, CBLAS_TRANSPOSE op_A, const MatrixView<T>& A, const VectorView<T>& x) { 
+    template <typename T> inline Vector<T>& mult (Vector<T>& y, T alpha, CBLAS_TRANSPOSE op_A, const Matrix<T>& A, const Vector<T>& x) { 
       y.allocate (A.rows());
-      mult (y.view(), T(0.0), alpha, op_A, A, x);
+      mult (y, T(0.0), alpha, op_A, A, x);
       return (y); 
     }
 
     //! Computes the matrix-vector product \a y = \a A \a x
     /** \param y the target vector
-     * \param A a MatrixView
+     * \param A a Matrix
      * \param x a Vector
      * \return a reference to the target vector \a y
      */
-    template <typename T> inline VectorView<T>& mult (VectorView<T>& y, const MatrixView<T>& A, const VectorView<T>& x) {
-      return (mult (y, T(0.0), T(1.0), CblasNoTrans, A, x)); 
-    }
-
-    //! Simplified matrix-vector multiplication with allocation
-    //! Computes the matrix-vector product \a y = \a A \a x, allocating storage for \a y
-    /** \param y the target vector
-     * \param A a MatrixView
-     * \param x a Vector
-     * \note this version will perform the appropriate allocation for \a y. 
-     * \return a reference to the target vector \a y
-     */
-    template <typename T> inline Vector<T>& mult (Vector<T>& y, const MatrixView<T>& A, const VectorView<T>& x) {
-      y.allocate (A.rows());
-      mult (y.view(), A, x); 
-      return (y);
+    template <typename T> inline Vector<T>& mult (Vector<T>& y, const Matrix<T>& A, const Vector<T>& x) {
+      return (mult (y, T(1.0), CblasNoTrans, A, x)); 
     }
 
     /** @} */
@@ -522,15 +498,15 @@ namespace MR {
      *    - CblasNoTrans: use A
      *    - CblasTrans: use transpose of A
      *    - CblasConjTrans: use conjugate transpose of A
-     * \param A a MatrixView
+     * \param A a Matrix
      * \param op_B determines how to use the matrix \a B:
      *    - CblasNoTrans: use B
      *    - CblasTrans: use transpose of B
      *    - CblasConjTrans: use conjugate transpose of B
-     * \param B a MatrixView
+     * \param B a Matrix
      * \return a reference to the target matrix \a C
      */
-    template <typename T> inline MatrixView<T>& mult (MatrixView<T>& C, T beta, T alpha, CBLAS_TRANSPOSE op_A, const MatrixView<T>& A, CBLAS_TRANSPOSE op_B, const MatrixView<T>& B) 
+    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, T beta, T alpha, CBLAS_TRANSPOSE op_A, const Matrix<T>& A, CBLAS_TRANSPOSE op_B, const Matrix<T>& B) 
     { 
       gemm (op_A, op_B, alpha, A, B, beta, C);
       return (C); 
@@ -543,42 +519,29 @@ namespace MR {
      *    - CblasNoTrans: use A
      *    - CblasTrans: use transpose of A
      *    - CblasConjTrans: use conjugate transpose of A
-     * \param A a MatrixView
+     * \param A a Matrix
      * \param op_B determines how to use the matrix \a B:
      *    - CblasNoTrans: use B
      *    - CblasTrans: use transpose of B
      *    - CblasConjTrans: use conjugate transpose of B
-     * \param B a MatrixView
+     * \param B a Matrix
      * \note this version will perform the appropriate allocation for \a C 
      * \return a reference to the target matrix \a C
      */
-    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, T alpha, CBLAS_TRANSPOSE op_A, const MatrixView<T>& A, CBLAS_TRANSPOSE op_B, const MatrixView<T>& B) { 
+    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, T alpha, CBLAS_TRANSPOSE op_A, const Matrix<T>& A, CBLAS_TRANSPOSE op_B, const Matrix<T>& B) { 
       C.allocate (A.rows(), B.columns());
-      mult (C.view(), T(0.0), alpha, op_A, A, op_B, B);
+      mult (C, T(0.0), alpha, op_A, A, op_B, B);
       return (C); 
     }
 
     //! computes the simplified general matrix-matrix multiplication \a C = \a A \a B
     /** \param C the target matrix
-     * \param A a MatrixView
-     * \param B a MatrixView
+     * \param A a Matrix
+     * \param B a Matrix
      * \return a reference to the target matrix \a C
      */
-    template <typename T> inline MatrixView<T>& mult (MatrixView<T>& C, const MatrixView<T>& A, const MatrixView<T>& B) { 
-      return (mult (C, T(0.0), T(1.0), CblasNoTrans, A, CblasNoTrans, B)); 
-    }
-
-    //! computes the simplified general matrix-matrix multiplication \a C = \a A \a B, allocating storage for \a C
-    /** \param C the target matrix
-     * \param A a MatrixView
-     * \param B a MatrixView
-     * \note this version will perform the appropriate allocation for \a C 
-     * \return a reference to the target matrix \a C
-     */
-    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, const MatrixView<T>& A, const MatrixView<T>& B) { 
-      C.allocate (A.rows(), B.columns());
-      mult (C.view(), A, B); 
-      return (C);
+    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, const Matrix<T>& A, const Matrix<T>& B) { 
+      return (mult (C, T(1.0), CblasNoTrans, A, CblasNoTrans, B)); 
     }
 
     /** @} */
@@ -603,11 +566,11 @@ namespace MR {
      * \param uplo determines which part of \a A will be used:
      *    - CblasUpper: upper triangle and diagonal of A
      *    - CblasLower: lower triangle and diagonal of A
-     * \param A a MatrixView
-     * \param B a MatrixView
+     * \param A a Matrix
+     * \param B a Matrix
      * \return a reference to the target matrix \a C
      */
-    template <typename T> inline MatrixView<T>& mult (MatrixView<T>& C, CBLAS_SIDE side, T beta, T alpha, CBLAS_UPLO uplo, const MatrixView<T>& A, const MatrixView<T>& B) 
+    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, CBLAS_SIDE side, T beta, T alpha, CBLAS_UPLO uplo, const Matrix<T>& A, const Matrix<T>& B) 
     { 
       symm (side, uplo, alpha, A, B, beta, C);
       return (C);
@@ -624,12 +587,12 @@ namespace MR {
      * \param uplo determines which part of \a A will be used:
      *    - CblasUpper: upper triangle and diagonal of A
      *    - CblasLower: lower triangle and diagonal of A
-     * \param A a MatrixView
-     * \param B a MatrixView
+     * \param A a Matrix
+     * \param B a Matrix
      * \note this version will perform the appropriate allocation for \a C 
      * \return a reference to the target matrix \a C
      */
-    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, CBLAS_SIDE side, T alpha, CBLAS_UPLO uplo, const MatrixView<T>& A, const MatrixView<T>& B) 
+    template <typename T> inline Matrix<T>& mult (Matrix<T>& C, CBLAS_SIDE side, T alpha, CBLAS_UPLO uplo, const Matrix<T>& A, const Matrix<T>& B) 
     { 
       C.allocate (A);
       symm (side, uplo, alpha, A, B, T(0.0), C);
@@ -642,7 +605,7 @@ namespace MR {
 
     //! solve for \a y in the triangular matrix problem: \a op_A(\a A) \a y = \a x
     /** \param x the target vector. On input, \a x should be set to \a y
-     * \param A a MatrixView
+     * \param A a Matrix
      * \param uplo determines which part of \a A will be used:
      *    - CblasUpper: upper triangle and diagonal of A
      *    - CblasLower: lower triangle and diagonal of A
@@ -655,7 +618,7 @@ namespace MR {
      *    - CblasNonUnit: diagonal elements of \a A are used
      * \return a reference to the target vector \a x
      */
-    template <typename T> inline VectorView<T>& solve_triangular (VectorView<T>& x, const MatrixView<T>& A, CBLAS_UPLO uplo = CblasUpper, CBLAS_TRANSPOSE op_A = CblasNoTrans, CBLAS_DIAG diag = CblasNonUnit)
+    template <typename T> inline Vector<T>& solve_triangular (Vector<T>& x, const Matrix<T>& A, CBLAS_UPLO uplo = CblasUpper, CBLAS_TRANSPOSE op_A = CblasNoTrans, CBLAS_DIAG diag = CblasNonUnit)
     {
       trsv (uplo, op_A, diag, A, x);
       return (x);
@@ -663,12 +626,12 @@ namespace MR {
 
     //! rank-1 update: \a A = \a alpha \a x \a y^T + \a A
     /** \param A the target matrix. 
-     * \param x a VectorView
-     * \param y a VectorView
+     * \param x a Vector
+     * \param y a Vector
      * \param alpha used to scale the product
      * \return a reference to the target matrix \a A
      */
-    template <typename T> inline MatrixView<T>& rank1_update (MatrixView<T>& A, const VectorView<T>& x, const VectorView<T>& y, T alpha = 1.0)
+    template <typename T> inline Matrix<T>& rank1_update (Matrix<T>& A, const Vector<T>& x, const Vector<T>& y, T alpha = 1.0)
     {
       ger (alpha, x, y, A);
       return (A);
@@ -676,14 +639,14 @@ namespace MR {
 
     //! symmetric rank-1 update: \a A = \a alpha \a x \a x^T + \a A, for symmetric \a A
     /** \param A the symmetric target matrix. 
-     * \param x a VectorView
+     * \param x a Vector
      * \param alpha used to scale the product
      * \param uplo determines which part of \a A will be used:
      *    - CblasUpper: upper triangle and diagonal of A
      *    - CblasLower: lower triangle and diagonal of A
      * \return a reference to the target matrix \a A
      */
-    template <typename T> inline MatrixView<T>& sym_rank1_update (MatrixView<T>& A, const VectorView<T>& x, T alpha = 1.0, CBLAS_UPLO uplo = CblasUpper)
+    template <typename T> inline Matrix<T>& sym_rank1_update (Matrix<T>& A, const Vector<T>& x, T alpha = 1.0, CBLAS_UPLO uplo = CblasUpper)
     {
       syr (uplo, alpha, x, A);
       return (A);
@@ -691,7 +654,7 @@ namespace MR {
 
     //! symmetric rank-N update: \a C = \a alpha \a op_A(\a A) op_A(\a A)^T + \a beta C, for symmetric \a A
     /** \param C the target matrix. If \a beta is non-zero, \a C should be symmetric
-     * \param A a MatrixView
+     * \param A a Matrix
      * \param op_A determines how to use the matrix \a A:
      *    - CblasNoTrans: use A
      *    - CblasTrans: use transpose of A
@@ -703,7 +666,7 @@ namespace MR {
      * \param beta used to add a multiple of \a C to the final product
      * \return a reference to the target matrix \a C
      */
-    template <typename T> inline MatrixView<T>& rankN_update (MatrixView<T>& C, const MatrixView<T>& A, CBLAS_TRANSPOSE op_A = CblasNoTrans, CBLAS_UPLO uplo = CblasUpper, T alpha = 1.0, T beta = 0.0)
+    template <typename T> inline Matrix<T>& rankN_update (Matrix<T>& C, const Matrix<T>& A, CBLAS_TRANSPOSE op_A = CblasNoTrans, CBLAS_UPLO uplo = CblasUpper, T alpha = 1.0, T beta = 0.0)
     {
       syrk (uplo, op_A, alpha, A, beta, C);
       return (C);
@@ -711,29 +674,18 @@ namespace MR {
 
     //! compute transpose \a A = \a B^T
     /** \param A the target matrix. 
-     * \param B a MatrixView to be transposed
+     * \param B a Matrix to be transposed
      * \return a reference to the target matrix \a A
+     * \note this version will perform the appropriate allocation for \a A 
      */
-    template <typename T> inline MatrixView<T>& transpose (MatrixView<T>& A, const MatrixView<T>& B)
+    template <typename T> inline Matrix<T>& transpose (Matrix<T>& A, const Matrix<T>& B)
     {
+      A.allocate (B.columns(), B.rows());
       for (size_t i = 0; i < B.rows(); i++)
         for (size_t j = 0; j < B.columns(); j++)
           A(j,i) = B(i,j);
       return (A);
     }
-
-    //! compute transpose \a A = \a B^T
-    /** \param A the target matrix. 
-     * \param B a MatrixView to be transposed
-     * \return a reference to the target matrix \a A
-     * \note this version will perform the appropriate allocation for \a A 
-     */
-    template <typename T> inline Matrix<T>& transpose (Matrix<T>& A, const MatrixView<T>& B)
-    {
-      A.allocate (B.columns(), B.rows());
-      transpose (A.view(), B);
-      return (A);
-    }
     /** @} */
 
 
@@ -746,50 +698,50 @@ namespace MR {
     namespace {
       // double definitions:
 
-      inline void gemm (CBLAS_TRANSPOSE op_A, CBLAS_TRANSPOSE op_B, double alpha, const MatrixView<double>& A, const MatrixView<double>& B, double beta, MatrixView<double>& C)
-      { gsl_blas_dgemm (op_A, op_B, alpha, &A, &B, beta, &C); }
+      inline void gemm (CBLAS_TRANSPOSE op_A, CBLAS_TRANSPOSE op_B, double alpha, const Matrix<double>& A, const Matrix<double>& B, double beta, Matrix<double>& C)
+      { gsl_blas_dgemm (op_A, op_B, alpha, A.gsl(), B.gsl(), beta, C.gsl()); }
 
-      inline void gemv (CBLAS_TRANSPOSE op_A, double alpha, const MatrixView<double>& A, const VectorView<double>& x, double beta, VectorView<double>& y)
-      { gsl_blas_dgemv (op_A, alpha, &A, &x, beta, &y); }
+      inline void gemv (CBLAS_TRANSPOSE op_A, double alpha, const Matrix<double>& A, const Vector<double>& x, double beta, Vector<double>& y)
+      { gsl_blas_dgemv (op_A, alpha, A.gsl(), x.gsl(), beta, y.gsl()); }
 
-      inline void symm (CBLAS_SIDE side, CBLAS_UPLO uplo, double alpha, const MatrixView<double>& A, const MatrixView<double>& B, double beta, MatrixView<double>& C)
-       { gsl_blas_dsymm (side, uplo, alpha, &A, &B, beta, &C); }
+      inline void symm (CBLAS_SIDE side, CBLAS_UPLO uplo, double alpha, const Matrix<double>& A, const Matrix<double>& B, double beta, Matrix<double>& C)
+       { gsl_blas_dsymm (side, uplo, alpha, A.gsl(), B.gsl(), beta, C.gsl()); }
 
-      inline void trsv (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, CBLAS_DIAG diag, const MatrixView<double>& A, VectorView<double>& x)
-      { gsl_blas_dtrsv (uplo, op_A, diag, &A, &x); }
+      inline void trsv (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, CBLAS_DIAG diag, const Matrix<double>& A, Vector<double>& x)
+      { gsl_blas_dtrsv (uplo, op_A, diag, A.gsl(), x.gsl()); }
 
-      inline void ger (double alpha, const VectorView<double>& x, const VectorView<double>& y, MatrixView<double>& A)
-      { gsl_blas_dger (alpha, &x, &y, &A); }
+      inline void ger (double alpha, const Vector<double>& x, const Vector<double>& y, Matrix<double>& A)
+      { gsl_blas_dger (alpha, x.gsl(), y.gsl(), A.gsl()); }
 
-      inline void syr (CBLAS_UPLO uplo, double alpha, const VectorView<double>& x, MatrixView<double>& A)
-      { gsl_blas_dsyr (uplo, alpha, &x, &A); }
+      inline void syr (CBLAS_UPLO uplo, double alpha, const Vector<double>& x, Matrix<double>& A)
+      { gsl_blas_dsyr (uplo, alpha, x.gsl(), A.gsl()); }
 
-      inline void syrk (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, double alpha, const MatrixView<double>& A, double beta, MatrixView<double>& C)
-      { gsl_blas_dsyrk (uplo, op_A, alpha, &A, beta, &C); }
+      inline void syrk (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, double alpha, const Matrix<double>& A, double beta, Matrix<double>& C)
+      { gsl_blas_dsyrk (uplo, op_A, alpha, A.gsl(), beta, C.gsl()); }
 
       // float definitions:
 
 
-      inline void gemm (CBLAS_TRANSPOSE op_A, CBLAS_TRANSPOSE op_B, float alpha, const MatrixView<float>& A, const MatrixView<float>& B, float beta, MatrixView<float>& C)
-      { gsl_blas_sgemm (op_A, op_B, alpha, &A, &B, beta, &C); }
+      inline void gemm (CBLAS_TRANSPOSE op_A, CBLAS_TRANSPOSE op_B, float alpha, const Matrix<float>& A, const Matrix<float>& B, float beta, Matrix<float>& C)
+      { gsl_blas_sgemm (op_A, op_B, alpha, A.gsl(), B.gsl(), beta, C.gsl()); }
 
-      inline void gemv (CBLAS_TRANSPOSE op_A, float alpha, const MatrixView<float>& A, const VectorView<float>& x, float beta, VectorView<float>& y)
-      { gsl_blas_sgemv (op_A, alpha, &A, &x, beta, &y); }
+      inline void gemv (CBLAS_TRANSPOSE op_A, float alpha, const Matrix<float>& A, const Vector<float>& x, float beta, Vector<float>& y)
+      { gsl_blas_sgemv (op_A, alpha, A.gsl(), x.gsl(), beta, y.gsl()); }
 
-      inline void symm (CBLAS_SIDE side, CBLAS_UPLO uplo, float alpha, const MatrixView<float>& A, const MatrixView<float>& B, float beta, MatrixView<float>& C)
-      { gsl_blas_ssymm (side, uplo, alpha, &A, &B, beta, &C); }
+      inline void symm (CBLAS_SIDE side, CBLAS_UPLO uplo, float alpha, const Matrix<float>& A, const Matrix<float>& B, float beta, Matrix<float>& C)
+      { gsl_blas_ssymm (side, uplo, alpha, A.gsl(), B.gsl(), beta, C.gsl()); }
 
-      inline void trsv (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, CBLAS_DIAG diag, const MatrixView<float>& A, VectorView<float>& x)
-      { gsl_blas_strsv (uplo, op_A, diag, &A, &x); }
+      inline void trsv (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, CBLAS_DIAG diag, const Matrix<float>& A, Vector<float>& x)
+      { gsl_blas_strsv (uplo, op_A, diag, A.gsl(), x.gsl()); }
 
-      inline void ger (float alpha, const VectorView<float>& x, const VectorView<float>& y, MatrixView<float>& A)
-      { gsl_blas_sger (alpha, &x, &y, &A); }
+      inline void ger (float alpha, const Vector<float>& x, const Vector<float>& y, Matrix<float>& A)
+      { gsl_blas_sger (alpha, x.gsl(), y.gsl(), A.gsl()); }
 
-      inline void syr (CBLAS_UPLO uplo, float alpha, const VectorView<float>& x, MatrixView<float>& A)
-      { gsl_blas_ssyr (uplo, alpha, &x, &A); }
+      inline void syr (CBLAS_UPLO uplo, float alpha, const Vector<float>& x, Matrix<float>& A)
+      { gsl_blas_ssyr (uplo, alpha, x.gsl(), A.gsl()); }
 
-      inline void syrk (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, float alpha, const MatrixView<float>& A, float beta, MatrixView<float>& C)
-      { gsl_blas_ssyrk (uplo, op_A, alpha, &A, beta, &C); }
+      inline void syrk (CBLAS_UPLO uplo, CBLAS_TRANSPOSE op_A, float alpha, const Matrix<float>& A, float beta, Matrix<float>& C)
+      { gsl_blas_ssyrk (uplo, op_A, alpha, A.gsl(), beta, C.gsl()); }
 
     }
 
diff --git a/lib/math/permutation.h b/lib/math/permutation.h
index ea6f73f493..e4703ce382 100644
--- a/lib/math/permutation.h
+++ b/lib/math/permutation.h
@@ -31,9 +31,9 @@
 namespace MR {
   namespace Math {
     namespace {
-      template <typename T> inline void permute (const gsl_permutation* p, VectorView<T>& V);
-      template <> inline void permute (const gsl_permutation* p, VectorView<double>& V) { gsl_permute_vector (p, &V); }
-      template <> inline void permute (const gsl_permutation* p, VectorView<float>& V) { gsl_permute_vector_float (p, &V); }
+      template <typename T> inline void permute (const gsl_permutation* p, Vector<T>& V);
+      template <> inline void permute (const gsl_permutation* p, Vector<double>& V) { gsl_permute_vector (p, V.gsl()); }
+      template <> inline void permute (const gsl_permutation* p, Vector<float>& V) { gsl_permute_vector_float (p, V.gsl()); }
     }
     
     class Permutation {
@@ -44,7 +44,7 @@ namespace MR {
         size_t operator[] (const size_t i) const { return (gsl_permutation_get (p, i)); }
         size_t& operator[] (const size_t i)       { return (gsl_permutation_data (p)[i]); }
 
-        template <typename T> VectorView<T>& apply (VectorView<T>& V) const { permute (p, V); return (V); }
+        template <typename T> Vector<T>& apply (Vector<T>& V) const { permute (p, V); return (V); }
         bool   valid () const { return (gsl_permutation_valid (p)); }
 
 	gsl_permutation* operator& () { return (p); }
diff --git a/lib/math/rician.h b/lib/math/rician.h
index d0ff52b34f..62f47fdc50 100644
--- a/lib/math/rician.h
+++ b/lib/math/rician.h
@@ -93,7 +93,7 @@ namespace MR {
 
 
 
-      template <typename T> inline T lnP (const VectorView<T>& measured, const VectorView<T>& actual, const T one_over_noise_squared, VectorView<T>& dP_dactual, T& dP_dN) 
+      template <typename T> inline T lnP (const Vector<T>& measured, const Vector<T>& actual, const T one_over_noise_squared, Vector<T>& dP_dactual, T& dP_dN) 
       {
         assert (one_over_noise_squared > 0.0);
         assert (measured.size() == actual.size());
diff --git a/lib/math/rng.h b/lib/math/rng.h
index f98f0f61de..35801ed62d 100644
--- a/lib/math/rng.h
+++ b/lib/math/rng.h
@@ -57,7 +57,7 @@ namespace MR {
           return (sqrt (amplitude*amplitude + imag*imag));
         }
 
-        template <typename T> void shuffle (VectorView<T>& V) { gsl_ran_shuffle (generator, (&V)->data, V.size(), sizeof (T)); }
+        template <typename T> void shuffle (Vector<T>& V) { gsl_ran_shuffle (generator, V->ptr(), V.size(), sizeof (T)); }
         template <class T> void shuffle (std::vector<T>& V) { gsl_ran_shuffle (generator, &V[0], V.size(), sizeof (T)); }
     };
 
diff --git a/lib/math/vector.h b/lib/math/vector.h
index eaeb55a63b..67607c2c5e 100644
--- a/lib/math/vector.h
+++ b/lib/math/vector.h
@@ -67,14 +67,39 @@ namespace MR {
     //! provides access to data as a vector
     /** \note this class is not capable of managing its own data allocation. 
       See the Vector class for a more general interface. */
-    template <typename T> class VectorView : protected GSLVector<T>
+    template <typename T> class Vector : protected GSLVector<T>
     {
       public:
-        template <typename U> friend class VectorView;
+        template <typename U> friend class Vector;
+
+	//! construct empty vector
+        Vector () throw () {
+	  GSLVector<T>::size = GSLVector<T>::stride = 1; 
+          GSLVector<T>::data = NULL; GSLVector<T>::block = NULL; GSLVector<T>::owner = 1;
+        }
+
+	//! copy constructor
+        /*! \note the new instance will reference the same data as the
+         * original, but will not be responsible for (or even capable of) any
+         * operation requiring data allocation or deallocation. */
+        Vector (const Vector& V) throw () {
+	  GSLVector<T>::size = V.GSLVector<T>::size;
+          GSLVector<T>::stride = V.GSLVector<T>::stride; 
+          GSLVector<T>::data = V.GSLVector<T>::data; 
+          GSLVector<T>::block = NULL; GSLVector<T>::owner = 0;
+        }
+
+	//! construct vector of size \a nelements
+	/** \note the elements of the vector are left uninitialised. */
+        Vector (size_t nelements) {   
+	  GSLVector<T>::block = GSLBlock<T>::alloc (nelements);
+          if (!GSLVector<T>::block) throw Exception ("Failed to allocate memory for Vector data");  
+	  GSLVector<T>::size = nelements; GSLVector<T>::stride = 1; 
+          GSLVector<T>::data = GSLVector<T>::block->data; GSLVector<T>::owner = 1;
+        }
 
 	//! construct from existing data array
-        explicit VectorView (T* vector_data, size_t nelements, size_t skip = 1)
-        { 
+        Vector (T* vector_data, size_t nelements, size_t skip = 1) { 
 	  GSLVector<T>::size = nelements;
 	  GSLVector<T>::stride = skip;
 	  GSLVector<T>::set (vector_data);
@@ -82,6 +107,73 @@ namespace MR {
 	  GSLVector<T>::owner = 0;
         }
 
+	//! construct a vector by reading from the text file \a filename
+        Vector (const std::string& file) { load (file); } 
+
+	//! destructor
+        ~Vector () { if (GSLVector<T>::block) { assert (GSLVector<T>::owner); GSLBlock<T>::free (GSLVector<T>::block); } }
+
+	//! deallocate the vector data
+        Vector<T>& clear () {
+          if (GSLVector<T>::block) {
+            assert (GSLVector<T>::owner);
+            GSLBlock<T>::free (GSLVector<T>::block);
+          }
+	  GSLVector<T>::size = GSLVector<T>::stride = 0; 
+          GSLVector<T>::data = NULL; 
+          GSLVector<T>::block = NULL; 
+          GSLVector<T>::owner = 1; 
+          return (*this); 
+        }
+
+	//! allocate the vector to have the same size as \a V
+        Vector& allocate (const Vector<T>& V) { return (allocate (V.size())); } 
+
+	//! allocate the vector to have the same size as \a V
+        template <typename U> Vector& allocate (const Vector<U>& V) { return (allocate (V.size())); }
+
+	//! allocate the vector to have size \a nelements
+        Vector& allocate (size_t nelements) {
+          if (nelements == size()) return (*this);
+          if (GSLVector<T>::block) {
+            assert (GSLVector<T>::owner);
+            if (GSLVector<T>::block->size < nelements) {
+	      GSLBlock<T>::free (GSLVector<T>::block);
+	      GSLVector<T>::block = NULL;
+            }
+          }
+          if (!GSLVector<T>::block && nelements) {
+            GSLVector<T>::block = GSLBlock<T>::alloc (nelements);
+            if (!GSLVector<T>::block) throw Exception ("Failed to allocate memory for Vector data");
+          }
+	  GSLVector<T>::size = nelements; GSLVector<T>::stride = 1; GSLVector<T>::owner = 1;
+          GSLVector<T>::data = GSLVector<T>::block ? GSLVector<T>::block->data : NULL;
+          return (*this);
+        }
+
+	//! resize the vector to have size \a nelements
+        /** \note no bounds checking and no data
+         * allocation/deallocation/copying is performed by this function. */
+        Vector& resize (size_t nelements) { GSLVector<T>::size = nelements; return (*this); }
+
+	//! resize the vector to refer to the subvector specified
+        /** \note no bounds checking and no data
+         * allocation/deallocation/copying is performed by this function. */
+        Vector& resize (size_t from, size_t to) { 
+          ptr() += from*GSLVector<T>::stride; 
+          GSLVector<T>::size = to-from; 
+          return (*this); 
+        }
+
+	//! read vector data from the text file \a filename
+        Vector<T>& load (const std::string& filename) { 
+          std::ifstream in (filename.c_str()); 
+          if (!in) throw Exception ("cannot open matrix file \"" + filename + "\": " + strerror (errno)); 
+          try { in >> *this; }
+          catch (Exception& E) { throw Exception (E, "error loading matrix file \"" + filename + "\""); }
+          return (*this);
+        }
+
 	//! write to text file \a filename
         void save (const std::string& filename) const {   
           std::ofstream out (filename.c_str()); 
@@ -90,225 +182,109 @@ namespace MR {
         }
 
 	//! used to obtain a pointer to the underlying GSL structure
-	GSLVector<T>* operator& () { return (this); }
+	GSLVector<T>* gsl () { return (this); }
 	//! used to obtain a pointer to the underlying GSL structure
-        const GSLVector<T>* operator& () const { return (this); }
+        const GSLVector<T>* gsl () const { return (this); }
 
 	//! true if vector points to existing data
-        bool       is_set () const throw () { return (ptr()); }
+        bool is_set () const throw () { return (ptr()); }
 	//! returns number of elements of vector
-        size_t     size () const throw ()  { return (GSLVector<T>::size); }
+        size_t size () const throw ()  { return (GSLVector<T>::size); }
 
 	//! returns a reference to the element at \a i
-        T&        operator[] (size_t i) throw ()          { return (ptr()[i*GSLVector<T>::stride]); }
+        T& operator[] (size_t i) throw ()          { return (ptr()[i*GSLVector<T>::stride]); }
+
 	//! returns a reference to the element at \a i
-        const T&  operator[] (size_t i) const throw ()    { return (ptr()[i*GSLVector<T>::stride]); }
+        const T& operator[] (size_t i) const throw ()    { return (ptr()[i*GSLVector<T>::stride]); }
+
+	//! return a pointer to the underlying data
+        T* ptr () throw () { return ((T*) (GSLVector<T>::data)); }
+
+	//! return a pointer to the underlying data
+        const T* ptr () const throw () { return ((const T*) (GSLVector<T>::data)); }
+
+	//! return the stride of the vector
+        size_t stride () const throw () { return (GSLVector<T>::stride); }
 
 	//! assign the specified \a value to all elements of the vector
-        VectorView& operator= (T value) throw () { LOOP (operator[](i) = value); return (*this); }   
+        Vector& operator= (T value) throw () { LOOP (operator[](i) = value); return (*this); }   
+
 	//! assign the values in \a V to the corresponding elements of the vector
-        VectorView& operator= (const VectorView& V) { assert (GSLVector<T>::size == V.GSLVector<T>::size); LOOP (operator[](i) = V[i]); return (*this); }
+        Vector& operator= (const Vector& V) { 
+          assert (GSLVector<T>::size == V.GSLVector<T>::size); 
+          LOOP (operator[](i) = V[i]); 
+          return (*this);
+        }
+
 	//! assign the values in \a V to the corresponding elements of the vector
-        template <typename U> VectorView& operator= (const VectorView<U>& V) { assert (GSLVector<T>::size == V.GSLVector<T>::size); LOOP (operator[](i) = V[i]); return (*this); }
+        template <typename U> Vector& operator= (const Vector<U>& V) { 
+          assert (GSLVector<T>::size == V.GSLVector<T>::size); 
+          LOOP (operator[](i) = V[i]); 
+          return (*this);
+        }
 
 
 	//! set all elements of vector to zero
-        VectorView& zero () throw () { LOOP (operator[](i) = 0.0); return (*this); }
+        Vector& zero () throw () { LOOP (operator[](i) = 0.0); return (*this); }
+
 	//! swap contents with \a V without copying
-        void swap (VectorView& V) throw () { 
-          char c [sizeof (VectorView)];
-          memcpy (&c, this, sizeof (VectorView));
-          memcpy (this, &V, sizeof (VectorView));
-          memcpy (&V, &c, sizeof (VectorView));
+        void swap (Vector& V) throw () { 
+          char c [sizeof (Vector)];
+          memcpy (&c, this, sizeof (Vector));
+          memcpy (this, &V, sizeof (Vector));
+          memcpy (&V, &c, sizeof (Vector));
         } 
 
 	//! add \a value to all elements of the vector
-        VectorView& operator+= (T value) throw () { LOOP (operator[](i) += value); return (*this); } 
+        Vector& operator+= (T value) throw () { LOOP (operator[](i) += value); return (*this); } 
 	//! subtract \a value from all elements of the vector
-        VectorView& operator-= (T value) throw () { LOOP (operator[](i) -= value); return (*this); }   
+        Vector& operator-= (T value) throw () { LOOP (operator[](i) -= value); return (*this); }   
 	//! multiply all elements of the vector by \a value
-        VectorView& operator*= (T value) throw () { LOOP (operator[](i) *= value); return (*this); }   
+        Vector& operator*= (T value) throw () { LOOP (operator[](i) *= value); return (*this); }   
 	//! divide all elements of the vector by \a value
-        VectorView& operator/= (T value) throw () { LOOP (operator[](i) /= value); return (*this); }   
+        Vector& operator/= (T value) throw () { LOOP (operator[](i) /= value); return (*this); }   
 
 	//! add each element of \a V to the corresponding element of the vector
-        VectorView& operator+= (const VectorView& V) throw () { LOOP (operator[](i) += V[i]); return (*this); }   
+        Vector& operator+= (const Vector& V) throw () { LOOP (operator[](i) += V[i]); return (*this); }   
 	//! subtract each element of \a V from the corresponding element of the vector
-        VectorView& operator-= (const VectorView& V) throw () { LOOP (operator[](i) -= V[i]); return (*this); }
+        Vector& operator-= (const Vector& V) throw () { LOOP (operator[](i) -= V[i]); return (*this); }
 	//! multiply each element of \a V by the corresponding element of the vector
-        VectorView& operator*= (const VectorView& V) throw () { LOOP (operator[](i) *= V[i]); return (*this); }
+        Vector& operator*= (const Vector& V) throw () { LOOP (operator[](i) *= V[i]); return (*this); }
 	//! divide each element of \a V by the corresponding element of the vector
-        VectorView& operator/= (const VectorView& V) throw () { LOOP (operator[](i) /= V[i]); return (*this); }
+        Vector& operator/= (const Vector& V) throw () { LOOP (operator[](i) /= V[i]); return (*this); }
 
-	//! return a VectorView corresponding to a subvector of the vector
-        VectorView view (size_t from, size_t to) throw () { 
+	//! return a subvector of the vector
+        Vector sub (size_t from, size_t to) throw () { 
           assert (from <= to && to <= size());
-          return (VectorView (ptr() + from*GSLVector<T>::stride, to-from, GSLVector<T>::stride));
+          return (Vector (ptr() + from*GSLVector<T>::stride, to-from, GSLVector<T>::stride));
         }
 
-	//! return a VectorView corresponding to a subvector of the vector
-        const VectorView view (size_t from, size_t to) const throw () { 
+	//! return a subvector of the vector
+        const Vector sub (size_t from, size_t to) const throw () { 
           assert (from <= to && to <= size());
-          return (VectorView (ptr() + from*GSLVector<T>::stride, to-from, GSLVector<T>::stride));
+          return (Vector (ptr() + from*GSLVector<T>::stride, to-from, GSLVector<T>::stride));
         }
 
-	//! return a VectorView corresponding to a subvector of the vector
-        VectorView view (size_t from, size_t to, size_t skip) throw () { 
+	//! return a subvector of the vector
+        Vector sub (size_t from, size_t to, size_t skip) throw () { 
           assert (from <= to && to <= size());
-          return (VectorView (ptr() + from*GSLVector<T>::stride, ceil<size_t> ((to-from)/float(skip)), GSLVector<T>::stride*skip));
+          return (Vector (ptr() + from*GSLVector<T>::stride, ceil<size_t> ((to-from)/float(skip)), GSLVector<T>::stride*skip));
         }
 
-	//! return a VectorView corresponding to a subvector of the vector
-        const VectorView view (size_t from, size_t to, size_t skip) const throw () { 
+	//! return a subvector of the vector
+        const Vector sub (size_t from, size_t to, size_t skip) const throw () { 
           assert (from <= to && to <= size());
-          return (VectorView (ptr() + from*GSLVector<T>::stride, ceil<size_t> ((to-from)/float(skip)), GSLVector<T>::stride*skip));
+          return (Vector (ptr() + from*GSLVector<T>::stride, ceil<size_t> ((to-from)/float(skip)), GSLVector<T>::stride*skip));
         }
 
-	//! resize the vector to have size \a nelements
-	/** \note no bounds checking and no data allocation/deallocation/copying is performed by this function. */
-        VectorView& resize (size_t nelements) { GSLVector<T>::size = nelements; return (*this); }
-	//! resize the vector to refer to the subvector specified
-	/** \note no bounds checking and no data allocation/deallocation/copying is performed by this function. */
-        VectorView& resize (size_t from, size_t to) { ptr() += from*GSLVector<T>::stride; GSLVector<T>::size = to-from; return (*this); }
-
-	//! return a pointer to the underlying data
-        T* ptr () throw () { return ((T*) (GSLVector<T>::data)); }
-	//! return a pointer to the underlying data
-        const T* ptr () const throw () { return ((const T*) (GSLVector<T>::data)); }
-	//! return the stride of the vector
-        size_t stride () const throw () { return (GSLVector<T>::stride); }
-
 	//! write the vector \a V to \a stream as text
-        friend std::ostream& operator<< (std::ostream& stream, const VectorView& V)  {
+        friend std::ostream& operator<< (std::ostream& stream, const Vector& V)  {
           for (size_t i = 0; i < V.size(); i++) stream << V[i] << "\n"; 
           return (stream); 
         }
 
-      protected:
-        VectorView () { 
-	  GSLVector<T>::size = 0;
-	  GSLVector<T>::stride = 0;
-	  GSLVector<T>::data = NULL;
-	  GSLVector<T>::block = NULL;
-	  GSLVector<T>::owner = 0;
-        }
-
-    };
-
-
-    //! The main vector class 
-    template <typename T> class Vector : public VectorView<T> { 
-      public:
-        template <typename U> friend class Vector; 
-
-	//! construct empty vector
-        explicit Vector () throw () { } 
-	//! construct vector of size \a nelements
-	/** \note the elements of the vector are left uninitialised. */
-        explicit Vector (size_t nelements) {   
-	  GSLVector<T>::block = GSLBlock<T>::alloc (nelements);
-          if (!GSLVector<T>::block) throw Exception ("Failed to allocate memory for Vector data");  
-	  GSLVector<T>::size = nelements; GSLVector<T>::stride = 1; GSLVector<T>::data = GSLVector<T>::block->data; GSLVector<T>::owner = 1;
-        }
-	//! construct a vector by copying the elements of \a V
-        explicit Vector (const Vector& V) {   
-          if (V.size()) {
-	    GSLVector<T>::block = GSLBlock<T>::alloc (V.size());
-            if (!GSLVector<T>::block) throw Exception ("Failed to allocate memory for Vector data");  
-	    GSLVector<T>::data = GSLVector<T>::block->data; GSLVector<T>::size = V.size(); GSLVector<T>::stride = 1; GSLVector<T>::owner = 1;
-            view() = V.view();
-          }
-        }
-	//! construct a vector by copying the elements of \a V
-        template <typename U> Vector (const Vector<U>& V) {
-          if (V.size()) {
-	    GSLVector<T>::block = GSLBlock<T>::alloc (V.size());
-            if (!GSLVector<T>::block) throw Exception ("Failed to allocate memory for Vector data");  
-	    GSLVector<T>::data = GSLVector<T>::block->data; GSLVector<T>::size = V.size(); GSLVector<T>::stride = 1; GSLVector<T>::owner = 1;
-            view() = V.view();
-          }
-        }
-
-	//! construct a vector by reading from the text file \a filename
-        explicit Vector (const std::string& file) { load (file); } 
-	//! destructor
-        ~Vector () { delete [] GSLVector<T>::block; } 
-
-	//! assignment operator: allocate vector to have the same size as \a V, and copy contents from \a V
-        Vector&   operator= (const Vector<T>& V) { return (operator= (V.view())); }
-	//! assignment operator: allocate vector to have the same size as \a V, and copy contents from \a V
-        Vector&   operator= (const VectorView<T>& V) { allocate (V); view() = V; return (*this); }
-	//! assignment operator: allocate vector to have the same size as \a V, and copy contents from \a V
-        template <typename U> Vector& operator= (const Vector<U>& V) { return (operator= (V.view())); }
-	//! assignment operator: allocate vector to have the same size as \a V, and copy contents from \a V
-        template <typename U> Vector& operator= (const VectorView<U>& V) { allocate (V); view() = V; return (*this); }
-
-	//! allocate the vector to have the same size as \a V
-        Vector& allocate (const VectorView<T>& V) { return (allocate (V.size())); } 
-	//! allocate the vector to have the same size as \a V
-        template <typename U> Vector& allocate (const Vector<U>& V) { return (allocate (V.size())); }
-	//! allocate the vector to have size \a nelements
-        Vector& allocate (size_t nelements) {
-          if (GSLVector<T>::block) {
-            if (GSLVector<T>::block->size < nelements) {
-	      GSLBlock<T>::free (GSLVector<T>::block);
-	      GSLVector<T>::block = NULL;
-            }
-          }
-          if (!GSLVector<T>::block && nelements) {
-            GSLVector<T>::block = GSLBlock<T>::alloc (nelements);
-            if (!GSLVector<T>::block) throw Exception ("Failed to allocate memory for Vector data");
-          }
-	  GSLVector<T>::size = nelements; GSLVector<T>::stride = 1;
-          GSLVector<T>::data = GSLVector<T>::block ? GSLVector<T>::block->data : NULL;
-          return (*this);
-        }
-
-	//! resize vector to have size \a nelements
-        Vector<T>& resize (size_t nelements) {
-          if (nelements <= GSLVector<T>::size) { GSLVector<T>::size = nelements; return (*this); }
-          if (GSLVector<T>::block)
-            if ((nelements-1)*GSLVector<T>::stride + 1 <= GSLVector<T>::block->size) { GSLVector<T>::size = nelements; return (*this); }
-          Vector<T> V (nelements);
-          size_t n = MIN(GSLVector<T>::size, nelements);
-          V.VectorView<T>::view(0,n) = VectorView<T>::view(0,n);
-          clear();
-          VectorView<T>::swap (V);
-          return (*this);
-        }
-
-	//! deallocate the vector data
-        Vector<T>& clear () {
-          delete [] GSLVector<T>::block;
-	  GSLVector<T>::size = GSLVector<T>::stride = 0; GSLVector<T>::data = NULL; GSLVector<T>::block = NULL; GSLVector<T>::owner = 0; return (*this); 
-        }
-
-	//! return as VectorView class
-        VectorView<T>& view () { return (*this); }
-	//! return as VectorView class
-        const VectorView<T>& view () const { return (*this); }
-
-	//! return a VectorView corresponding to a subvector of the vector
-        VectorView<T> view (size_t from, size_t to, size_t skip = 1) throw () { return (VectorView<T>::view (from, to, skip)); }
-	//! return a VectorView corresponding to a subvector of the vector
-        const VectorView<T> view (size_t from, size_t to, size_t skip = 1) const throw () { return (VectorView<T>::view (from, to, skip)); }
-
-	//! read vector data from the text file \a filename
-        Vector<T>& load (const std::string& filename) { 
-          std::ifstream in (filename.c_str()); 
-          if (!in) throw Exception ("cannot open matrix file \"" + filename + "\": " + strerror (errno)); 
-          try { 
-            in >> *this; 
-          }
-          catch (Exception& E) { 
-            throw Exception (E, "error loading matrix file \"" + filename + "\""); 
-          }
-          return (*this);
-        }
-
 	//! read the vector data from \a stream and assign to the vector \a V
-        friend std::istream& operator>> (std::istream& stream, Vector<T>& V)  
-        {
+        friend std::istream& operator>> (std::istream& stream, Vector<T>& V)  {
           std::vector<T> vec;
           while (true) {
             T val;
@@ -321,8 +297,8 @@ namespace MR {
           for (size_t n = 0; n < V.size(); n++) V[n] = vec[n];
           return (stream);
         }
-    };
 
+    };
 
 
 
@@ -337,13 +313,13 @@ namespace MR {
     }
 
     //! compute the squared 2-norm of a vector
-    template <typename T> inline T norm2 (const VectorView<T>& V) { return (norm2 (V.ptr(), V.size(), V.stride())); }
+    template <typename T> inline T norm2 (const Vector<T>& V) { return (norm2 (V.ptr(), V.size(), V.stride())); }
 
     //! compute the 2-norm of a vector
     template <typename T> inline T norm (const T* V, size_t size = 3, size_t stride = 1) { return (sqrt(norm2(V, size, stride))); }
 
     //! compute the 2-norm of a vector
-    template <typename T> inline T norm (const VectorView<T>& V) { return (norm(V.ptr(), V.size(), V.stride())); }
+    template <typename T> inline T norm (const Vector<T>& V) { return (norm(V.ptr(), V.size(), V.stride())); }
 
     //! compute the squared 2-norm of the difference between two vectors
     template <typename T> inline T norm_diff2 (const T* x, const T* y, size_t size = 3, size_t x_stride = 1, size_t y_stride = 1) 
@@ -354,7 +330,7 @@ namespace MR {
     }
 
     //! compute the squared 2-norm of the difference between two vectors
-    template <typename T> inline T norm_diff2 (const VectorView<T>& x, const VectorView<T>& y) {
+    template <typename T> inline T norm_diff2 (const Vector<T>& x, const Vector<T>& y) {
       return (norm_diff2 (x.ptr(), y.ptr(), x.size(), x.stride(), y.stride()));
     }
 
@@ -366,7 +342,7 @@ namespace MR {
     }
 
     //! compute the mean of the elements of a vector
-    template <typename T> inline T mean (const VectorView<T>& V) { return (mean(V.ptr(), V.size(), V.stride())); }
+    template <typename T> inline T mean (const Vector<T>& V) { return (mean(V.ptr(), V.size(), V.stride())); }
 
     //! normalise a vector to have unit 2-norm
     template <typename T> inline void normalise (T* V, size_t size = 3, size_t stride = 1) {
@@ -375,7 +351,10 @@ namespace MR {
     }
 
     //! normalise a vector to have unit 2-norm
-    template <typename T> inline VectorView<T> normalise (VectorView<T>& V) { normalise (V.ptr(), V.size(), V.stride()); return (V); }
+    template <typename T> inline Vector<T>& normalise (Vector<T>& V) {
+      normalise (V.ptr(), V.size(), V.stride());
+      return (V); 
+    }
 
     //! compute the dot product between two vectors
     template <typename T> inline T dot (const T* x, const T* y, size_t size = 3, size_t x_stride = 1, size_t y_stride = 1) 
@@ -386,7 +365,9 @@ namespace MR {
     }
 
     //! compute the dot product between two vectors
-    template <typename T> inline T dot (const VectorView<T>& x, const VectorView<T>& y) { return (dot (x.ptr(), y.ptr(), x.size(), x.stride(), y.stride())); }
+    template <typename T> inline T dot (const Vector<T>& x, const Vector<T>& y) {
+      return (dot (x.ptr(), y.ptr(), x.size(), x.stride(), y.stride())); 
+    }
 
     //! compute the cross product between two vectors
     template <typename T> inline void cross (T* c, const T* x, const T* y, size_t c_stride = 1, size_t x_stride = 1, size_t y_stride = 1) 
@@ -397,14 +378,14 @@ namespace MR {
     }
 
     //! compute the cross product between two vectors
-    template <typename T> inline VectorView<T> cross (VectorView<T> c, const VectorView<T>& x, const VectorView<T>& y) { 
+    template <typename T> inline Vector<T> cross (Vector<T> c, const Vector<T>& x, const Vector<T>& y) { 
       cross (c.ptr(), x.ptr(), y.ptr(), c.stride(), x.stride(), y.stride()); 
       return (c); 
     }
 
 
     //! find the maximum value of any elements within a vector
-    template <typename T> inline T max (const VectorView<T>& V, size_t& i) 
+    template <typename T> inline T max (const Vector<T>& V, size_t& i) 
     {
       T val (V[0]);
       i = 0;
@@ -415,7 +396,7 @@ namespace MR {
     }
 
     //! find the minimum value of any elements within a vector
-    template <typename T> inline T min (const VectorView<T>& V, size_t& i) 
+    template <typename T> inline T min (const Vector<T>& V, size_t& i) 
     {
       T val (V[0]);
       i = 0;
@@ -427,7 +408,7 @@ namespace MR {
 
 
     //! find the maximum absolute value of any elements within a vector
-    template <typename T> inline T absmax (const VectorView<T>& V, size_t& i) 
+    template <typename T> inline T absmax (const Vector<T>& V, size_t& i) 
     {
       T val (abs(V[0]));
       i = 0;
diff --git a/lib/mrtrix.cpp b/lib/mrtrix.cpp
index ae34b666ba..0ca6d3942a 100644
--- a/lib/mrtrix.cpp
+++ b/lib/mrtrix.cpp
@@ -64,30 +64,30 @@ namespace MR {
 
 
 
-  std::vector<ssize_t> parse_ints (const std::string& spec, ssize_t last)
+  std::vector<int> parse_ints (const std::string& spec, int last)
   {
-    std::vector<ssize_t> V;
+    std::vector<int> V;
     if (!spec.size()) throw Exception ("integer sequence specifier is empty"); 
     std::string::size_type start = 0, end;
-    ssize_t num[3];
-    ssize_t i = 0;
+    int num[3];
+    int i = 0;
     try {
       do {
         end = spec.find_first_of (",:", start);
         std::string str (strip (spec.substr (start, end-start)));
         lowercase (str);
         if (str == "end") {
-          if (last == std::numeric_limits<ssize_t>::max()) 
+          if (last == std::numeric_limits<int>::max()) 
             throw Exception ("value of \"end\" is not known in number sequence \"" + spec + "\"");
           num[i] = last;
         }
-        else num[i] = to<ssize_t> (spec.substr (start, end-start));
+        else num[i] = to<int> (spec.substr (start, end-start));
 
         char last_char = end < spec.size() ? spec[end] : '\0';
         if (last_char == ':') { i++; if (i > 2) throw Exception ("invalid number range in number sequence \"" + spec + "\""); }
         else {
           if (i) {
-            ssize_t inc, last;
+            int inc, last;
             if (i == 2) { inc = num[1]; last = num[2]; }
             else { inc = 1; last = num[1]; }
             if (inc * (last - num[0]) < 0) inc = -inc;
diff --git a/lib/mrtrix.h b/lib/mrtrix.h
index 3e4e52832d..ab4cad0c7d 100644
--- a/lib/mrtrix.h
+++ b/lib/mrtrix.h
@@ -46,7 +46,7 @@
 
 #define MRTRIX_MAJOR_VERSION 0
 #define MRTRIX_MINOR_VERSION 3
-#define MRTRIX_MICRO_VERSION 2
+#define MRTRIX_MICRO_VERSION 3
 
 
 /** Prints the file and line number. Useful for debugging purposes. */
@@ -217,7 +217,7 @@ namespace MR {
   }
 
   std::vector<float> parse_floats (const std::string& spec);
-  std::vector<ssize_t>   parse_ints (const std::string& spec, ssize_t last = std::numeric_limits<ssize_t>::max());
+  std::vector<int>   parse_ints (const std::string& spec, int last = std::numeric_limits<int>::max());
 
   inline int round (float x) { return (int (x + (x > 0.0 ? 0.5 : -0.5))); }
 
diff --git a/src/dwi/sdeconv/constrained.h b/src/dwi/sdeconv/constrained.h
index dd66bdd819..a47ab4f1ee 100644
--- a/src/dwi/sdeconv/constrained.h
+++ b/src/dwi/sdeconv/constrained.h
@@ -81,8 +81,8 @@ namespace MR {
                 HR_trans *= neg_lambda * T(fconv.rows()) * response[0] / T(HR_trans.rows());
 
                 M.allocate (DW_dirs.rows(), HR_trans.columns());
-                M.view (0,M.rows(),0,fconv.columns()) = fconv.view();
-                M.view (0,M.rows(),fconv.columns(),M.columns()) = 0.0;
+                M.sub (0,M.rows(),0,fconv.columns()) = fconv;
+                M.sub (0,M.rows(),fconv.columns(),M.columns()) = 0.0;
                 Mt_M.allocate (M.columns(), M.columns());
                 rankN_update (Mt_M, M, CblasTrans);
 
@@ -106,13 +106,25 @@ namespace MR {
             norm_lambda = NORM_LAMBDA_MULTIPLIER * P.norm_lambda * P.Mt_M(0,0);
           }
 
+        CSDeconv (const CSDeconv& c) :
+          P (c.P), 
+          work (P.Mt_M.rows(), P.Mt_M.columns()),
+          HR_T (P.HR_trans.rows(), P.HR_trans.columns()), 
+          F (P.HR_trans.columns()),
+          init_F (P.rconv.rows()),
+          HR_amps (P.HR_trans.rows()),
+          Mt_b (P.HR_trans.columns()),
+          old_neg (P.HR_trans.rows()) { 
+            norm_lambda = NORM_LAMBDA_MULTIPLIER * P.norm_lambda * P.Mt_M(0,0);
+          }
+
         ~CSDeconv() { }
 
         void set (const Math::Vector<T>& DW_signals)
         {
           Math::mult (init_F, P.rconv, DW_signals);
-          F.view (0, init_F.size()) = init_F.view();
-          F.view (init_F.size(), F.size()) = 0.0;
+          F.sub (0, init_F.size()) = init_F;
+          F.sub (init_F.size(), F.size()) = 0.0;
           old_neg.assign (P.HR_trans.rows(), -1);
           threshold = P.threshold * init_F[0] * P.HR_trans(0,0);
           computed_once = false;
@@ -144,7 +156,7 @@ namespace MR {
             rankN_update (work, HR_T, CblasTrans, CblasLower, T(1.0), T(1.0));
           }
 
-          F.view() = Mt_b.view();
+          F = Mt_b;
 
           Math::Cholesky::decomp (work); 
           Math::Cholesky::solve (F, work);