init

.gitignore

+*build
+cpubuild/
+cudabuild/

BiasFieldEstimator.h

+/*! \file BiasFieldEstimator.h
+    \brief Contains declaration of class for estimation of a bias field
+
+    \author Jesper Andersson
+    \version 1.0b, April, 2017.
+*/
+// 
+// BiasFieldEstimator.h
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2017 University of Oxford 
+//
+
+#ifndef BiasFieldEstimator_h
+#define BiasFieldEstimator_h
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <memory>
+#include <time.h>
+#include "newimage/newimageall.h"
+#include "ECScanClasses.h"
+
+namespace EDDY {
+
+class BiasFieldEstimatorImpl;
+
+/****************************************************************//**
+*
+* \brief Class used to estimate a receive bias-field.
+*
+* Class used to estimate a receive bias-field based on data
+* from a subject that moves around within that bias-field. 
+* It is implemented using the "Pimpl idiom" which means that this class 
+* only implements an interface whereas the actual work is being performed
+* by the BiasFieldEstimatorImpl class which is declared and defined in
+* BiasFieldEstimatorImpl.cpp or cuda/BiasFieldEstimatorImpl.cu depending on 
+* what platform the code is compiled for.
+*
+* The way you would use it is by first constructing the object and
+* then calling `AddScan()` for all volumes and finally call `GetField()`
+* to retrieve the estimated field
+*
+********************************************************************/ 
+class BiasFieldEstimator
+{
+public:
+  BiasFieldEstimator();
+  ~BiasFieldEstimator();
+  /// Set ref scan
+  void SetRefScan(const NEWIMAGE::volume<float>& mask, const EDDY::ImageCoordinates& coords);
+  /// Add another scan
+  void AddScan(const NEWIMAGE::volume<float>& predicted, const NEWIMAGE::volume<float>& observed, 
+	       const NEWIMAGE::volume<float>& mask, const EDDY::ImageCoordinates& coords);
+  /// Calculate and return "direct" representation of the bias-field
+  NEWIMAGE::volume<float> GetField(double lambda) const;
+  /// Calculate and return spline basis-field representation of the bias-field
+  NEWIMAGE::volume<float> GetField(double ksp, double lambda) const;
+  /// Caclulate and return At matrix for debug purposes
+  MISCMATHS::SpMat<float> GetAtMatrix(const EDDY::ImageCoordinates&  coords, 
+				      const NEWIMAGE::volume<float>& predicted, 
+				      const NEWIMAGE::volume<float>& mask) const;
+  /// Calculate and return At* matrix for debug purposes
+  MISCMATHS::SpMat<float> GetAtStarMatrix(const EDDY::ImageCoordinates&  coords, const NEWIMAGE::volume<float>& mask) const;
+  ///
+  NEWIMAGE::volume<float> GetATimesField(const EDDY::ImageCoordinates&  coords, 
+					 const NEWIMAGE::volume<float>& predicted, 
+					 const NEWIMAGE::volume<float>& mask,
+					 const NEWIMAGE::volume<float>& field); 
+  /// Write out current state for debug purposes
+  void Write(const std::string& basename) const;
+private:
+  BiasFieldEstimatorImpl* _pimpl;
+};
+
+} // End namespace EDDY
+#endif // End #ifndef BiasFieldEstimator_h

BiasFieldEstimatorImpl.cpp

