splinterpolator.h

//
//  splinterpolator.h
//
// Jesper Andersson, FMRIB Image Analysis Group
//
// Copyright (C) 2008 University of Oxford
//
//     CCOPYRIGHT
//

#ifndef splinterpolator_h
#define splinterpolator_h

#include <vector>
#include <string>
#include <cmath>
#include <thread>
#include <iomanip>
#include "armawrap/newmat.h"
#include "utils/threading.h"
#include "miscmaths/miscmaths.h"

namespace SPLINTERPOLATOR {

enum ExtrapolationType {Zeros, Constant, Mirror, Periodic};

class SplinterpolatorException: public std::exception
{
public:
  SplinterpolatorException(const std::string& msg) noexcept : m_msg(std::string("Splinterpolator::") + msg) {}
  ~SplinterpolatorException() noexcept {}
  virtual const char *what() const noexcept { return(m_msg.c_str()); }
private:
  std::string m_msg;
};

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
// Class Splinterpolator:
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

template<class T>
class Splinterpolator
{
public:
  // Constructors
  Splinterpolator()
    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(1) {}
  Splinterpolator(const T *data_or_coefs,
                  const std::vector<unsigned int>& dim,
                  const std::vector<ExtrapolationType>& et,
                  unsigned int order=3,
                  bool copy_low_order=true,
                  Utilities::NoOfThreads nthr=Utilities::NoOfThreads(1),
                  double prec=1e-8,
                  bool data_are_coefs=false)
    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr._n)
  {
    common_construction(data_or_coefs,dim,order,prec,et,copy_low_order,data_are_coefs);
  }
  Splinterpolator(const T *data_or_coefs,
                  const std::vector<unsigned int>& dim,
                  ExtrapolationType et=Zeros,
                  unsigned int order=3,
                  bool copy_low_order=true,
                  Utilities::NoOfThreads nthr=Utilities::NoOfThreads(1),
                  double prec=1e-8,
                  bool data_are_coefs=false)
    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr._n)
  {
    std::vector<ExtrapolationType>   ett(dim.size(),et);
    common_construction(data_or_coefs,dim,order,prec,ett,copy_low_order,data_are_coefs);
  }
  // Copy construction. May be removed in future
  Splinterpolator(const Splinterpolator<T>& src)
    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0) { assign(src); }

  // Destructor
  ~Splinterpolator() { if(_own_coef) delete [] _coef; }

  // Assignment. May be removed in future
  Splinterpolator& operator=(const Splinterpolator& src)
  { if(_own_coef) delete [] _coef; assign(src); return(*this); }

  // Copy the spline coefficients into dest.
  void Copy(std::vector<T>& dest)
  {
    unsigned int N = 1;
    for (auto d : _dim) {
      N *= d;
    }
    dest.resize(N);

    auto coefs = coef_ptr();

    for (unsigned int i = 0; i < N; i++) {
       dest[i] = coefs[i];
    }
  }

  // Set new data in Splinterpolator.
  void Set(const T *data_or_coefs,
           const std::vector<unsigned int>& dim,
           const std::vector<ExtrapolationType>& et,
           unsigned int order=3,
           bool copy_low_order=true,
           double prec=1e-8,
           bool data_are_coefs=false)
  {
    if (_own_coef) delete [] _coef;
    common_construction(data_or_coefs,dim,order,prec,et,copy_low_order,data_are_coefs);
  }
  void Set(const T *data_or_coefs,
           const std::vector<unsigned int>& dim,
           ExtrapolationType et,
           unsigned int order=3,
           bool copy_low_order=true,
           double prec=1e-8,
           bool data_are_coefs=false)
  {
    std::vector<ExtrapolationType>  vet(dim.size(),Zeros);
    Set(data_or_coefs,dim,vet,order,copy_low_order,prec,data_are_coefs);
  }

  // Return interpolated value
  T operator()(const std::vector<float>&  coord) const;
  T operator()(double x, double y=0, double z=0, double t=0) const
  {
   if (!_valid) throw SplinterpolatorException("operator(): Cannot interpolate un-initialized object");
   if (_ndim>4 || (t && _ndim<4) || (z && _ndim<3) || (y && _ndim<2)) throw SplinterpolatorException("operator(): input has wrong dimensionality");
    double coord[5] = {x,y,z,t,0.0};
    return(static_cast<T>(value_at(coord)));
  }

  // Return interpolated value along with first derivative in one direction (useful for distortion correction)
  T operator()(const std::vector<float>& coord, unsigned int dd, T *dval) const;
  T operator()(double x, double y, double z, unsigned int dd, T *dval) const;
  T operator()(double x, double y, unsigned int dd, T *dval) const { return((*this)(x,y,0.0,dd,dval)); }
  T operator()(double x, T *dval) const { return((*this)(x,0.0,0.0,0,dval)); }

  // Return interpolated value along with selected derivatives
  T ValAndDerivs(const std::vector<float>& coord, const std::vector<unsigned int>& deriv, std::vector<T>& rderiv) const;
  T ValAndDerivs(const std::vector<float>& coord, std::vector<T>& rderiv) const
  {
    std::vector<unsigned int>   deriv(_ndim,1);
    return(ValAndDerivs(coord,deriv,rderiv));
  }
  T ValAndDerivs(double x, double y, double z, std::vector<T>& rderiv) const;

  // Return continous derivative at voxel centres (only works for order>1)
  T Deriv(const std::vector<unsigned int>& indx, unsigned int ddir) const;
  T Deriv1(const std::vector<unsigned int>& indx) const {return(Deriv(indx,0));}
  T Deriv2(const std::vector<unsigned int>& indx) const {return(Deriv(indx,1));}
  T Deriv3(const std::vector<unsigned int>& indx) const {return(Deriv(indx,2));}
  T Deriv4(const std::vector<unsigned int>& indx) const {return(Deriv(indx,3));}
  T Deriv5(const std::vector<unsigned int>& indx) const {return(Deriv(indx,4));}
  T DerivXYZ(unsigned int i, unsigned int j, unsigned int k, unsigned int dd) const;
  T DerivX(unsigned int i, unsigned int j, unsigned int k) const {return(DerivXYZ(i,j,k,0));}
  T DerivY(unsigned int i, unsigned int j, unsigned int k) const {return(DerivXYZ(i,j,k,1));}
  T DerivZ(unsigned int i, unsigned int j, unsigned int k) const {return(DerivXYZ(i,j,k,2));}
  void Grad3D(unsigned int i, unsigned int j, unsigned int k, T *xg, T *yg, T *zg) const;
  void Grad(const std::vector<unsigned int>& indx, std::vector<T>& grad) const;

  // Return continous addition (since previous voxel) of integral at voxel centres
  T IntX() const;
  T IntY() const;
  T IntZ() const;

  //
  // The "useful" functionality pretty much ends here.
  // Remaining functions are mainly for debugging/diagnostics.
  //
  unsigned int NDim() const { return(_ndim); }
  unsigned int Order() const { return(_order); }
  bool Valid() const { return(_valid); }
  ExtrapolationType Extrapolation(unsigned int dim) const
  {
    if (dim >= _ndim) throw SplinterpolatorException("Extrapolation: Invalid dimension");
    return(_et[dim]);
  }
  const std::vector<unsigned int>& Size() const { return(_dim); }
  unsigned int Size(unsigned int dim) const { if (dim > 4) return(0); else return(_dim[dim]);}
  T Coef(unsigned int x, unsigned int y=0, unsigned int z=0) const
  {
    std::vector<unsigned int> indx(3,0);
    indx[0] = x; indx[1] = y; indx[2] = z;
    return(Coef(indx));
  }
  T Coef(std::vector<unsigned int> indx) const;
  NEWMAT::ReturnMatrix CoefAsNewmatMatrix() const;
  NEWMAT::ReturnMatrix KernelAsNewmatMatrix(double sp=0.1, unsigned int deriv=0) const;

  //
  // Here we declare nested helper-class SplineColumn
  //
  class SplineColumn
  {
  public:
    // Constructor