+/*! \file BiasFieldEstimatorImpl.cpp
+    \brief Contains one implementation of class for estimation of a bias field
+
+    \author Jesper Andersson
+    \version 1.0b, December, 2017.
+*/
+//
+// BiasFieldEstimatorImpl.h
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2017 University of Oxford
+//
+
+#ifndef BiasFieldEstimatorImpl_h
+#define BiasFieldEstimatorImpl_h
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <memory>
+#include <time.h>
+#include <boost/shared_ptr.hpp>
+#include "armawrap/newmat.h"
+#include "basisfield/basisfield.h"
+#include "basisfield/splinefield.h"
+#include "newimage/newimageall.h"
+#include "miscmaths/SpMat.h"
+#include "miscmaths/SpMatMatrices.h"
+#include "EddyHelperClasses.h"
+#include "BiasFieldEstimator.h"
+
+namespace EDDY {
+
+class BiasFieldEstimatorImpl
+{
+public:
+  BiasFieldEstimatorImpl() : _nima(0) {}
+  /// Set ref scan
+  void SetRefScan(const NEWIMAGE::volume<float>& mask, const EDDY::ImageCoordinates& coords);
+  /// Add scan
+  void AddScan(const NEWIMAGE::volume<float>& predicted, const NEWIMAGE::volume<float>& observed,
+	       const NEWIMAGE::volume<float>& mask, const EDDY::ImageCoordinates& coords);
+  /// Calculate and return "direct" representation of the bias-field
+  NEWIMAGE::volume<float> GetField(double lambda) const;
+  /// Calculate and return spline basis-field representation of the bias-field
+  NEWIMAGE::volume<float> GetField(double ksp, double lambda) const;
+  /// Caclulate and return At matrix, for debug purposes
+  MISCMATHS::SpMat<float> GetAtMatrix(const EDDY::ImageCoordinates&  coords,
+				      const NEWIMAGE::volume<float>& predicted,
+				      const NEWIMAGE::volume<float>& mask) const;
+  /// Calculate and return At* matrix, for debug purposes
+  MISCMATHS::SpMat<float> GetAtStarMatrix(const EDDY::ImageCoordinates&  coords, const NEWIMAGE::volume<float>& mask) const;
+  /// Calculate and return A*b, where b is the field, for debug purposes
+  NEWIMAGE::volume<float> GetATimesField(const EDDY::ImageCoordinates&  coords,
+					 const NEWIMAGE::volume<float>& predicted,
+					 const NEWIMAGE::volume<float>& mask,
+					 const NEWIMAGE::volume<float>& field) const;
+  /// Write out current state for debug purposes
+  void Write(const std::string& basename) const EddyTry {
+    MISCMATHS::write_ascii_matrix(_Aty,"BiasFieldEstimatorImpl_" + basename + "_Aty.txt");
+    _AtA.Print("BiasFieldEstimatorImpl_" + basename + "_AtA.txt");
+  } EddyCatch
+private:
+  unsigned int              _nima;    /// Number of loaded images
+  std::vector<unsigned int> _isz;     /// Image matrix size
+  std::vector<double>       _vxs;     /// Voxel size
+  MISCMATHS::SpMat<float>   _rAt;     /// A^T matrix for reference volume
+  MISCMATHS::SpMat<float>   _AtA;     /// \Sum A^T A in paper
+  NEWMAT::ColumnVector      _Aty;     /// \Sum A^t f-\hat{f} in paper
+
+  struct PointersForSpMat {
+    PointersForSpMat(unsigned int N) : nz(0),
+				       irp(std::unique_ptr<unsigned int[]>(new unsigned int[8*N])),
+				       jcp(std::unique_ptr<unsigned int[]>(new unsigned int[N+1])),
+				       sp(std::unique_ptr<double[]>(new double[8*N])) {}
+    PointersForSpMat(PointersForSpMat&& rhs) : nz(rhs.nz),              // Moving is fine
+					       irp(std::move(rhs.irp)),
+					       jcp(std::move(rhs.jcp)),
+					       sp(std::move(rhs.sp)) {}
+    PointersForSpMat(const PointersForSpMat& rhs) = delete;             // No copying
+    PointersForSpMat& operator=(const PointersForSpMat& rhs) = delete;  // No copying
+    unsigned int                    nz;
+    std::unique_ptr<unsigned int[]> irp;
+    std::unique_ptr<unsigned int[]> jcp;
+    std::unique_ptr<double[]>       sp;
+  };
+
+  unsigned int nvox() const { return(_isz[0]*_isz[1]*_isz[2]); }
+  unsigned int ijk2indx(unsigned int i, unsigned int j, unsigned int k, const EDDY::ImageCoordinates& c) const {
+    return(k*c.NX()*c.NY() + j*c.NX() + i);
+  }
+  void get_wgts(const double& dx, const double& dy, const double& dz,
+		double& w000, double& w100, double& w010, double& w110,
+		double& w001, double& w101, double& w011, double& w111) const
+  {
+    w000 = w100 = w010 = w110 = 1.0 - dz;
+    w001 = w101 = w011 = w111 = dz;
+    w000 *= 1.0 - dy; w100 *= 1.0 - dy; w001 *= 1.0 - dy; w101 *= 1.0 - dy;
+    w010 *= dy; w110 *= dy; w011 *= dy; w111 *= dy;
+    w000 *= 1.0 - dx; w010 *= 1.0 - dx; w001 *= 1.0 - dx; w011 *= 1.0 - dx;
+    w100 *= dx; w110 *= dx; w101 *= dx; w111 *= dx;
+  }
+  PointersForSpMat make_At_star_CSC(// Input
+				    const EDDY::ImageCoordinates&  coords,
+				    const NEWIMAGE::volume<float>& mask) const;
+  void multiply_At_star_CSC_by_image(const NEWIMAGE::volume<float>& ima,
+				     const NEWIMAGE::volume<float>& mask,
+				     unsigned int* const            jcp,
+				     double* const                  sp) const;
+
+};
+
+BiasFieldEstimatorImpl::PointersForSpMat BiasFieldEstimatorImpl::make_At_star_CSC(// Input
+										  const EDDY::ImageCoordinates&  coords,
+										  const NEWIMAGE::volume<float>& mask) const EddyTry
+{
+  // Make A^T in a Compressed Column Storage representation
+  PointersForSpMat ptrs(coords.N());
+  unsigned int ii = 0; // Index intp irp and sp
+  unsigned int ji = 0; // Index into jcp
+  NEWMAT::ColumnVector vmask = mask.vec();
+  for (unsigned int i=0; i<coords.N(); i++) {
+    ptrs.jcp[ji++] = ii;
+    if (vmask(i+1) > 0.0 && coords.IsInBounds(i)) { // If voxel falls within volume
+      unsigned int xi = floor(coords.x(i));
+      unsigned int yi = floor(coords.y(i));
+      unsigned int zi = floor(coords.z(i));
+      double dx = coords.x(i) - xi;
+      double dy = coords.y(i) - yi;
+      double dz = coords.z(i) - zi;
+      if (dx < 1e-6 && dy < 1e-6 && dz < 1e-6) { // Voxel falls "exactly" on original voxel centre
+	ptrs.irp[ii] = this->ijk2indx(xi,yi,zi,coords);
+	ptrs.sp[ii++] = 1.0;
+      }
+      else {
+	double w000, w100, w010, w110, w001, w101, w011, w111;
+	get_wgts(dx,dy,dz,w000,w100,w010,w110,w001,w101,w011,w111);
+	if (w000 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi,yi,zi,coords); ptrs.sp[ii++] = w000; }
+	if (w100 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi+1,yi,zi,coords); ptrs.sp[ii++] = w100; }
+	if (w010 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi,yi+1,zi,coords); ptrs.sp[ii++] = w010; }
+	if (w110 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi+1,yi+1,zi,coords); ptrs.sp[ii++] = w110; }
+	if (w001 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi,yi,zi+1,coords); ptrs.sp[ii++] = w001; }
+	if (w101 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi+1,yi,zi+1,coords); ptrs.sp[ii++] = w101; }
+	if (w011 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi,yi+1,zi+1,coords); ptrs.sp[ii++] = w011; }
+	if (w111 > 1e-6) { ptrs.irp[ii] = this->ijk2indx(xi+1,yi+1,zi+1,coords); ptrs.sp[ii++] = w111; }
+      }
+    }
+  }
+  ptrs.jcp[ji] = ii;
+  ptrs.nz = ii;
+  return(ptrs);
+} EddyCatch
+
+void BiasFieldEstimatorImpl::multiply_At_star_CSC_by_image(const NEWIMAGE::volume<float>& ima,
+							   const NEWIMAGE::volume<float>& mask,
+							   unsigned int* const            jcp,
+							   double* const                  sp) const EddyTry
+{
+  unsigned int indx = 0;
+  for (int k=0; k<ima.zsize(); k++) {
+    for (int j=0; j<ima.ysize(); j++) {
+      for (int i=0; i<ima.xsize(); i++) {
+	if (mask(i,j,k)) { for (unsigned int ii=jcp[indx]; ii<jcp[indx+1]; ii++) sp[ii] *= ima(i,j,k); }
+	indx++;
+      }
+    }
+  }
+  return;
+} EddyCatch
+
+/*!
+ * Creates what is effectively a "resampling matrix" (A*) for a "reference"
+ * location. It is recommended that this is the average location of all volumes,
+ * but it can in principle also be a "zero| position in which case the
+ * resampling is simply an identity matrix.
+ * \throw EDDY::EddyException if there is a mismatch between mask an previous reference
+ */
+void BiasFieldEstimatorImpl::SetRefScan(const NEWIMAGE::volume<float>& mask,    //<! [in] mask
+					const EDDY::ImageCoordinates&  coords)  //<! [in] coordinates of desired reference
+EddyTry
+{
+  if (!_rAt.Nrows()) { // If no ref scan has been set before
+    _isz.resize(3);
+    _isz[0] = static_cast<unsigned int>(mask.xsize());
+    _isz[1] = static_cast<unsigned int>(mask.ysize());
+    _isz[2] = static_cast<unsigned int>(mask.zsize());
+    _vxs.resize(3);
+    _vxs[0] = static_cast<double>(mask.xdim());
+    _vxs[1] = static_cast<double>(mask.ydim());
+    _vxs[2] = static_cast<double>(mask.zdim());
+  }
+  else { // Check that new ref scan is compatible with old
+    if (static_cast<unsigned int>(mask.xsize()) != _isz[0] ||
+	static_cast<unsigned int>(mask.ysize()) != _isz[1] ||
+	static_cast<unsigned int>(mask.zsize()) != _isz[2]) {
+      throw EddyException("BiasFieldEstimatorImpl::SetRefScan: Size mismatch between new and previously set ref image");
+    }
+  }
+  PointersForSpMat ptrs = make_At_star_CSC(coords,mask);
+  _rAt = MISCMATHS::SpMat<float>(coords.N(),coords.N(),ptrs.irp.get(),ptrs.jcp.get(),ptrs.sp.get());
+  return;
+} EddyCatch
+
+/*!
+ * Used to add a scan to _Aty and _AtA.
+ * \throw EDDY::EddyException if there is a mismatch between input images
+ */
+void BiasFieldEstimatorImpl::AddScan(const NEWIMAGE::volume<float>& predicted, //<! [in] Predicted image
+				     const NEWIMAGE::volume<float>& observed,  //<! [in] Observed image
+				     const NEWIMAGE::volume<float>& mask,      //<! [in] Mask
+				     const EDDY::ImageCoordinates&  coords)    //<! [in] Coordinates where observed image was "observed"
+EddyTry
+{
+  static int cnt = 0;
+
+  if (predicted.xsize() != observed.xsize() ||
+      predicted.ysize() != observed.ysize() ||
+      predicted.zsize() != observed.zsize()) {
+    throw EddyException("BiasFieldEstimatorImpl::AddScan: Size mismatch between predicted and observed image");
+  }
+  if (!_rAt.Nrows()) { // If no ref scan has been set yet
+    throw EddyException("BiasFieldEstimatorImpl::AddScan: Attempting to add scan before ref scan has been set");
+  }
+  if (static_cast<unsigned int>(predicted.xsize()) != _isz[0] ||
+      static_cast<unsigned int>(predicted.ysize()) != _isz[1] ||
+      static_cast<unsigned int>(predicted.zsize()) != _isz[2]) {
+    throw EddyException("BiasFieldEstimatorImpl::AddScan: Size mismatch between predicted and previously set ref image");
+  }
+
+  NEWMAT::ColumnVector v_predicted = predicted.vec();
+  NEWMAT::ColumnVector v_observed = observed.vec();
+  NEWMAT::ColumnVector v_mask = mask.vec();
+  _nima++;
+  cnt++;
+  PointersForSpMat ptrs = make_At_star_CSC(coords,mask);
+  MISCMATHS::SpMat<float> At(coords.N(),coords.N(),ptrs.irp.get(),ptrs.jcp.get(),ptrs.sp.get());
+  At -= _rAt;
+  At.MultiplyColumns(NEWMAT::SP(v_predicted,v_mask));
+  if (!_Aty.Nrows()) _Aty = At * NEWMAT::SP((v_observed - v_predicted),v_mask);
+  else _Aty += At * NEWMAT::SP((v_observed - v_predicted),v_mask);
+  if (!_AtA.Nrows()) _AtA = At*At.t();
+  else _AtA += At*At.t();
+} EddyCatch
+
+/*!
+ * Fits a "direct" representation of the bias-field to the data and return it in newimage format
+ * \return A "direct" representation of the field in newimage format
+ * \throw EDDY::EddyException if object not ready to estimate field
+ */
+NEWIMAGE::volume<float> BiasFieldEstimatorImpl::GetField(double lambda) const  //<! [in] Weight of Laplacian regularisation when fitting the field
+EddyTry
+{
+  if (_nima==0) {
+    throw EddyException("BiasFieldEstimatorImpl::GetField(double): The field cannot be estimated until images have been loaded");
+  }
+  // Make (approximate) Bending energy regularisation matrix
+  MISCMATHS::SpMat<float> StS = MISCMATHS::Sparse3DBendingEnergyHessian(_isz,_vxs,MISCMATHS::PERIODIC);
+  // Solve for field
+  NEWMAT::ColumnVector field = (_AtA + (lambda/static_cast<double>(nvox())) * StS).SolveForx(_Aty,MISCMATHS::SYM_POSDEF,1.0e-6,1000);
+  // Retreive field as NEWIMAGE volume
+  NEWIMAGE::volume<float> fima(static_cast<int>(_isz[0]),static_cast<int>(_isz[1]),static_cast<int>(_isz[2]));
+  fima.setxdim(_vxs[0]); fima.setydim(_vxs[1]); fima.setzdim(_vxs[2]);
+  fima.insert_vec(field);
+  return(fima);
+} EddyCatch
+
+/*!
+ *  Fits a spline basis-field of the bias-field to the data and return it in newimage format
+ *  \return A basis-field representation of the field in newimage format
+ * \throw EDDY::EddyException if object not ready to estimate field
+ */
+NEWIMAGE::volume<float> BiasFieldEstimatorImpl::GetField(double ksp,           //<! [in] Knot-spacing of spline field in mm (approximate)
+							 double lambda) const  //<! [in] Weight of Bending Energy regularisation when fitting the field
+EddyTry
+{
+  if (_nima==0) {
+    throw EddyException("BiasFieldEstimatorImpl::GetField(double, double): The field cannot be estimated until images have been loaded");
+  }
+  // knot-spacing mm->voxels
+  std::vector<unsigned int> iksp(3);
+  iksp[0] = static_cast<unsigned int>((ksp / _vxs[0]) + 0.5);
+  iksp[1] = static_cast<unsigned int>((ksp / _vxs[1]) + 0.5);
+  iksp[2] = static_cast<unsigned int>((ksp / _vxs[2]) + 0.5);
+  // Make splinefield
+  BASISFIELD::splinefield spf(_isz,_vxs,iksp);
+
+  _AtA.Save("sum_AtA.txt");
+  MISCMATHS::write_ascii_matrix(_Aty,"sum_Aty.txt");
+
+  // Calculate B^T(\Sum A^T A)B and B^T(\Sum A^T \hat{f}) from paper
+  std::shared_ptr<MISCMATHS::SpMat<float> > B = spf.J();
+  MISCMATHS::SpMat<float> BtAtAB = B->t() * _AtA * *B;
+  BtAtAB += (lambda/static_cast<double>(nvox())) * (*spf.BendEnergyHessAsSpMat());
+  NEWMAT::ColumnVector BtAtf = B->t() * _Aty;
+  // Solve for spline coefficients of field
+  NEWMAT::ColumnVector b = BtAtAB.SolveForx(BtAtf,MISCMATHS::SYM_POSDEF,1.0e-6,200);
+  // Set coefficients in field
+  spf.SetCoef(b);
+  // Retrieve field as NEWIMAGE volume
+  NEWIMAGE::volume<float> field(static_cast<int>(_isz[0]),static_cast<int>(_isz[1]),static_cast<int>(_isz[2]));
+  field.setxdim(_vxs[0]); field.setydim(_vxs[1]); field.setzdim(_vxs[2]);
+  spf.AsVolume(field);
+
+  return(field);
+} EddyCatch
+
+MISCMATHS::SpMat<float> BiasFieldEstimatorImpl::GetAtMatrix(const EDDY::ImageCoordinates&  coords,
+							    const NEWIMAGE::volume<float>& predicted,
+							    const NEWIMAGE::volume<float>& mask) const EddyTry
+{
+  PointersForSpMat ptrs = make_At_star_CSC(coords,mask);
+  multiply_At_star_CSC_by_image(predicted,mask,ptrs.jcp.get(),ptrs.sp.get());
+  MISCMATHS::SpMat<float> At(coords.N(),coords.N(),ptrs.irp.get(),ptrs.jcp.get(),ptrs.sp.get());
+  return(At);
+} EddyCatch
+
+MISCMATHS::SpMat<float> BiasFieldEstimatorImpl::GetAtStarMatrix(const EDDY::ImageCoordinates&  coords,
+								const NEWIMAGE::volume<float>& mask) const EddyTry
+{
+  PointersForSpMat ptrs = make_At_star_CSC(coords,mask);
+  MISCMATHS::SpMat<float> At(coords.N(),coords.N(),ptrs.irp.get(),ptrs.jcp.get(),ptrs.sp.get());
+  return(At);
+} EddyCatch
+
+NEWIMAGE::volume<float> BiasFieldEstimatorImpl::GetATimesField(const EDDY::ImageCoordinates&  coords,
+							       const NEWIMAGE::volume<float>& predicted,
+							       const NEWIMAGE::volume<float>& mask,
+							       const NEWIMAGE::volume<float>& field) const EddyTry
+{
+  if (!_rAt.Nrows()) { // If no ref scan has been set yet
+    throw EddyException("BiasFieldEstimatorImpl::GetATimesField: No ref scan set yet");
+  }
+  PointersForSpMat ptrs = make_At_star_CSC(coords,mask);
+  MISCMATHS::SpMat<float> At(coords.N(),coords.N(),ptrs.irp.get(),ptrs.jcp.get(),ptrs.sp.get());
+  At -= _rAt;
+  At.MultiplyColumns(NEWMAT::SP(predicted.vec(),mask.vec()));
+  MISCMATHS::SpMat<float> A = At.t();
+  NEWMAT::ColumnVector fvec = field.vec();
+  NEWMAT::ColumnVector Ab = A*fvec;
+  NEWIMAGE::volume<float> Ab_ima = predicted;
+  Ab_ima.insert_vec(Ab);
+  return(Ab_ima);
+} EddyCatch
+
+BiasFieldEstimator::BiasFieldEstimator() EddyTry
+{
+  _pimpl = new BiasFieldEstimatorImpl();
+} EddyCatch
+
+BiasFieldEstimator::~BiasFieldEstimator()  { delete _pimpl; }
+
+void BiasFieldEstimator::SetRefScan(const NEWIMAGE::volume<float>& mask, const EDDY::ImageCoordinates& coords) EddyTry
+{
+  _pimpl->SetRefScan(mask,coords);
+} EddyCatch
+
+void BiasFieldEstimator::AddScan(const NEWIMAGE::volume<float>& predicted, const NEWIMAGE::volume<float>& observed,
+				 const NEWIMAGE::volume<float>& mask, const EDDY::ImageCoordinates& coords) EddyTry
+{
+  _pimpl->AddScan(predicted,observed,mask,coords);
+} EddyCatch
+
+NEWIMAGE::volume<float> BiasFieldEstimator::GetField(double lambda) const EddyTry { return(_pimpl->GetField(lambda)); } EddyCatch
+
+NEWIMAGE::volume<float> BiasFieldEstimator::GetField(double ksp, double lambda) const EddyTry { return(_pimpl->GetField(ksp,lambda)); } EddyCatch
+
+MISCMATHS::SpMat<float> BiasFieldEstimator::GetAtMatrix(const EDDY::ImageCoordinates&  coords,
+							const NEWIMAGE::volume<float>& predicted,
+							const NEWIMAGE::volume<float>& mask) const EddyTry
+{
+  return(_pimpl->GetAtMatrix(coords,predicted,mask));
+} EddyCatch
+
+MISCMATHS::SpMat<float> BiasFieldEstimator::GetAtStarMatrix(const EDDY::ImageCoordinates&  coords,
+							    const NEWIMAGE::volume<float>& mask) const EddyTry
+{
+  return(_pimpl->GetAtStarMatrix(coords,mask));
+} EddyCatch
+
+NEWIMAGE::volume<float> BiasFieldEstimator::GetATimesField(const EDDY::ImageCoordinates&  coords,
+							   const NEWIMAGE::volume<float>& predicted,
+							   const NEWIMAGE::volume<float>& mask,
+							   const NEWIMAGE::volume<float>& field) EddyTry
+{
+  return(_pimpl->GetATimesField(coords,predicted,mask,field));
+} EddyCatch
+
+void BiasFieldEstimator::Write(const std::string& basename) const EddyTry { _pimpl->Write(basename); } EddyCatch
+
+} // End namespace EDDY
+
+#endif // End #ifndef BiasFieldEstimatorImpl_h

CPUStackResampler.cpp

+/*! \file CPUStackResampler.cpp
+    \brief Contains definitions of a class for spline/tri-linear resampling of irregularly sampled columns in the z-direction.
+
+    \author Jesper Andersson
+    \version 1.0b, May, 2021.
+*/
+//
+// CPUStackResampler.cpp
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2021 University of Oxford 
+//
+
+#include <algorithm>
+#include "CPUStackResampler.h"
+
+namespace EDDY {
+
+CPUStackResampler::CPUStackResampler(const NEWIMAGE::volume<float>&  stack,
+				     const NEWIMAGE::volume<float>&  zcoord,
+				     const NEWIMAGE::volume<float>&  pred,
+				     const NEWIMAGE::volume<float>&  mask,
+				     double                          lambda)
+{
+  _ovol = stack;
+  _ovol = 0.0;
+  _omask = _ovol;
+  spline_interpolate_slice_stack(stack,zcoord,mask,lambda,&pred,_ovol,_omask);
+}
+
+CPUStackResampler::CPUStackResampler(const NEWIMAGE::volume<float>&  stack,
+				     const NEWIMAGE::volume<float>&  zcoord,
+				     const NEWIMAGE::volume<float>&  mask,
+				     NEWIMAGE::interpolation         interp,
+				     double                          lambda)
+{
+  _ovol = stack;
+  _ovol = 0.0;
+  _omask = _ovol;
+  if (interp==NEWIMAGE::spline) {
+    spline_interpolate_slice_stack(stack,zcoord,mask,lambda,nullptr,_ovol,_omask);
+  }
+  else if (interp==NEWIMAGE::trilinear) {
+    linear_interpolate_slice_stack(stack,zcoord,mask,_ovol,_omask);
+  }
+}
+
+void CPUStackResampler::spline_interpolate_slice_stack(// Input
+						       const NEWIMAGE::volume<float>&   slice_stack,
+						       const NEWIMAGE::volume<float>&   z_coord,
+						       const NEWIMAGE::volume<float>&   stack_mask,
+						       double                           lambda,
+						       // Optional input
+						       const NEWIMAGE::volume<float>    *pred_ptr,
+						       // Output
+						       NEWIMAGE::volume<float>&         ovol,
+						       NEWIMAGE::volume<float>&         omask) EddyTry
+{
+  // Get regularisation and regular sampline spline matrices once and for all
+  arma::Mat<float> StS = get_StS(ovol.zsize(),lambda);
+  arma::Mat<float> W = get_regular_W(ovol.zsize());
+  for (int i=0; i<ovol.xsize(); i++) {
+    for (int j=0; j<ovol.ysize(); j++) {
+      arma::Mat<float> Wir = get_Wir(z_coord,i,j);
+      arma::Col<float> y = get_y(slice_stack,i,j);
+      arma::Col<float> interpolated_column;
+      if (pred_ptr == nullptr) { // If we don't use predictions for support
+	interpolated_column = W * solve(Wir.t()*Wir + StS,Wir.t()*y);
+      }
+      else { // If we want to use predictions
+	std::vector<float> sorted_zcoords = sort_zcoord(z_coord,i,j);
+	arma::Mat<float> PW = get_prediction_weights(sorted_zcoords);
+	arma::Mat<float> WirW = arma::join_cols(Wir,PW*W);
+	arma::Col<float> y_pred = arma::join_cols(y,PW*get_y(*pred_ptr,i,j));
+	interpolated_column = W * solve(WirW.t()*WirW + StS,WirW.t()*y_pred);
+      }
+      for (int k=0; k<ovol.zsize(); k++) ovol(i,j,k) = interpolated_column[k]; // Insert column
+    }
+  }
+  omask = stack_mask; // Revisit
+} EddyCatch
+
+void CPUStackResampler::linear_interpolate_slice_stack(// Input
+						       const NEWIMAGE::volume<float>&   slice_stack,
+						       const NEWIMAGE::volume<float>&   z_coord,
+						       const NEWIMAGE::volume<float>&   stack_mask,
+						       // Output
+						       NEWIMAGE::volume<float>&         ovol,
+						       NEWIMAGE::volume<float>&         omask) EddyTry
+{
+  struct triplet {
+    triplet(float zz, float ii, float mm) : z(zz), i(ii), m(mm) {}
+    triplet() : z(0.0), i(0.0), m(0.0) {}
+    float z, i, m; // z-coord, ima-value, vaild_mask
+  };
+  // Allocate vector for sorting
+  std::vector<triplet> z_col(ovol.zsize()); 
+  // Do the interpolation
+  for (int j=0; j<ovol.ysize(); j++) {
+    for (int i=0; i<ovol.xsize(); i++) {
+      // Repack z-column into vector and sort if needed
+      z_col[0] = triplet(z_coord(i,j,0),slice_stack(i,j,0),stack_mask(i,j,0));
+      bool needs_sorting = false;
+      for (int k=1; k<ovol.zsize(); k++) {
+	z_col[k] = triplet(z_coord(i,j,k),slice_stack(i,j,k),stack_mask(i,j,k));
+	if (z_col[k].z < z_col[k-1].z) needs_sorting = true;
+      }
+      if (needs_sorting) std::sort(z_col.begin(),z_col.end(),[](const triplet& a, const triplet& b) { return(a.z < b.z); });
+      // Here starts the actual interpolation
+      for (int k=0; k<ovol.zsize(); k++) {
+	int kk=0;
+	for (kk=0; kk<ovol.zsize(); kk++) if (z_col[kk].z > k) break;
+	if (kk==0) { 
+	  if (z_col[kk].z < 0.5 && z_col[kk].m) { ovol(i,j,k) = z_col[kk].i; omask(i,j,k) = 1; }
+	  else { ovol(i,j,k) = 0; omask(i,j,k) = 0; }
+	}
+	else if (kk==ovol.zsize()) { 
+	  if (z_col[kk-1].z > ovol.zsize() - 0.5 && z_col[kk-1].m) { ovol(i,j,k) = z_col[kk-1].i; omask(i,j,k) = 1; }
+	  else { ovol(i,j,k) = 0; omask(i,j,k) = 0; }
+	}
+	else {
+	  if (z_col[kk-1].m && z_col[kk].m) {
+	    ovol(i,j,k) = z_col[kk-1].i + (k-z_col[kk-1].z) * (z_col[kk].i-z_col[kk-1].i) / (z_col[kk].z-z_col[kk-1].z);
+	    omask(i,j,k) = 1;
+	  }
+	  else { ovol(i,j,k) = 0; omask(i,j,k) = 0; }
+	}
+      }
+    }
+  }  
+  return;
+} EddyCatch
+
+arma::Mat<float> CPUStackResampler::get_StS(int sz, float lambda) const EddyTry
+{
+  arma::Mat<float> StS(sz,sz,arma::fill::zeros);
+  StS(0,0) = 6.0*lambda; StS(0,1) = -4.0*lambda; StS(0,2) = lambda; StS(0,sz-2) = lambda; StS(0,sz-1) = -4.0*lambda; 
+  StS(1,0) = -4.0*lambda; StS(1,1) = 6.0*lambda; StS(1,2) = -4.0*lambda; StS(1,3) = lambda; StS(1,sz-1) = lambda;
+  for (int i=2; i<(sz-2); i++) {
+    StS(i,i-2) = lambda; StS(i,i-1) = -4.0*lambda; StS(i,i) = 6.0*lambda; StS(i,i+1) = -4.0*lambda; StS(i,i+2) = lambda;
+  }
+  StS(sz-2,sz-4) = lambda; StS(sz-2,sz-3) = -4.0*lambda; StS(sz-2,sz-2) = 6.0*lambda; StS(sz-2,sz-1) = -4.0*lambda; StS(sz-2,0) = lambda;
+  StS(sz-1,sz-3) = lambda; StS(sz-1,sz-2) = -4.0*lambda; StS(sz-1,sz-1) = 6.0*lambda; StS(sz-1,0) = -4.0*lambda; StS(sz-1,1) = lambda;
+  return(StS);
+} EddyCatch
+
+arma::Mat<float> CPUStackResampler::get_regular_W(int sz) const EddyTry
+{
+  arma::Mat<float> W(sz,sz,arma::fill::zeros);
+  W(0,0) = 5.0/6.0; W(0,1) = 1.0/6.0;
+  for (int i=1; i<(sz-1); i++) {
+    W(i,i-1) = 1.0/6.0; W(i,i) = 4.0/6.0; W(i,i+1) = 1.0/6.0;
+  }
+  W(sz-1,sz-2) = 1.0/6.0; W(sz-1,sz-1) = 5.0/6.0;
+  return(W);
+} EddyCatch
+
+arma::Mat<float> CPUStackResampler::get_Wir(const NEWIMAGE::volume<float>& zcoord,
+					    int i, int j) const EddyTry
+{
+  arma::Mat<float> Wir(zcoord.zsize(),zcoord.zsize(),arma::fill::zeros);
+  for (int k=0; k<zcoord.zsize(); k++) {
+    if (zcoord(i,j,k)>=0 && zcoord(i,j,k)<=(zcoord.zsize()-1)) { // If in valid range
+      int iz = static_cast<int>(zcoord(i,j,k));
+      for (int c=iz-2; c<iz+3; c++) {
+	Wir(k,std::min(std::max(0,c),static_cast<int>(zcoord.zsize()-1))) += wgt_at(zcoord(i,j,k)-static_cast<float>(c));
+      }
+    }
+  }
+  return(Wir);
+} EddyCatch
+
+float CPUStackResampler::wgt_at(float x) const EddyTry
+{
+  float wgt = 0.0;
+  x = (x<0.0) ? -x : x;
+  if (x < 1) wgt = 2.0/3.0 + 0.5*x*x*(x-2.0);
+  else if (x < 2) wgt = (1.0/6.0) * (2.0-x)*(2.0-x)*(2.0-x);
+
+  return(wgt);
+} EddyCatch
+
+std::vector<float> CPUStackResampler::sort_zcoord(const NEWIMAGE::volume<float>& zcoord,
+						  int i, int j) const EddyTry
+{
+  std::vector<float> ovec(zcoord.zsize());
+  bool needs_sorting = false;
+  ovec[0] = zcoord(i,j,0);
+  for (int k=1; k<zcoord.zsize(); k++) {
+    ovec[k] = zcoord(i,j,k);
+    if (ovec[k] < ovec[k-1]) needs_sorting = true;
+  }
+  if (needs_sorting) std::sort(ovec.begin(),ovec.end());
+
+  return(ovec);
+} EddyCatch
+
+arma::Mat<float> CPUStackResampler::get_prediction_weights(const std::vector<float> zcoord) const EddyTry
+{
+  arma::Mat<float> wgts(zcoord.size(),zcoord.size(),arma::fill::zeros);
+  for (unsigned int k=0; k<zcoord.size(); k++) {
+    unsigned int i=0;
+    for (i=0; i<zcoord.size(); i++) if (zcoord[i] > k) break;
+    if (i==0) {
+      wgts(k,k) = std::min(zcoord[i]-k,1.0f);
+    }
+    else if (i==zcoord.size()) {
+      wgts(k,k) = std::min(k-zcoord[i-1],1.0f);
+    }
+    else {
+      float gap = zcoord[i]-zcoord[i-1];
+      if (gap < 1.0) { // If gap < one voxel
+	wgts(k,k) = 0.0;
+      }
+      else if (gap < 2.0 && std::max(k-zcoord[i-1],zcoord[i]-k) < 1.0) { // If gap < 2 voxels and only one prediction in gap
+	wgts(k,k) = gap - 1.0f;
+      }
+      else { // If there is more than one prediction in gap
+	wgts(k,k) = std::min(1.0f,std::min(k-zcoord[i-1],zcoord[i]-k));
+      }
+    }
+    if (wgts(k,k) > 1e-12) wgts(k,k) = std::sqrt(wgts(k,k)); // Avoid taking sqrt of very small value
+  }
+  return(wgts);
+} EddyCatch
+
+} // End namespace EDDY

CPUStackResampler.h

+/*! \file CPUStackResampler.h
+    \brief Contains declaration of a class for spline/tri-linear resampling of irregularly sampled columns in the z-direction.
+
+    \author Jesper Andersson
+    \version 1.0b, May, 2021.
+*/
+//
+// CPUStackResampler.h
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2021 University of Oxford
+//
+
+#ifndef CPUStackResampler_h
+#define CPUStackResampler_h
+
+#include <cstdlib>
+#include <vector>
+#include <armadillo>
+#include "newimage/newimageall.h"
+#include "EddyHelperClasses.h"
+
+namespace EDDY {
+
+class CPUStackResampler
+{
+public:
+  /// Constructor. Performs the actual work so that when the object is created there is already a resampled image ready.
+
+  /// This version of the constructor uses a predicted volume and the laplacian for regularisation
+  CPUStackResampler(const NEWIMAGE::volume<float>&  stack,
+		    const NEWIMAGE::volume<float>&  zcoord,
+		    const NEWIMAGE::volume<float>&  pred,
+		    const NEWIMAGE::volume<float>&  mask,
+		    double                          lambda=0.005);
+
+  /// This version of the constructor uses either splines and Laplacian regularisation or linear interpolation.
+  CPUStackResampler(const NEWIMAGE::volume<float>&  stack,
+		    const NEWIMAGE::volume<float>&  zcoord,
+		    const NEWIMAGE::volume<float>&  mask,
+		    NEWIMAGE::interpolation         interp=NEWIMAGE::spline,
+		    double                          lambda=0.005);
+
+  ~CPUStackResampler() {}
+
+  // Returns interpolated image
+  const NEWIMAGE::volume<float>& GetImage() const EddyTry { return(_ovol); } EddyCatch
+  // Returns mask
+  const NEWIMAGE::volume<float>& GetMask() const EddyTry { return(_omask); } EddyCatch
+
+private:
+  NEWIMAGE::volume<float> _ovol;
+  NEWIMAGE::volume<float> _omask;
+  void spline_interpolate_slice_stack(// Input
+				      const NEWIMAGE::volume<float>&   slice_stack,
+				      const NEWIMAGE::volume<float>&   z_coord,
+				      const NEWIMAGE::volume<float>&   stack_mask,
+				      double                           lambda,
+				      // Optional input
+				      const NEWIMAGE::volume<float>    *pred_ptr,
+				      // Output
+				      NEWIMAGE::volume<float>&         ovol,
+				      NEWIMAGE::volume<float>&         omask);
+  void linear_interpolate_slice_stack(// Input
+				      const NEWIMAGE::volume<float>&   slice_stack,
+				      const NEWIMAGE::volume<float>&   z_coord,
+				      const NEWIMAGE::volume<float>&   stack_mask,
+				      // Output
+				      NEWIMAGE::volume<float>&         ovol,
+				      NEWIMAGE::volume<float>&         omask);
+  arma::Mat<float> get_StS(int sz, float lambda) const;
+  arma::Mat<float> get_regular_W(int sz) const;
+  arma::Mat<float> get_Wir(const NEWIMAGE::volume<float>& zcoord,
+					    int i, int j) const;
+  arma::Col<float> get_y(const NEWIMAGE::volume<float>& stack,
+			 int i, int j) const EddyTry {
+    arma::Col<float> y(stack.zsize());
+    for (int k=0; k<stack.zsize(); k++) y[k] = stack(i,j,k);
+    return(y);
+  } EddyCatch
+  float wgt_at(float x) const;
+  std::vector<float> sort_zcoord(const NEWIMAGE::volume<float>& zcoord,
+				 int i, int j) const;
+  arma::Mat<float> get_prediction_weights(const std::vector<float> zcoord) const;
+};
+
+} // End namespace EDDY
+
+#endif // End #ifndef CPUStackResampler_h

DWIPredictionMaker.h

+/*! \file DWIPredictionMaker.h
+    \brief Contains declaration of virtual base class for making predictions about DWI data.
+
+    \author Jesper Andersson
+    \version 1.0b, Sep., 2012.
+*/
+// Declarations of virtual base class for
+// making predictions about DWI data.
+//
+// DWIPredictionMaker.h
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2011 University of Oxford
+//
+
+#ifndef DWIPredictionMaker_h
+#define DWIPredictionMaker_h
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cmath>
+#include "armawrap/newmat.h"
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+
+namespace EDDY {
+
+
+/****************************************************************//**
+*
+* \brief Virtual base class for classes used to make predictions
+* about diffusion data.
+*
+* The idea of the prediction makers is to be able to provide them
+* with some set of data (for some set of diffusion gradients) and
+* then use it to make predictions about what a certain data point
+* should be. The predictions could be about unobserved points (in
+* which case it would perform an interpolation) or about observed
+* points (which would amount to smoothing).
+* The virtual class provides a minimal interface for actual prediction
+* makers. These could be based e.g. on the diffusion tensor model
+* or on Gaussian processes.
+********************************************************************/
+class DWIPredictionMaker
+{
+public:
+  DWIPredictionMaker() {}
+  virtual ~DWIPredictionMaker() {}
+  /// Returns prediction for point given by indx.
+  virtual NEWIMAGE::volume<float> Predict(unsigned int indx, bool exclude=false) const = 0;
+  /// Returns prediction for point given by indx.
+  virtual NEWIMAGE::volume<float> Predict(unsigned int indx, bool exclude=false) = 0;
+  /// Returns prediction for point given by indx. This is used only as a means of directly comparing CPU and GPU outputs.
+  virtual NEWIMAGE::volume<float> PredictCPU(unsigned int indx, bool exclude=false) = 0;
+  /// Returns predictions for a set of points given by indicies
+  virtual std::vector<NEWIMAGE::volume<float> > Predict(const std::vector<unsigned int>& indicies, bool exclude=false) = 0;
+  /// Returns input data for point given by indx
+  virtual NEWIMAGE::volume<float> InputData(unsigned int indx) const = 0;
+  /// Returns input data for points given by indicies
+  virtual std::vector<NEWIMAGE::volume<float> > InputData(const std::vector<unsigned int>& indicies) const = 0;
+  /// Returns variance of prediction for point given by indx.
+  virtual double PredictionVariance(unsigned int indx, bool exclude=false) = 0;
+  /// Returns measurement-error variance for point given by indx.
+  virtual double ErrorVariance(unsigned int indx) const = 0;
+  /// Returns true if all data has been loaded
+  virtual bool IsPopulated() const = 0;
+  /// Indicates if it is ready to make predictions.
+  virtual bool IsValid() const = 0;
+  /// Specify the # of points we plan to put into the predictor.
+  virtual void SetNoOfScans(unsigned int n) = 0;
+  /*
+  /// Adds a new point to the end of the current list. This function is NOT thread safe.
+  virtual void AddScan(const NEWIMAGE::volume<float>& scan,  // NOT thread safe
+                       const DiffPara&                dp) = 0;
+  */
+  /// Set a point given by indx. This function is thread safe as long as different threads set different points.
+  virtual void SetScan(const NEWIMAGE::volume<float>& scan,  // May be thread safe if used "sensibly"
+                       const DiffPara&                dp,
+                       unsigned int                   indx,
+		       unsigned int                   sess) = 0;
+  /// Set a point given by indx. This function is thread safe as long as different threads set different points.
+  virtual void SetScan(const NEWIMAGE::volume<float>& scan,  // May be thread safe if used "sensibly"
+                       const DiffPara&                dp,
+                       unsigned int                   indx) = 0;
+  /// Returns the number of hyperparameters for the model
+  virtual unsigned int NoOfHyperPar() const = 0;
+  /// Returns the hyperparameters for the model
+  virtual std::vector<double> GetHyperPar() const = 0;
+  /// Evaluates the model so as to make the predictor ready to make predictions.
+  virtual void EvaluateModel(const NEWIMAGE::volume<float>& mask, bool verbose=false) = 0;
+  /// Evaluates the model so as to make the predictor ready to make predictions.
+  virtual void EvaluateModel(const NEWIMAGE::volume<float>& mask, float fwhm, bool verbose=false) = 0;
+  /// Writes internal content to disk for debug purposes
+  virtual void WriteImageData(const std::string& fname) const = 0;
+  virtual void WriteMetaData(const std::string& fname) const = 0;
+  virtual void Write(const std::string& fname) const = 0;
+};
+
+} // End namespace EDDY
+
+#endif // End #ifndef DWIPredictionMaker_h
+
+/*!
+  \fn NEWIMAGE::volume<float> DWIPredictionMaker::Predict(unsigned int indx) const = 0
+  Returns a prediction for the scan given by indx, where indx pertains to the indx that was used when setting a given
+  scan with a call to SetScan. So if for example there was a call PM.SetScan(scan,my_dp,5) then PM.Predict(5) will return
+  a prediction for the diffusion weighting specified in my_dp.
+  \param indx Specifies which scan, and indirectly what diffusion weighting, we want the prediction for. It is zero-offset so should be in the range 0--n-1 where n has been set by an earlier call to DWIPredictionMaker::SetNoOfScans(unsigned int n).
+  \return A predicted image volume.
+*/
+
+/*!
+  \fn NEWIMAGE::volume<float> DWIPredictionMaker::Predict(const DiffPara& dpar) const = 0
+  Returns a prediction for a scan with diffusion parameters given by dpar. This may pertain to an observed scan (smoothing) or an unobserved scan (interpolation).
+  \param dpar Specifies the diffusion weighting.
+  \return A predicted image volume.
+*/
+
+/*!
+  \fn NEWIMAGE::volume4D<float> DWIPredictionMaker::PredictAll()
+  Returns predictions for all images/diffusion weightings that have been set in the prediction maker.
+  \return A 4D volume with as many volumes as there are images in the prediction maker (as set by an earlier call to DWIPredictionMaker::SetNoOfScans(unsigned int n)).
+*/
+
+/*!
+  \fn bool DWIPredictionMaker::IsPopulated() const
+  Returns true if valid scans have been set for all "slots" as
+  defined by a call to DiffusionGP::SetNoOfScans(unsigned int).
+*/
+
+/*!
+  \fn bool DWIPredictionMaker::IsValid() const
+  Will return true if the object is ready to make predictions. For this to be true it must have been fully populated.
+*/
+
+/*!
+  \fn void DWIPredictionMaker::SetNoOfScans(unsigned int n)
+  Specifies the number of scans to use for the prediction maker. It can override an earlier call to increase or shrink the number of scans.
+  \param n The number of scans
+*/
+
+/*!
+  \fn void DWIPredictionMaker::AddScan(const NEWIMAGE::volume<float>& scan, const DiffPara& dp)
+  Adds a scan to the end of the current list of scans in the prediction maker.
+  \param scan The scan that should be added.
+  \param dp Diffusion parameters pertaining to scan.
+*/
+
+/*!
+  \fn void DWIPredictionMaker::SetScan(const NEWIMAGE::volume<float>& scan, const DiffPara& dp, unsigned int indx)
+  Inserts a scan in the slot given by indx into the list of scans in the prediction maker.
+  \param scan The scan that should be inserted.
+  \param dp Diffusion parameters pertaining to scan.
+  \param indx Slot in scan list into which to put scan. It is zero-offset so should be in the range 0--n-1 where n has been set by an earlier call to DWIPredictionMaker::SetNoOfScans(unsigned int n).
+*/
+
+/*!
+  \fn void DWIPredictionMaker::EvaluateModel(const NEWIMAGE::volume<float>& mask)
+  Performs the calculations that need to be done after all the data has been loaded and before predictions can be made. The details of this will depend on the specific derived class. If for example the derived class is based on the diffusion tensor model it will entail calculating the diffusion tensor for each voxel.
+  \param mask Binary mask used to limit the calculations to where the voxels are set to one in the mask.
+*/

DiffusionGP.cpp

+/*! \file DiffusionGP.cpp
+    \brief Contains definitions for class for making Gaussian process based predictions about DWI data.
+
+    \author Jesper Andersson
+    \version 1.0b, Sep., 2013.
+*/
+// Definitions of class to make Gaussian-Process
+// based predictions about diffusion data.
+//
+// DiffusionGP.cpp
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2013 University of Oxford
+//
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cmath>
+#include "armawrap/newmat.h"
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "EddyHelperClasses.h"
+#include "EddyUtils.h"
+#include "KMatrix.h"
+#include "HyParEstimator.h"
+#include "DiffusionGP.h"
+
+using namespace EDDY;
+
+/****************************************************************//**
+*
+* Constructs a DiffusionGP object given files containing input data.
+* When using this constructor the object is immediately ready to do
+* predictions.
+* \param scans_fname Name of file containing multiple diffusion
+*        weighted volumes.
+* \param var_mask_fname
+* \param dpars Vector of objects of DiffPara type that specifies
+*        diffusion weighting and direction for the volumes in
+*        scans_fname.
+*
+*
+********************************************************************/
+DiffusionGP::DiffusionGP(const std::shared_ptr<const KMatrix>&        Kmat,
+			 const std::shared_ptr<const HyParEstimator>& hpe,
+			 const std::string&                           scans_fname,
+			 const std::string&                           brain_mask_fname,
+                         const std::vector<DiffPara>&                 dpars,
+			 float                                        fwhm,
+			 bool                                         verbose) EddyTry
+: _Kmat(Kmat->Clone()), _hpe(hpe->Clone()), _pop(false), _mc(false)
+{
+  NEWIMAGE::volume4D<float> scans;
+  EddyUtils::read_DWI_volume4D(scans,scans_fname,dpars);
+  for (int s=0; s<scans.tsize(); s++) _sptrs.push_back(std::shared_ptr<NEWIMAGE::volume<float> >(new NEWIMAGE::volume<float>(scans[s])));
+  _pop = true;
+  _dpars = EddyUtils::GetDWIDiffParas(dpars);
+  NEWIMAGE::volume<float> brain_mask; NEWIMAGE::read_volume(brain_mask,brain_mask_fname);
+  _Kmat->SetDiffusionPar(_dpars);
+  std::vector<std::vector<unsigned int> > mi = _Kmat->GetMeanIndicies();
+  mean_correct(mi);
+  DataSelector vd(_sptrs,brain_mask,_hpe->GetNVox(),fwhm,_hpe->RndInit());
+  _hpe->SetData(vd.GetData());
+  _hpe->Estimate(_Kmat,verbose);
+  _Kmat->SetHyperPar(hpe->GetHyperParameters());
+  _Kmat->CalculateInvK();
+} EddyCatch
+
+NEWIMAGE::volume<float> DiffusionGP::Predict(unsigned int indx,
+					     bool         exclude) EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::Predict:non-const: Not yet fully populated");
+  if (!IsValid()) throw EddyException("DiffusionGP::Predict:non-const: Not yet ready for predictions");
+  NEWMAT::RowVector pv = _Kmat->PredVec(indx,exclude); // Calls non-const version
+  NEWIMAGE::volume<float> pi = *_sptrs[0]; pi = 0.0;
+  #ifdef COMPILE_GPU
+  predict_image_gpu(indx,exclude,pv,pi);
+  #else
+  predict_image_cpu(indx,exclude,pv,pi);
+  #endif
+  return(pi);
+} EddyCatch
+
+NEWIMAGE::volume<float> DiffusionGP::PredictCPU(unsigned int indx,
+						bool         exclude) EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::Predict:non-const: Not yet fully populated");
+  if (!IsValid()) throw EddyException("DiffusionGP::Predict:non-const: Not yet ready for predictions");
+  NEWMAT::RowVector pv = _Kmat->PredVec(indx,exclude); // Calls non-const version
+  NEWIMAGE::volume<float> pi = *_sptrs[0]; pi = 0.0;
+  predict_image_cpu(indx,exclude,pv,pi);
+  return(pi);
+} EddyCatch
+
+NEWIMAGE::volume<float> DiffusionGP::Predict(unsigned int indx,
+					     bool         exclude) const EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::Predict:const: Not yet fully populated");
+  if (!IsValid()) throw EddyException("DiffusionGP::Predict:const: Not yet ready for predictions");
+  NEWMAT::RowVector pv = _Kmat->PredVec(indx,exclude); // Calls const version
+  NEWIMAGE::volume<float> pi = *_sptrs[0]; pi = 0.0;
+  #ifdef COMPILE_GPU
+  predict_image_gpu(indx,exclude,pv,pi);
+  #else
+  predict_image_cpu(indx,exclude,pv,pi);
+  #endif
+  return(pi);
+} EddyCatch
+
+std::vector<NEWIMAGE::volume<float> > DiffusionGP::Predict(const std::vector<unsigned int>& indicies,
+							   bool                             exclude) EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::Predict: Not yet fully populated");
+  if (!IsValid()) throw EddyException("DiffusionGP::Predict: Not yet ready for predictions");
+  std::vector<NEWIMAGE::volume<float> > pi(indicies.size());
+  std::vector<NEWMAT::RowVector> pvecs(indicies.size());
+  for (unsigned int i=0; i<indicies.size(); i++) pvecs[i] = _Kmat->PredVec(indicies[i],exclude);
+  #ifdef COMPILE_GPU
+  predict_images_gpu(indicies,exclude,pvecs,pi);
+  #else
+  predict_images_cpu(indicies,exclude,pvecs,pi);
+  #endif
+  return(pi);
+} EddyCatch
+
+NEWIMAGE::volume<float> DiffusionGP::InputData(unsigned int indx) const EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::InputData: Not yet fully populated");
+  if (indx >= _sptrs.size()) throw EddyException("DiffusionGP::InputData: indx out of range");
+  return(*(_sptrs[indx]) + *(_mptrs[which_mean(indx)]));
+} EddyCatch
+
+std::vector<NEWIMAGE::volume<float> > DiffusionGP::InputData(const std::vector<unsigned int>& indx) const EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::InputData: Not yet fully populated");
+  for (unsigned int i=0; i<indx.size(); i++) if (indx[i] >= _sptrs.size()) throw EddyException("DiffusionGP::InputData: indx out of range");
+  std::vector<NEWIMAGE::volume<float> > rval(indx.size());
+  for (unsigned int i=0; i<indx.size(); i++) rval[i] = *(_sptrs[indx[i]]) + *(_mptrs[which_mean(indx[i])]);
+  return(rval);
+} EddyCatch
+
+double DiffusionGP::PredictionVariance(unsigned int indx,
+				       bool         exclude) EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::PredictionVariance:const: Not yet fully populated");
+  if (!IsValid()) throw EddyException("DiffusionGP::PredictionVariance:const: Not yet ready for predictions");
+  double pv = _Kmat->PredVar(indx,exclude);
+  return(pv);
+} EddyCatch
+
+double DiffusionGP::ErrorVariance(unsigned int indx) const EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::ErrorVariance:const: Not yet fully populated");
+  if (!IsValid()) throw EddyException("DiffusionGP::ErrorVariance:const: Not yet ready for predictions");
+  double ev = _Kmat->ErrVar(indx);
+  return(ev);
+} EddyCatch
+
+void DiffusionGP::SetNoOfScans(unsigned int n) EddyTry
+{
+  if (n == _sptrs.size()) return; // No change
+  else if (n > _sptrs.size()) {   // If increasing size
+    _sptrs.resize(n,std::shared_ptr<NEWIMAGE::volume<float> >()); // New elements populated by NULL
+    _dpars.resize(n); // New elements populated with (arbitrary) [1 0 0] bvec
+    _Kmat->Reset();
+    _pop = false;
+  }
+  else { // Decreasing size not allowed
+    throw EddyException("DiffusionGP::SetNoOfScans: Decreasing size not allowed");
+  }
+  return;
+} EddyCatch
+
+void DiffusionGP::SetScan(const NEWIMAGE::volume<float>& scan,
+                          const DiffPara&                dp,
+                          unsigned int                   indx) EddyTry
+{
+  if (int(indx) > (int(_sptrs.size())-1)) throw EddyException("DiffusionGP::SetScan: Invalid image index");
+  if (_sptrs.size() && _sptrs[0] && !NEWIMAGE::samesize(*_sptrs[0],scan)) throw EddyException("DiffusionGP::SetScan: Wrong image dimension");
+  if (_sptrs[indx] && dp != _dpars[indx]) throw EddyException("DiffusionGP::SetScan: You cannot change shell or direction of scan");
+  _sptrs[indx] = std::shared_ptr<NEWIMAGE::volume<float> >(new NEWIMAGE::volume<float>(scan));
+  _dpars[indx] = dp;
+  _pop = is_populated();
+} EddyCatch
+
+void DiffusionGP::EvaluateModel(const NEWIMAGE::volume<float>& mask,
+				float                          fwhm,
+				bool                           verbose) EddyTry
+{
+  if (!IsPopulated()) throw EddyException("DiffusionGP::EvaluateModel: Predictor not fully populated");
+  _Kmat->Reset(); _Kmat->SetDiffusionPar(_dpars);
+  std::vector<std::vector<unsigned int> > mi = _Kmat->GetMeanIndicies();
+  mean_correct(mi);
+  DataSelector vd(_sptrs,mask,_hpe->GetNVox(),fwhm,_hpe->RndInit());
+  _hpe->SetData(vd.GetData());
+  _hpe->Estimate(_Kmat,verbose);
+  _Kmat->SetHyperPar(_hpe->GetHyperParameters());
+  _Kmat->CalculateInvK();
+} EddyCatch
+
+void DiffusionGP::WriteImageData(const std::string& fname) const EddyTry
+{
+  char ofname[256];
+  if (!IsPopulated()) throw EddyException("DiffusionGP::WriteImageData: Not yet fully populated");
+  // For practical reasons the volumes are written individually
+  for (unsigned int i=0; i<_sptrs.size(); i++) {
+    sprintf(ofname,"%s_%03d",fname.c_str(),i);
+    NEWIMAGE::write_volume(*(_sptrs[i]),ofname);
+  }
+  for (unsigned int i=0; i<_mptrs.size(); i++) {
+    sprintf(ofname,"%s_mean_%03d",fname.c_str(),i);
+    NEWIMAGE::write_volume(*(_mptrs[i]),ofname);
+  }
+} EddyCatch
+
+void DiffusionGP::mean_correct(const std::vector<std::vector<unsigned int> >& mi) EddyTry
+{
+  if (!_mc) {
+    // First, calculate mean images
+    if (_mptrs.size() != mi.size()) _mptrs.resize(mi.size());
+    for (unsigned int m=0; m<mi.size(); m++) {
+      if (_mptrs[m]==nullptr || !samesize(*_sptrs[0],*_mptrs[m])) {
+	_mptrs[m] = std::shared_ptr<NEWIMAGE::volume<float> >(new NEWIMAGE::volume<float>(*_sptrs[0]));
+      }
+      *_mptrs[m] = 0.0;
+      for (unsigned int s=0; s<mi[m].size(); s++) {
+	*_mptrs[m] += *_sptrs[mi[m][s]];
+      }
+      *_mptrs[m] /= mi[m].size();
+    }
+    // Next, correct individual images
+    for (unsigned int m=0; m<mi.size(); m++) {
+      for (unsigned int s=0; s<mi[m].size(); s++) {
+	*_sptrs[mi[m][s]] -= *_mptrs[m];
+      }
+    }
+    _mc = true;
+  }
+} EddyCatch
+
+bool DiffusionGP::is_populated() const EddyTry
+{
+  bool pop = true;
+  for (unsigned int i=0; i<_sptrs.size(); i++) {
+    if (!_sptrs[i]) { pop = false; break; }
+  }
+  return(pop);
+} EddyCatch
+
+void DiffusionGP::predict_image_cpu(// Input
+				    unsigned int             indx,
+				    bool                     exclude,
+				    const NEWMAT::RowVector& pv,
+				    // Output
+				    NEWIMAGE::volume<float>& pi) const EddyTry
+{
+  unsigned int mi = which_mean(indx);
+  unsigned int ys = (exclude) ? _sptrs.size()-1 : _sptrs.size();
+  NEWMAT::ColumnVector y(ys);
+  pi = *_mptrs[mi];
+  for (int k=0; k<pi.zsize(); k++) {
+    for (int j=0; j<pi.ysize(); j++) {
+      for (int i=0; i<pi.xsize(); i++) {
+	if (get_y(i,j,k,indx,exclude,y)) pi(i,j,k) += static_cast<float>((pv*y).AsScalar());
+	else pi(i,j,k) = 0.0;
+      }
+    }
+  }
+} EddyCatch
+
+void DiffusionGP::predict_images_cpu(// Input
+				     const std::vector<unsigned int>&       indicies,
+				     bool                                   exclude,
+				     const std::vector<NEWMAT::RowVector>&  pvecs,
+				     // Output
+				     std::vector<NEWIMAGE::volume<float> >& pi) const EddyTry
+{
+  if (indicies.size() != pvecs.size() || indicies.size() != pi.size()) {
+    throw EDDY::EddyException("DiffusionGP::predict_images_cpu: mismatch among indicies, pvecs and pi");
+  }
+  for (unsigned int i=0; i<indicies.size(); i++) predict_image_cpu(indicies[i],exclude,pvecs[i],pi[i]);
+  return;
+} EddyCatch
+
+unsigned int DiffusionGP::which_mean(unsigned int indx) const EddyTry
+{
+  std::vector<std::vector<unsigned int> > mis = _Kmat->GetMeanIndicies();
+  if (mis.size() == 1) return(0);
+  else {
+    for (unsigned int m=0; m<mis.size(); m++) {
+      unsigned int j;
+      for (j=0; j<mis[m].size(); j++) if (mis[m][j] == indx) break;
+      if (j<mis[m].size()) return(m);
+    }
+    throw EddyException("DiffusionGP::which_mean: Invalid indx");
+  }
+} EddyCatch
+
+bool DiffusionGP::get_y(// Input
+			unsigned int           i,
+			unsigned int           j,
+			unsigned int           k,
+			unsigned int           indx,
+			bool                   exclude,
+			// Output
+			NEWMAT::ColumnVector&  y) const EddyTry
+{
+  for (unsigned int t=0, tt=1; t<_sptrs.size(); t++) {
+    if (!exclude || t!=indx) {
+      if (!(*_sptrs[t])(i,j,k)) return(false);
+      else y(tt++) = (*_sptrs[t])(i,j,k);
+    }
+  }
+  return(true);
+} EddyCatch

DiffusionGP.h

+/*! \file DiffusionGP.h
+    \brief Contains declaration of class for making Gaussian process based predictions about DWI data.
+
+    \author Jesper Andersson
+    \version 1.0b, Sep., 2012.
+*/
+// Declarations of class to make Gaussian-Process
+// based predictions about diffusion data.
+//
+// DiffusionGP.h
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2011 University of Oxford
+//
+
+#ifndef DiffusionGP_h
+#define DiffusionGP_h
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <memory>
+#include "armawrap/newmat.h"
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "EddyHelperClasses.h"
+#include "DWIPredictionMaker.h"
+#include "KMatrix.h"
+#include "HyParEstimator.h"
+
+namespace EDDY {
+
+/****************************************************************//**
+*
+* \brief Class used to make Gaussian process based predictions
+* about diffusion data.
+*
+* Will make predictions about observed (smoothing) and unobserved
+* (interpolation) scans using Gaussian Processes. It works by first
+* creating an object, setting the # of scans one wants it to contain
+* and then populate it with scans. Once that is done it is ready
+* to start making predictions. If you think of the signal (in a
+* given voxel) as points on a surface the Gaussian process will make
+* predictions from the assumption that that surface should be smooth.
+*
+********************************************************************/
+class DiffusionGP : public DWIPredictionMaker
+{
+public:
+  /// Default constructor
+  DiffusionGP(const std::shared_ptr<const KMatrix>&        Kmat,
+	      const std::shared_ptr<const HyParEstimator>& hpe) EddyTry : _Kmat(Kmat->Clone()), _hpe(hpe->Clone()), _pop(true), _mc(false) {} EddyCatch
+  /// Constructor that takes filenames from which to load data
+  DiffusionGP(const std::shared_ptr<const KMatrix>&        Kmat,
+	      const std::shared_ptr<const HyParEstimator>& hpe,
+	      const std::string&                           scans_fname,
+              const std::string&                           var_mask_fname,
+              const std::vector<DiffPara>&                 dpars,
+	      float                                        fwhm=0.0,
+	      bool                                         verbose=false);
+  /// Returns prediction for point given by indx.
+  virtual NEWIMAGE::volume<float> Predict(unsigned int indx, bool exclude=false) const;
+  /// Returns prediction for point given by indx.
+  virtual NEWIMAGE::volume<float> Predict(unsigned int indx, bool exclude=false);
+  /// Returns prediction for point given by indx. This is used only as a means of directly comparing CPU and GPU outputs.
+  virtual NEWIMAGE::volume<float> PredictCPU(unsigned int indx, bool exclude=false);
+  /// Returns predictions for points given by indicies
+  virtual std::vector<NEWIMAGE::volume<float> > Predict(const std::vector<unsigned int>& indicies, bool exclude=false);
+  /// Returns input data for point given by indx
+  virtual NEWIMAGE::volume<float> InputData(unsigned int indx) const;
+  /// Returns input data for points given by indicies
+  virtual std::vector<NEWIMAGE::volume<float> > InputData(const std::vector<unsigned int>& indicies) const;
+  /// Returns variance of prediction for point given by indx.
+  virtual double PredictionVariance(unsigned int indx, bool exclude=false);
+  /// Returns measurement-error variance for point given by indx.
+  virtual double ErrorVariance(unsigned int indx) const;
+  /// Returns true if all data has been loaded
+  virtual bool IsPopulated() const { return(_pop); }
+  /// Indicates if it is ready to make predictions.
+  virtual bool IsValid() const EddyTry { return(IsPopulated() && _mc && _Kmat->IsValid()); } EddyCatch
+  /// Specify the # of points we plan to put into the predictor.
+  virtual void SetNoOfScans(unsigned int n);
+  /// Set a point given by indx. This function is thread safe as long as different threads set different points.
+  virtual void SetScan(const NEWIMAGE::volume<float>& scan,  // May be thread safe if used "sensibly"
+                       const DiffPara&                dp,
+                       unsigned int                   indx,
+		       unsigned int                   sess) { SetScan(scan,dp,indx); }
+  /// Set a point given by indx. This function is thread safe as long as different threads set different points.
+  virtual void SetScan(const NEWIMAGE::volume<float>& scan,  // May be thread safe if used "sensibly"
+		       const DiffPara&                dp,
+		       unsigned int                   indx);
+  /// Returns the number of hyperparameters for the model
+  virtual unsigned int NoOfHyperPar() const EddyTry { return(_Kmat->NoOfHyperPar()); } EddyCatch
+  /// Returns the hyperparameters for the model
+  virtual std::vector<double> GetHyperPar() const EddyTry { return(_Kmat->GetHyperPar()); } EddyCatch
+  /// Evaluates the model so as to make the predictor ready to make predictions.
+  virtual void EvaluateModel(const NEWIMAGE::volume<float>& mask, float fwhm, bool verbose=false);
+  /// Evaluates the model so as to make the predictor ready to make predictions.
+  virtual void EvaluateModel(const NEWIMAGE::volume<float>& mask, bool verbose=false) EddyTry { EvaluateModel(mask,0.0,verbose); } EddyCatch
+  /// Writes internal content to disk for debug purposes
+  virtual void WriteImageData(const std::string& fname) const;
+  virtual void WriteMetaData(const std::string& fname) const EddyTry {
+    if (!IsPopulated()) throw EddyException("DiffusionGP::WriteMetaData: Not yet fully populated"); _Kmat->Write(fname);
+  } EddyCatch
+  virtual void Write(const std::string& fname) const EddyTry { WriteImageData(fname); WriteMetaData(fname); } EddyCatch
+private:
+  static const unsigned int nvoxhp = 500;
+
+  std::vector<std::shared_ptr<NEWIMAGE::volume<float> > > _sptrs;   // Pointers to the scans
+  std::shared_ptr<KMatrix>                                _Kmat;    // K-matrix for Gaussian Process
+  std::vector<DiffPara>                                   _dpars;   // Diffusion parameters
+  std::shared_ptr<HyParEstimator>                         _hpe;     // Hyperparameter estimator
+  bool                                                    _pop;     // True if all data is present
+  bool                                                    _mc;      // True if data has been mean-corrected
+  std::vector<std::shared_ptr<NEWIMAGE::volume<float> > > _mptrs;   // Pointers to mean images
+
+  void mean_correct(const std::vector<std::vector<unsigned int> >& mi);
+  bool is_populated() const;
+  void predict_image_cpu(unsigned int             indx,
+			 bool                     excl,
+			 const NEWMAT::RowVector& pv,
+			 NEWIMAGE::volume<float>& ima) const;
+  void predict_images_cpu(// Input
+			  const std::vector<unsigned int>&       indicies,
+			  bool                                   exclude,
+			  const std::vector<NEWMAT::RowVector>&  pvecs,
+			  // Output
+			  std::vector<NEWIMAGE::volume<float> >& pi) const;
+  #ifdef COMPILE_GPU
+  void predict_image_gpu(unsigned int             indx,
+			 bool                     excl,
+			 const NEWMAT::RowVector& pv,
+			 NEWIMAGE::volume<float>& ima) const;
+  void predict_images_gpu(// Input
+			  const std::vector<unsigned int>&       indicies,
+			  bool                                   exclude,
+			  const std::vector<NEWMAT::RowVector>&  pvecs,
+			  // Output
+			  std::vector<NEWIMAGE::volume<float> >& pi) const;
+  #endif
+  unsigned int which_mean(unsigned int indx) const;
+
+  bool get_y(// Input
+	     unsigned int           i,
+	     unsigned int           j,
+	     unsigned int           k,
+	     unsigned int           indx,
+	     bool                   exclude,
+	     // Output
+	     NEWMAT::ColumnVector&  y) const;
+};
+
+} // End namespace EDDY
+
+#endif // End #ifndef DiffusionGP_h

ECModels.cpp

+// Declarations of classes that implements a hirearchy
+// of models for fields from eddy currents induced by
+// diffusion gradients.
+//
+// ECModels.cpp
+//
+// Jesper Andersson, FMRIB Image Analysis Group
+//
+// Copyright (C) 2011 University of Oxford
+//
+
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <time.h>
+#include "armawrap/newmat.h"
+#ifndef EXPOSE_TREACHEROUS
+#define EXPOSE_TREACHEROUS           // To allow us to use .set_sform etc
+#endif
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "topup/topup_file_io.h"
+#include "EddyHelperClasses.h"
+#include "ECModels.h"
+
+using namespace EDDY;
+
+void DerivativeInstructions::SetSecondary(unsigned int               i,
+					  unsigned int               index,
+					  const SliceDerivModulator& sdm) EddyTry
+{
+  this->check_index(i,"SetSecondary");
+  if (_mt != ModulationType::SliceWise) throw EddyException("DerivativeInstructions::SetSecondary: Wrong modulator for slice object");
+  _scnd[i]._index = index;
+  _scnd[i]._slmod = sdm;
+  _set[i] = true;
+  return;
+} EddyCatch
+
+void DerivativeInstructions::SetSecondary(unsigned int                 i,
+					  unsigned int                 index,
+					  const SpatialDerivModulator& sdm) EddyTry
+{
+  this->check_index(i,"SetSecondary");
+  if (_mt != ModulationType::Spatial) throw EddyException("DerivativeInstructions::SetSecondary: Wrong modulator for spatial object");
+  _scnd[i]._index = index;
+  _scnd[i]._spmod = sdm;
+  _set[i] = true;
+  return;
+} EddyCatch
+
+EDDY::DerivativeInstructions ScanMovementModel::GetCompoundDerivInstructions(unsigned int                 indx,
+									     const EDDY::MultiBandGroups& mbg) const EddyTry
+{
+  if (indx>5) throw EddyException("ScanMovementModel::GetCompoundDerivInstructions: indx out of range");
+  float scale = (indx<3) ? 1e-2 : 1e-5;
+  EDDY::DerivativeInstructions di(indx*(_order+1),scale,ModulationType::SliceWise,_order);
+  NEWMAT::Matrix X = get_design(mbg.NGroups());
+  for (unsigned int i=0; i<_order; i++) {
+    std::vector<float> sw_wgt(mbg.NSlices()); // Slice wise weight
+    for (unsigned int j=0; j<mbg.NGroups(); j++) {
+      std::vector<unsigned int> satp = mbg.SlicesAtTimePoint(j);
+      for (unsigned int k=0; k<satp.size(); k++) sw_wgt[satp[k]] = X(j+1,i+2);
+    }
+    di.SetSecondary(i,di.GetPrimaryIndex()+1+i,EDDY::SliceDerivModulator(sw_wgt));
+  }
+  return(di);
+} EddyCatch
+
+NEWMAT::Matrix ScanMovementModel::ForwardMovementMatrix(const NEWIMAGE::volume<float>& scan) const EddyTry
+{
+  if (_order) return(TOPUP::MovePar2Matrix(get_zero_order_mp(),scan));
+  else return(TOPUP::MovePar2Matrix(_mp,scan));
+} EddyCatch
+
+NEWMAT::Matrix ScanMovementModel::ForwardMovementMatrix(const NEWIMAGE::volume<float>& scan,
+							unsigned int                   grp,
+							unsigned int                   ngrp) const EddyTry
+{
+  if (grp>=ngrp) throw EddyException("ScanMovementModel::ForwardMovementMatrix: grp has to be smaller than ngrp");
+  NEWMAT::ColumnVector gmp = get_gmp(grp,ngrp);
+  return(TOPUP::MovePar2Matrix(gmp,scan));
+} EddyCatch
+
+NEWMAT::Matrix ScanMovementModel::RestrictedForwardMovementMatrix(const NEWIMAGE::volume<float>&       scan,
+								  const std::vector<unsigned int>&     rindx) const EddyTry
+{
+  NEWMAT::ColumnVector rmp;
+  if (_order) rmp = get_zero_order_mp();
+  else rmp = _mp;
+  for (unsigned int i=0; i<rindx.size(); i++) {
+    if (rindx[i] > 5) throw EddyException("ScanMovementModel::RestrictedForwardMovementMatrix: rindx has to be less than 6");
+    else rmp(rindx[i]+1) = 0.0;
+  }
+  return(TOPUP::MovePar2Matrix(rmp,scan));
+} EddyCatch
+
+NEWMAT::Matrix ScanMovementModel::RestrictedForwardMovementMatrix(const NEWIMAGE::volume<float>&       scan,
+								  unsigned int                         grp,
+								  unsigned int                         ngrp,
+								  const std::vector<unsigned int>&     rindx) const EddyTry
+{
+  if (grp>=ngrp) throw EddyException("ScanMovementModel::RestrictedForwardMovementMatrix: grp has to be smaller than ngrp");
+  NEWMAT::ColumnVector rgmp = get_gmp(grp,ngrp);
+  for (unsigned int i=0; i<rindx.size(); i++) {
+    if (rindx[i] > 5) throw EddyException("ScanMovementModel::RestrictedForwardMovementMatrix: rindx has to be less than 6");
+    else rgmp(rindx[i]+1) = 0.0;
+  }
+  return(TOPUP::MovePar2Matrix(rgmp,scan));
+} EddyCatch
+
+NEWIMAGE::volume<float> LinearScanECModel::ECField(const NEWIMAGE::volume<float>& scan) const EddyTry
+{
+  NEWIMAGE::volume<float> field = scan;
+  field = _ep(4);
+  for (int k=0; k<field.zsize(); k++) {
+    float zcomp = _ep(3)*field.zdim()*(k-(field.zsize()-1)/2.0);
+    for (int j=0; j<field.ysize(); j++) {
+      float ycomp = _ep(2)*field.ydim()*(j-(field.ysize()-1)/2.0);
+      for (int i=0; i<field.xsize(); i++) {
+        field(i,j,k) += _ep(1)*field.xdim()*(i-(field.xsize()-1)/2.0) + ycomp + zcomp;
+      }
+    }
+  }
+  return(field);
+} EddyCatch
+
+EDDY::DerivativeInstructions LinearScanECModel::GetCompoundDerivInstructions(unsigned int                     indx,
+									     const std::vector<unsigned int>& pev) const EddyTry
+{
+  if (indx<1) { // This has the constant phase as primary and is modulated in two (non-PE) directions
+    float scale = this->GetDerivScale(3,true);
+    EDDY::DerivativeInstructions di(3,scale,ModulationType::Spatial,2);
+    std::vector<unsigned int> zmod{0,0,1}; // Modulation in z
+    di.SetSecondary(0,2,EDDY::SpatialDerivModulator(zmod));
+    if (pev[0]) { // If phase encode in x
+      std::vector<unsigned int> ymod{0,1,0}; // Modulation in y
+      di.SetSecondary(1,1,EDDY::SpatialDerivModulator(ymod));
+    }
+    else if (pev[1]) { // If phase encode in y
+      std::vector<unsigned int> xmod{1,0,0}; // Modulation in x
+      di.SetSecondary(1,0,EDDY::SpatialDerivModulator(xmod));
+    }
+    return(di);
+  }
+  else if (indx<2) { // This is the linear phase in the PE direction
+    unsigned int pedir = (pev[0]) ? 0 : 1;
+    float scale = this->GetDerivScale(pedir);
+    EDDY::DerivativeInstructions di(pedir,scale,ModulationType::Spatial,0);
+    return(di);
+  }
+  else throw EddyException("LinearScanECModel::GetCompoundDerivInstructions: indx out of range");
+} EddyCatch
+
+NEWIMAGE::volume<float> QuadraticScanECModel::ECField(const NEWIMAGE::volume<float>& scan) const EddyTry
+{
+  NEWIMAGE::volume<float> field = scan;
+  field = _ep(10); // DC offset
+  for (int k=0; k<field.zsize(); k++) {
+    double z = field.zdim()*(k-(field.zsize()-1)/2.0);
+    double zcomp = _ep(3)*z;
+    double z2comp = _ep(6)*z*z;
+    for (int j=0; j<field.ysize(); j++) {
+      double y = field.ydim()*(j-(field.ysize()-1)/2.0);
+      double ycomp = _ep(2)*y;
+      double y2comp = _ep(5)*y*y;
+      double yzcomp = _ep(9)*y*z;
+      for (int i=0; i<field.xsize(); i++) {
+        double x = field.xdim()*(i-(field.xsize()-1)/2.0);
+        double xcomp = _ep(1)*x;
+        double x2comp = _ep(4)*x*x;
+	double xycomp = _ep(7)*x*y;
+	double xzcomp = _ep(8)*x*z;
+        field(i,j,k) += xcomp + ycomp + zcomp + x2comp + y2comp + z2comp + xycomp + xzcomp + yzcomp;
+      }
+    }
+  }
+  return(field);
+} EddyCatch
+
+EDDY::DerivativeInstructions QuadraticScanECModel::GetCompoundDerivInstructions(unsigned int                     indx,
+										const std::vector<unsigned int>& pev) const EddyTry
+{
+  if (indx<1) { // This has the constant phase and is modulated linearly and quadratically in two directions.
+    float scale = this->GetDerivScale(9,true);
+    EDDY::DerivativeInstructions di(9,scale,ModulationType::Spatial,5);
+    std::vector<unsigned int> mod{0,0,1}; // Linear moulation in z
+    di.SetSecondary(0,2,EDDY::SpatialDerivModulator(mod));
+    mod = {0,0,2}; // Quadratic modulation in z
+    di.SetSecondary(1,5,EDDY::SpatialDerivModulator(mod));
+    if (pev[0]) { // If phase encode in x
+      mod = {0,1,0}; // Linear modulation in y
+      di.SetSecondary(2,1,EDDY::SpatialDerivModulator(mod));
+      mod = {0,2,0}; // Quadratic modulation in y
+      di.SetSecondary(3,4,EDDY::SpatialDerivModulator(mod));
+      mod = {0,1,1}; // Linear modulation in y _and_ z
+      di.SetSecondary(4,8,EDDY::SpatialDerivModulator(mod));
+    }
+    else if (pev[1]) { // If phase encode in y
+      mod = {1,0,0}; // Linear modulation in x
+      di.SetSecondary(2,0,EDDY::SpatialDerivModulator(mod));
+      mod = {2,0,0}; // Quadratic modulation in x
+      di.SetSecondary(3,3,EDDY::SpatialDerivModulator(mod));
+      mod = {1,0,1}; // Linear modulation in x _and_ z
+      di.SetSecondary(4,7,EDDY::SpatialDerivModulator(mod));
+    }
+    return(di);
+  }
+  else if (indx<2) { // This has the linear phase in the PE-direction and is modulated linearly in two directions
+    unsigned int pedir = (pev[0]) ? 0 : 1;
+    float scale = this->GetDerivScale(pedir);
+    EDDY::DerivativeInstructions di(pedir,scale,ModulationType::Spatial,2);
+    if (pev[0]) { // If phase encode in x
+      std::vector<unsigned int> mod{0,1,0}; // Linear moulation in y
+      di.SetSecondary(0,6,EDDY::SpatialDerivModulator(mod));
+      mod = {0,0,1}; // Linear modulation in z
+      di.SetSecondary(1,7,EDDY::SpatialDerivModulator(mod));
+    }
+    else if (pev[1]) { // If phase encode in y
+      std::vector<unsigned int> mod{1,0,0}; // Linear moulation in x
+      di.SetSecondary(0,6,EDDY::SpatialDerivModulator(mod));
+      mod = {0,0,1}; // Linear modulation in z
+      di.SetSecondary(1,8,EDDY::SpatialDerivModulator(mod));
+    }
+    return(di);
+  }
+  else if (indx<3) { // This has the quadratic phase in the PE-direction and is not used for any other derivatives
+    unsigned int pedir = (pev[0]) ? 0 : 1;
+    float scale = this->GetDerivScale(3+pedir);
+    EDDY::DerivativeInstructions di(3+pedir,scale,ModulationType::Spatial,0);
+    return(di);
+  }
+  else throw EddyException("QuadraticScanECModel::GetCompoundDerivInstructions: indx out of range");
+} EddyCatch
+
+NEWIMAGE::volume<float> CubicScanECModel::ECField(const NEWIMAGE::volume<float>& scan) const EddyTry
+{
+  NEWIMAGE::volume<float> field = scan;
+  field = _ep(20); // DC offset
+  for (int k=0; k<field.zsize(); k++) {
+    double z = field.zdim()*(k-(field.zsize()-1)/2.0);
+    double z2 = z*z;
+    double zcomp = _ep(3)*z;
+    double z2comp = _ep(6)*z2;
+    double z3comp = _ep(12)*z*z2;
+    for (int j=0; j<field.ysize(); j++) {
+      double y = field.ydim()*(j-(field.ysize()-1)/2.0);
+      double y2 = y*y;
+      double ycomp = _ep(2)*y;
+      double y2comp = _ep(5)*y2;
+      double y3comp = _ep(11)*y*y2;
+      double yzcomp = _ep(9)*y*z;
+      double y2zcomp = _ep(17)*y2*z;
+      double yz2comp = _ep(19)*y*z2;
+      for (int i=0; i<field.xsize(); i++) {
+        double x = field.xdim()*(i-(field.xsize()-1)/2.0);
+	double x2 = x*x;
+        double xcomp = _ep(1)*x;
+        double x2comp = _ep(4)*x2;
+        double x3comp = _ep(10)*x*x2;
+	double xycomp = _ep(7)*x*y;
+	double xzcomp = _ep(8)*x*z;
+	double x2ycomp = _ep(13)*x2*y;
+	double x2zcomp = _ep(14)*x2*z;
+	double xyzcomp = _ep(15)*x*y*z;
+	double xy2comp = _ep(16)*x*y2;
+	double xz2comp = _ep(18)*x*z2;
+        field(i,j,k) += xcomp + ycomp + zcomp + x2comp + y2comp + z2comp + xycomp + xzcomp + yzcomp;
+        field(i,j,k) += x3comp + y3comp + z3comp + x2ycomp + x2zcomp + xyzcomp + xy2comp + y2zcomp + xz2comp + yz2comp;
+      }
+    }
+  }
+  return(field);
+} EddyCatch
+
+EDDY::DerivativeInstructions CubicScanECModel::GetCompoundDerivInstructions(unsigned int                     indx,
+									    const std::vector<unsigned int>& pev) const EddyTry
+{
+  if (indx<1) { // This has the constant phase and is modulated linearly, quadratically and cubically in two directions.
+    float scale = this->GetDerivScale(19,true);
+    EDDY::DerivativeInstructions di(19,scale,ModulationType::Spatial,9);
+    std::vector<unsigned int> mod{0,0,1}; // Linear modulation in z
+    di.SetSecondary(0,2,EDDY::SpatialDerivModulator(mod));
+    mod = {0,0,2}; // Quadratic modulation in z
+    di.SetSecondary(1,5,EDDY::SpatialDerivModulator(mod));
+    mod = {0,0,3}; // Cubic modulation in z
+    di.SetSecondary(2,11,EDDY::SpatialDerivModulator(mod));
+    if (pev[0]) { // If phase encode in x
+      mod = {0,1,0}; // Linear modulation in y
+      di.SetSecondary(3,1,EDDY::SpatialDerivModulator(mod));
+      mod = {0,2,0}; // Quadratic modulation in y
+      di.SetSecondary(4,4,EDDY::SpatialDerivModulator(mod));
+      mod = {0,3,0}; // Cubic modulation in y
+      di.SetSecondary(5,10,EDDY::SpatialDerivModulator(mod));
+      mod = {0,1,1}; // Linear modulation in y _and_ z
+      di.SetSecondary(6,8,EDDY::SpatialDerivModulator(mod));
+      mod = {0,2,1}; // Quadratic modulation in y and linear in z
+      di.SetSecondary(7,16,EDDY::SpatialDerivModulator(mod));
+      mod = {0,1,2}; // Linear modulation in y and quadratic in z
+      di.SetSecondary(8,18,EDDY::SpatialDerivModulator(mod));
+    }
+    else if (pev[1]) { // If phase encode in y
+      mod = {1,0,0}; // Linear modulation in x
+      di.SetSecondary(3,0,EDDY::SpatialDerivModulator(mod));
+      mod = {2,0,0}; // Quadratic modulation in x
+      di.SetSecondary(4,3,EDDY::SpatialDerivModulator(mod));
+      mod = {3,0,0}; // Cubic modulation in x
+      di.SetSecondary(5,9,EDDY::SpatialDerivModulator(mod));
+      mod = {1,0,1}; // Linear modulation in x _and_ z
+      di.SetSecondary(6,7,EDDY::SpatialDerivModulator(mod));
+      mod = {2,0,1}; // Quadratic modulation in x and linear in z
+      di.SetSecondary(7,13,EDDY::SpatialDerivModulator(mod));
+      mod = {1,0,2}; // Linear modulation in x and quadratic in z
+      di.SetSecondary(8,17,EDDY::SpatialDerivModulator(mod));
+    }
+    return(di);
+  }
+  else if (indx<2) { // This has the linear phase in the PE-direction and is modulated linearly and quadratically in two directions
+    unsigned int pedir = (pev[0]) ? 0 : 1;
+    float scale = this->GetDerivScale(pedir);
+    EDDY::DerivativeInstructions di(pedir,scale,ModulationType::Spatial,5);
+    if (pev[0]) { // If phase encode in x
+      std::vector<unsigned int> mod{0,1,0}; // Linear modulation in y
+      di.SetSecondary(0,6,EDDY::SpatialDerivModulator(mod));
+      mod = {0,0,1}; // Linear modulation in z
+      di.SetSecondary(1,7,EDDY::SpatialDerivModulator(mod));
+      mod = {0,2,0}; // Quadratic modulation in y
+      di.SetSecondary(2,15,EDDY::SpatialDerivModulator(mod));
+      mod ={0,1,1}; // Linear modulation in y _and_ z
+      di.SetSecondary(3,14,EDDY::SpatialDerivModulator(mod));
+      mod ={0,0,2}; // Quadratic modulation in z
+      di.SetSecondary(4,17,EDDY::SpatialDerivModulator(mod));
+    }
+    else if (pev[1]) { // If phase encode in y
+      std::vector<unsigned int> mod{1,0,0}; // Linear modulation in x
+      di.SetSecondary(0,6,EDDY::SpatialDerivModulator(mod));
+      mod ={0,0,1}; // Linear modulation in z
+      di.SetSecondary(1,8,EDDY::SpatialDerivModulator(mod));
+      mod ={2,0,0}; // Quadratic modulation in x
+      di.SetSecondary(2,12,EDDY::SpatialDerivModulator(mod));
+      mod ={1,0,1}; // Linear modulation in x _and_ z
+      di.SetSecondary(3,14,EDDY::SpatialDerivModulator(mod));
+      mod ={0,0,2}; // Quadratic modulation in z
+      di.SetSecondary(4,18,EDDY::SpatialDerivModulator(mod));
+    }
+    return(di);
+  }
+  else if (indx<3) { // This has the quadratic phase in the PE-direction and is used for linear modulation in two directions
+    unsigned int pedir = (pev[0]) ? 0 : 1;
+    float scale = this->GetDerivScale(3+pedir);
+    EDDY::DerivativeInstructions di(3+pedir,scale,ModulationType::Spatial,2);
+    if (pev[0]) { // If phase encode in x
+      std::vector<unsigned int> mod{0,1,0}; // Linear modulation in y
+      di.SetSecondary(0,12,EDDY::SpatialDerivModulator(mod));
+      mod = {0,0,1}; // Linear modulation in z
+      di.SetSecondary(1,13,EDDY::SpatialDerivModulator(mod));
+    }
+    else if (pev[1]) { // If phase encode in y
+      std::vector<unsigned int> mod{1,0,0}; // Linear modulation in x
+      di.SetSecondary(0,15,EDDY::SpatialDerivModulator(mod));
+      mod = {0,0,1}; // Linear modulation in z
+      di.SetSecondary(1,16,EDDY::SpatialDerivModulator(mod));
+    }
+    return(di);
+  }
+  else if (indx<4) { // This has the Cubic phase in the PE-direction and is not used for any other derivatives
+    unsigned int pedir = (pev[0]) ? 0 : 1;
+    float scale = this->GetDerivScale(9+pedir);
+    EDDY::DerivativeInstructions di(9+pedir,scale,ModulationType::Spatial,0);
+    return(di);
+  }
+  else throw EddyException("CubicScanECModel::GetCompoundDerivInstructions: indx out of range");
+} EddyCatch