Skip to content
Snippets Groups Projects
Commit 776a361e authored by Jesper Andersson's avatar Jesper Andersson
Browse files

Checking in for backup, though still not quite there yet

parent 4b1f6a5d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 641 additions and 40 deletions
splinterpolator.h
+ 641
− 40
View file @ 776a361e
......@@ -29,30 +29,74 @@ public:
SplinterpolatorException(const std::string& msg) throw(): m_msg(msg) {}
virtual const char *what() const throw() {
return string("FslSplines::" + m_msg).c_str();
return string("Splinterpolator::" + m_msg).c_str();
}
~SplinterpolatorException() throw() {}
};
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
// Class Splinterpolator:
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
template<class T>
class Splinterpolator
{
public:
// Constructors
Splinterpolator() : _valid(false) {}
Splinterpolator(const T *data, const std::vector<unsigned int>& dim, unsigned int order=3, std::vector<ExtrapolationType> et(3,Zeros))
Splinterpolator(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3, double prec=1e-8)
{
common_construction(data,dim,order,et);
common_construction(data,dim,order,prec,et);
}
Splinterpolator(const T *data, const std::vector<unsigned int>& dim, unsigned int order, ExtrapolationType et)
Splinterpolator(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et=Zeros, unsigned int order=3, double prec=1e-8)
{
std::vector<ExtrapolationType> ett(dim.size(),et);
common_construction(data,dim,order,ett);
common_construction(data,dim,order,prec,ett);
}
// Destructor
~Splinterpolator() { if(_valid) delete [] _coef; }
// Set new data in Splinterpolator.
void Set(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3)
{
if (_valid) delete [] _coef;
common_construction(data,dim,order,et);
}
void Set(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et, unsigned int order=3)
{
std::vector<ExtrapolationType> vet(dim.size(),Zeros);
Set(data,dim,vet,order);
}
// 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
{
double coord[5] = {x,y,z,t,0.0};
return(value_at(coord));
}
//
// The "useful" functionality pretty much ends here.
// Remaining functions are mainly for debugging/diagnostics.
//
unsigned int NDim() const { return(_ndim); }
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
//
......@@ -64,12 +108,12 @@ public:
// Destructor
~SplineColumn() { delete [] _col; }
// Extract a column from a volume
Get(const T *dp)
void Get(const T *dp)
{
for (unsigned int i=0; i<_sz; i++, dp+=_step) _col[i] = static_cast<double>(*dp);
}
// Insert column into volume
Set(T *dp)
void Set(T *dp) const
{
T test = 1.5;
if (test == 1) { // If T is not float or double
......@@ -80,35 +124,50 @@ public:
}
}
// Deconvolve column
Deconv(unsigned int order, ExtrapolationType et) {};
void Deconv(unsigned int order, ExtrapolationType et, double prec);
private:
unsigned int _sz;
unsigned int _step;
double *_col;
unsigned int get_poles(unsigned int order, double *z, unsigned int *sf) const;
double init_bwd_sweep(double z, double lv, ExtrapolationType et, double prec) const;
double init_fwd_sweep(double z, ExtrapolationType et, double prec) const;
SplineColumn(const SplineColumn& sc); // Don't allow copy-construction
SplineColumn& operator=(const SplineColumn& sc); // Dont allow assignment
};
//
// Here ends nested helper-class SplineColumn
//
void Set(const T *data, const std::vector<unsigned int>& dim, unsigned int order=3, std::vector<ExtrapolationType> et(3,Zeros))
{
if (_valid) delete [] _coef;
common_construction(data,dim,order,et);
}
private:
bool _valid; // Decides if neccessary information has been set or not
T *_coef; // Volume of spline coefficients
unsigned int _order; // Order of splines
unsigned int _ndim; // # of non-singleton dimensions
double _prec; // Precision when dealing with edges
std::vector<unsigned int> _dim; // Dimensions of data
std::vector<ExtrapolationType> _et; // How to do extrapolation
//
// Private helper-functions
//
void common_construction(const T *data, const std::vector<unsigned int>& dim, unsigned int order, const std::vector<ExtrapolationType>& et);
void common_construction(const T *data, const std::vector<unsigned int>& dim, unsigned int order, double prec, const std::vector<ExtrapolationType>& et);
void calc_coef(const T *data);
void deconv_along(unsigned int dim);
T coef(int *indx) const;
double value_at(const double *coord) const;
unsigned int get_start_indicies(const double *coord, int *sinds) const;
unsigned int get_wgts(const double *coord, const int *sinds, double **wgts) const;
double get_wgt(double x) const;
double get_dwgt(double x) const;
std::pair<double,double> range() const;
bool odd(unsigned int i) const {return(static_cast<bool>(i%2));}
bool even(unsigned int i) const {return(!odd(i));}
//
// Disallowed member functions
//
......@@ -116,20 +175,423 @@ private:
Splinterpolator& operator=(const Splinterpolator& s); // Don't allow assignment
};
/////////////////////////////////////////////////////////////////////
//
// Here starts public member functions for Splinterpolator
//
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
//
// Returns the value of the coefficient given by indx (zero-offset)
//
/////////////////////////////////////////////////////////////////////
template<class T>
T Splinterpolator<T>::Coef(std::vector<unsigned int> indx) const
{
if (!indx.size()) throw SplinterpolatorException("Coef: indx has zeros dimensions");
if (indx.size() > 5) throw SplinterpolatorException("Coef: indx has more than 5 dimensions");
for (unsigned int i=0; i<indx.size(); i++) if (indx[i] >= _dim[i]) throw SplinterpolatorException("Coef: indx out of range");
unsigned int lindx=indx[indx.size()-1];
for (int i=indx.size()-2; i>=0; i--) lindx = _dim[i]*lindx + indx[i];
return(_coef[lindx]);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the values of all coefficients as a Newmat matrix. If
// _ndim==1 it will return a row-vector, if _ndim==2 it will return
// a matrix, if _ndim==3 it will return a tiled matrix where the n
// first rows (where n is the number of rows in one slice) pertain to
// the first slice, the next n rows to the second slice, etc. And
// correspondingly for 4- and 5-D.
//
/////////////////////////////////////////////////////////////////////
template<class T>
NEWMAT::ReturnMatrix Splinterpolator<T>::CoefAsNewmatMatrix() const
{
NEWMAT::Matrix mat(_dim[1]*_dim[2]*_dim[3]*_dim[4],_dim[0]);
std::vector<unsigned int> cindx(5,0);
unsigned int r=0;
for (cindx[4]=0; cindx[4]<_dim[4]; cindx[4]++) {
for (cindx[3]=0; cindx[3]<_dim[3]; cindx[3]++) {
for (cindx[2]=0; cindx[2]<_dim[2]; cindx[2]++) {
for (cindx[1]=0; cindx[1]<_dim[1]; cindx[1]++, r++) {
for (cindx[0]=0; cindx[0]<_dim[0]; cindx[0]++) {
mat.element(r,cindx[0]) = Coef(cindx);
}
}
}
}
}
mat.Release();
return(mat);
}
/////////////////////////////////////////////////////////////////////
//
// Return the kernel matrix to verify its correctness.
//
/////////////////////////////////////////////////////////////////////
template<class T>
void Splinterpolator<T>::common_construction(const T *data, const std::vector<unsigned int>& dim, unsigned int order, const std::vector<ExtrapolationType>& et)
NEWMAT::ReturnMatrix Splinterpolator<T>::KernelAsNewmatMatrix(double sp, // Distance (in ksp) between points
unsigned int deriv) const // Derivative (only 0/1 implemented).
{
if (!dim.size()) throw SplinterpolatorException("Splinterpolator::common_construction: data has zeros dimensions");
if (!dim.size() > 5) throw SplinterpolatorException("Splinterpolator::common_construction: data cannot have more than 5 dimensions");
if (dim.size() != et.size()) throw SplinterpolatorException("Splinterpolator::common_construction: dim and et must have the same size");
for (unsigned int i=0; i<dim.size(); i++) if (!dim[i]) throw SplinterpolatorException("Splinterpolator::common_construction: data cannot have zeros size in any direction");
if (order > 7) throw SplinterpolatorException("Splinterpolator::common_construction: spline order must be lesst than 7");
if (!data) throw SplinterpolatorException("Splinterpolator::common_construction: zero data pointer");
if (deriv > 1) throw SplinterpolatorException("KernelAsNewmatMatrix: only 1st derivatives implemented");
std::pair<double,double> rng = range();
unsigned int i=0;
for (double x=rng.first; x<=rng.second; x+=sp, i++) ; // Intentional
NEWMAT::Matrix kernel(i,2);
for (double x=rng.first, i=0; x<=rng.second; x+=sp, i++) {
kernel.element(i,0) = x;
kernel.element(i,1) = (deriv) ? get_dwgt(x) : get_wgt(x);
}
kernel.Release();
return(kernel);
}
/////////////////////////////////////////////////////////////////////
//
// Here starts public member functions for SplineColumn
//
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
//
// This function implements the forward and backwards sweep
// as defined by equation 2.5 in Unsers paper:
//
// B-spline signal processing. II. Efficiency design and applications
//
/////////////////////////////////////////////////////////////////////
template<class T>
void Splinterpolator<T>::SplineColumn::Deconv(unsigned int order, ExtrapolationType et, double prec)
{
double z[3] = {0.0, 0.0, 0.0}; // Poles
unsigned int np = 0; // # of poles
unsigned int sf; // Scale-factor
np = get_poles(order,z,&sf);
for (unsigned int p=0; p<np; p++) {
_col[0] = init_fwd_sweep(z[p],et,prec);
double lv = _col[_sz-1];
// Forward sweep
double *ptr=&_col[1];
for (unsigned int i=1; i<_sz; i++, ptr++) *ptr += z[p] * *(ptr-1);
_col[_sz-1] = init_bwd_sweep(z[p],lv,et,prec);
// Backward sweep
ptr = &_col[_sz-2];
for (int i=_sz-2; i>=0; i--, ptr--) *ptr = z[p]*(*(ptr+1) - *ptr);
}
double *ptr=_col;
for (unsigned int i=0; i<_sz; i++, ptr++) *ptr *= sf;
}
/////////////////////////////////////////////////////////////////////
//
// Here starts private member functions for Splinterpolator
//
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
//
// Returns the interpolated value at location given by coord.
// coord must be a pointer to an array of indicies with _ndim
// values.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::value_at(const double *coord) const
{
double iwgt[8], jwgt[8], kwgt[8], lwgt[8], mwgt[8];
double *wgts[] = {iwgt, jwgt, kwgt, lwgt, mwgt};
int inds[5];
unsigned int ni = 0;
ni = get_start_indicies(coord,inds);
get_wgts(coord,inds,wgts);
double val=0.0;
for (int m=0, me=(_ndim>4)?ni:1; m<me; m++) {
for (int l=0, le=(_ndim>3)?ni:1; l<le; l++) {
for (int k=0, ke=(_ndim>2)?ni:1; k<ke; k++) {
double wgt1 = wgts[4][m]*wgts[3][l]*wgts[2][k];
for (int j=0, je=(_ndim>1)?ni:1; j<je; j++) {
double wgt2 = wgt1*wgts[1][j];
for (int i=0; i<static_cast<int>(ni); i++) {
int cindx[] = {inds[0]+i,inds[1]+j,inds[2]+k,inds[3]+l,inds[4]+m};
val += coef(cindx)*wgts[0][i]*wgt2;
}
}
}
}
}
return(val);
}
/////////////////////////////////////////////////////////////////////
//
// Returns (in sinds) the indicies of the first coefficient in all
// _ndim directions with a non-zero weight for the location given
// by coord. The caller is responsible for coord and sinds being
// valid pointers to arrays of 5 values.
// The return-value gives the total # of non-zero weights.
//
/////////////////////////////////////////////////////////////////////
template<class T>
unsigned int Splinterpolator<T>::get_start_indicies(const double *coord, int *sinds) const
{
unsigned int ni = _order+1;
for (unsigned int ts=1, i=0; i<dim.size(); i++) ts *= dim[i];
if (odd(ni)) {
for (unsigned int i=0; i<_ndim; i++) {
sinds[i] = static_cast<int>(coord[i]+0.5) - ni/2;
}
}
else {
for (unsigned int i=0; i<_ndim; i++) {
int ix = static_cast<int>(coord[i]+0.5);
if (ix < coord[i]) sinds[i] = ix - (ni-1)/2;
else sinds[i] = ix -ni/2;
}
}
for (unsigned int i=_ndim; i<5; i++) sinds[i] = 0;
return(ni);
}
/////////////////////////////////////////////////////////////////////
//
// Returns (in wgts) the weights for the coefficients given by sinds
// for the location given by coord.
//
/////////////////////////////////////////////////////////////////////
template<class T>
unsigned int Splinterpolator<T>::get_wgts(const double *coord, const int *sinds, double **wgts) const
{
unsigned int ni = _order+1;
for (unsigned int dim=0; dim<_ndim; dim++) {
for (unsigned int i=0; i<ni; i++) {
wgts[dim][i] = get_wgt(coord[dim]-(sinds[dim]+i));
}
}
for (unsigned int dim=_ndim; dim<5; dim++) wgts[dim][0] = 1.0;
return(ni);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the weight for a spline at coordinate x, where x is relative
// to the centre of the spline.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::get_wgt(double x) const
{
double val = 0.0;
double ax = abs(x); // Kernels all symmetric
switch (_order) {
case 0:
if (ax < 0.5) val = 1.0;
break;
case 1:
if (ax < 1) val = 1-ax;;
break;
case 2:
if (ax < 0.5) val = 0.75-ax*ax;
else if (ax < 1.5) val = 0.5*(1.5-ax)*(1.5-ax);
break;
case 3:
if (ax < 1) val = 2.0/3.0 + 0.5*ax*ax*(ax-2);
else if (ax < 2) { ax = 2-ax; val = (1.0/6.0)*(ax*ax*ax); }
break;
case 4:
if (ax < 0.5) { ax *= ax; val = (115.0/192.0) + ax*((2.0*ax-5.0)/8.0); }
else if (ax < 1.5) val = (55.0/96.0) + ax*(ax*(ax*((5.0-ax)/6.0) - 1.25) + 5.0/24.0);
else if (ax < 2.5) { ax -= 2.5; ax *= ax; val = (1.0/24.0)*ax*ax; }
break;
case 5:
if (ax < 1) { double xx = ax*ax; val = 0.55 + xx*(xx*((3.0-ax)/12.0) - 0.5); }
else if (ax < 2) val = 0.425 + ax*(ax*(ax*(ax*((ax-9.0)/24.0) + 1.25) - 1.75) + 0.625);
else if (ax < 3) { ax = 3-ax; double xx = ax*ax; val = (1.0/120.0)*ax*xx*xx; }
break;
case 6:
if (ax < 0.5) { ax *= ax; val = (5887.0/11520.0) + ax*(ax*((21.0-4.0*ax)/144.0) -77.0/192.0); }
else if (ax < 1.5) val = 7861.0/15360.0 + ax*(ax*(ax*(ax*(ax*((ax - 7.0)/48.0) + 0.328125) - 35.0/288.0) - 91.0/256.0) -7.0/768.0);
else if (ax < 2.5) val = 1379.0/7680.0 + ax*(ax*(ax*(ax*(ax*((14.0-ax)/120.0) - 0.65625) + 133.0/72.0) - 2.5703125) + 1267.0/960.0);
else if (ax < 3.5) { ax -= 3.5; ax *= ax*ax; val = (1.0/720.0) * ax*ax; }
break;
case 7:
if (ax < 1) { double xx = ax*ax; val = 151.0/315.0 + xx*(xx*(xx*((ax-4.0)/144.0) + 1.0/9.0) - 1.0/3.0); }
else if (ax < 2) val = 103.0/210.0 + ax*(ax*(ax*(ax*(ax*(ax*((12.0-ax)/240.0) -7.0/30.0) + 0.5) - 7.0/18.0) - 0.1) -7.0/90.0);
else if (ax < 3) val = ax*(ax*(ax*(ax*(ax*(ax*((ax-20.0)/720.0) + 7.0/30.0) - 19.0/18.0) + 49.0/18.0) - 23.0/6.0) + 217.0/90.0) - 139.0/630.0;
else if (ax < 4) { ax = 4-ax; double xxx=ax*ax*ax; val = (1.0/5040.0)*ax*xxx*xxx; }
break;
default:
throw SplinterpolatorException("get_wgt: invalid order spline");
break;
}
return(val);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the weight for the first derivative of a spline at
// coordinate x, where x is relative to the centre of the spline.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::get_dwgt(double x) const
{
double val = 0.0;
double ax = abs(x); // Kernels all anti-symmetric
int sign = (ax) ? x/ax : 1; // Arbitrary choice for when x=0
switch (_order) {
case 0:
throw SplinterpolatorException("get_dwgt: invalid order spline");
break;
case 1:
if (ax < 1) val = sign * -1;
break;
case 2:
if (ax < 0.5) val = sign * -2.0*ax;
else if (ax < 1.5) val = sign * (-1.5 + ax);
break;
case 3:
if (ax < 1) val = sign * (1.5*ax*ax - 2.0*ax);
else if (ax < 2) { ax = 2-ax; val = sign * -0.5*ax*ax; }
break;
case 4:
if (ax < 0.5) val = sign * (ax*ax*ax - 1.25*ax);
else if (ax < 1.5) val = sign * (5.0/24.0 - ax*(2.5 - ax*(2.5 - (2.0/3.0)*ax)));
else if (ax < 2.5) { ax -= 2.5; val = sign * (1.0/6.0)*ax*ax*ax; }
break;
case 5:
if (ax < 1) val = sign * ax*(ax*(ax*(1-(5.0/12.0)*ax)) - 1);
else if (ax < 2) val = sign * (0.625 - ax*(3.5 - ax*(3.75 - ax*(1.5 - (5.0/24.0)*ax))));
else if (ax < 3) { ax -= 3; ax = ax*ax; val = sign * (-1.0/24.0)*ax*ax; }
break;
case 6:
break;
case 7:
break;
default:
throw SplinterpolatorException("get_dwgt: invalid order spline");
break;
}
return(val);
}
template<class T>
std::pair<double,double> Splinterpolator<T>::range() const
{
std::pair<double,double> rng(0.0,0.0);
rng.second = static_cast<double>(_order+1.0)/2.0;
rng.first = - rng.second;
return(rng);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the value of the coefficient indexed by indx. Unlike the
// public Coef() this routine allows indexing outside the valid
// volume, returning values that are dependent on the extrapolation
// model when these are encountered.
//
// N.B. May change value of input index N.B.
//
/////////////////////////////////////////////////////////////////////
template<class T>
T Splinterpolator<T>::coef(int *indx) const
{
// First fix any outside-volume indicies
for (unsigned int i=0; i<_ndim; i++) {
if (indx[i] < 0) {
switch (_et[i]) {
case Zeros:
return(static_cast<T>(0));
break;
case Constant:
indx[i] = 0;
break;
case Mirror:
indx[i] = 1-indx[i];
break;
case Periodic:
indx[i] = _dim[i]+indx[i];
break;
default:
break;
}
}
else if (indx[i] >= static_cast<int>(_dim[i])) {
switch (_et[i]) {
case Zeros:
return(static_cast<T>(0));
break;
case Constant:
indx[i] = _dim[i]-1;
break;
case Mirror:
indx[i] = 2*_dim[i]-indx[i]-1;
break;
case Periodic:
indx[i] = indx[i]-_dim[i];
break;
default:
break;
}
}
}
// Now make linear index
unsigned int lindx=indx[_ndim-1];
for (int i=_ndim-2; i>=0; i--) lindx = _dim[i]*lindx + indx[i];
return(_coef[lindx]);
}
/////////////////////////////////////////////////////////////////////
//
// Takes care of the "common" tasks when constructing a
// Splinterpolator object. Called by constructors and by .Set()
//
/////////////////////////////////////////////////////////////////////
template<class T>
void Splinterpolator<T>::common_construction(const T *data, const std::vector<unsigned int>& dim, unsigned int order, double prec, const std::vector<ExtrapolationType>& et)
{
if (!dim.size()) throw SplinterpolatorException("common_construction: data has zeros dimensions");
if (!dim.size() > 5) throw SplinterpolatorException("common_construction: data cannot have more than 5 dimensions");
if (dim.size() != et.size()) throw SplinterpolatorException("common_construction: dim and et must have the same size");
for (unsigned int i=0; i<dim.size(); i++) if (!dim[i]) throw SplinterpolatorException("common_construction: data cannot have zeros size in any direction");
if (order > 7) throw SplinterpolatorException("common_construction: spline order must be lesst than 7");
if (!data) throw SplinterpolatorException("common_construction: zero data pointer");
unsigned int ts=1;
for (unsigned int i=0; i<dim.size(); i++) ts *= dim[i];
_coef = new T[ts];
_order = order;
_dim.resize(5); = dim;
_prec = prec;
_dim.resize(5);
_ndim = dim.size();
for (unsigned int i=0; i<5; i++) _dim[i] = (i < dim.size()) ? dim[i] : 1;
_et = et;
calc_coef(data);
......@@ -137,40 +599,52 @@ void Splinterpolator<T>::common_construction(const T *data, const std::vector<un
_valid = true;
}
/////////////////////////////////////////////////////////////////////
//
// Performs deconvolution, converting signal to spline coefficients.
//
/////////////////////////////////////////////////////////////////////
template<class T>
void Splinterpolator<T>::calc_coef(const T *data)
{
// Put the original data into _coef
//
for (unsigned int ts=1, i=0; i<dim.size(); i++) ts *= dim[i];
unsigned int ts=1;
for (unsigned int i=0; i<_dim.size(); i++) ts *= _dim[i];
memcpy(_coef,data,ts*sizeof(T));
if (_order < 2) return; // If nearest neighbour or linear, that's all we need
// Loop over all non-singleton dimensions and deconvolve along them
//
std::vector<unsigned int> tdim(dim.size()-1,0);
for (unsigned int cdir=0; cdir<dim.size(); cdir++) {
if (dim[cdir] > 1) deconv_along(cdir)
std::vector<unsigned int> tdim(_dim.size()-1,0);
for (unsigned int cdir=0; cdir<_dim.size(); cdir++) {
if (_dim[cdir] > 1) deconv_along(cdir);
}
return;
}
/////////////////////////////////////////////////////////////////////
//
// Performs deconvolution along one of the dimensions, visiting
// all points along the other dimensions.
//
/////////////////////////////////////////////////////////////////////
template<class T>
void Splinterpolator<T>::deconv_along(unsigned int dim);
void Splinterpolator<T>::deconv_along(unsigned int dim)
{
T *dp = _coef;
// Set up to reflect "missing" dimension
//
std::vector<unsigned int> rdim(4,1); // Sizes along remaining dimensions
std::vector<unsigned int> rstep(4,1); // Step-sizes (in "volume") of remaining dimensions
unsigned int mdim = 1; // Size along "missing" dimension
unsigned int mstep = 1; // Step-size along "missing" dimension
for (unsigned int i=0, j=0, s=1; i<5; i++) {
for (unsigned int i=0, j=0, ss=1; i<5; i++) {
if (i == dim) { // If it is our "missing" dimension
mdim = _dim(i);
mdim = _dim[i];
mstep = ss;
}
else {
......@@ -180,15 +654,16 @@ void Splinterpolator<T>::deconv_along(unsigned int dim);
ss *= _dim[i];
}
SplineColumn<T> col(mdim,mstep); // Column helps us do the job
for (unsigned int l=0; l<rdim[3]; l++, dp+=rstep[3]) {
for (unsigned int k=0; k<rdim[2]; k++, dp+=rstep[2]) {
for (unsigned int j=0; j<rdim[1]; j++, dp+=rstep[1]) {
SplineColumn col(mdim,mstep); // Column helps us do the job
for (unsigned int l=0; l<rdim[3]; l++) {
for (unsigned int k=0; k<rdim[2]; k++) {
for (unsigned int j=0; j<rdim[1]; j++) {
T *dp = _coef + l*rstep[3] + k*rstep[2] + j*rstep[1];
for (unsigned int i=0; i<rdim[0]; i++, dp+=rstep[0]) {
col.Get(dp); // Extract a column from the volume
col.Deconv(_order,_et[dim]); // Deconvolve it
col.Set(dp); // Put back the deconvolved column
col.Get(dp); // Extract a column from the volume
col.Deconv(_order,_et[dim],_prec); // Deconvolve it
col.Set(dp); // Put back the deconvolved column
}
}
}
......@@ -196,6 +671,132 @@ void Splinterpolator<T>::deconv_along(unsigned int dim);
return;
}
/////////////////////////////////////////////////////////////////////
//
// Here starts private member functions for SplineColumn
//
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
//
// This function returns the poles and scale-factors for splines
// of order 2--7. The values correspond to those found in
// table 1 in Unsers 1993 paper:
// B-spline signal processing. II. Efficiency design and applications
//
// The actual values have been taken from
// http://bigwww.epfl.ch/thevenaz/interpolation/coeff.c
//
/////////////////////////////////////////////////////////////////////
template<class T>
unsigned int Splinterpolator<T>::SplineColumn::get_poles(unsigned int order, double *z, unsigned int *sf) const
{
unsigned int np = 0; // # of poles
switch (order) {
case 2:
np = 1;
z[0] = 2.0*sqrt(2.0) - 3.0;
*sf = 8;
break;
case 3:
np = 1;
z[0] = sqrt(3.0) - 2.0;
*sf = 6;
break;
case 4:
np = 2;
z[0] = sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0;
z[1] = sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0;
*sf = 384;
break;
case 5:
np = 2;
z[0] = sqrt(135.0 / 2.0 - sqrt(17745.0 / 4.0)) + sqrt(105.0 / 4.0) - 13.0 / 2.0;
z[1] = sqrt(135.0 / 2.0 + sqrt(17745.0 / 4.0)) - sqrt(105.0 / 4.0) - 13.0 / 2.0;
*sf = 120;
break;
case 6:
np = 3;
z[0] = -0.48829458930304475513011803888378906211227916123938;
z[1] = -0.081679271076237512597937765737059080653379610398148;
z[2] = -0.0014141518083258177510872439765585925278641690553467;
*sf = 46080;
break;
case 7:
np = 3;
z[0] = -0.53528043079643816554240378168164607183392315234269;
z[1] = -0.12255461519232669051527226435935734360548654942730;
z[2] = -0.0091486948096082769285930216516478534156925639545994;
*sf = 5040;
break;
default:
throw SplinterpolatorException("SplineColumn::get_poles: invalid order of spline");
break;
}
return(np);
}
/////////////////////////////////////////////////////////////////////
//
// Initialises the first value for the forward sweep. The initialisation
// will always be an approximation (this is where the "infinite" in IIR
// breaks down) so the value will be calculated to a predefined precision.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::SplineColumn::init_fwd_sweep(double z, ExtrapolationType et, double prec) const
{
//
// Move logs away from here after debugging
//
unsigned int n = static_cast<unsigned int>((log(prec)/log(abs(z))) + 1.5);
n = (n > _sz) ? _sz : n;
double iv = _col[0];
if (et == Mirror) {
double *ptr=&_col[1];
double z2i=z;
for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) iv += z2i * *ptr;
}
else {
double *ptr=&_col[_sz-1];
double z2i=z;
for (unsigned int i=1; i<n; i++, ptr--, z2i*=z) iv += z2i * *ptr;
}
return(iv);
}
/////////////////////////////////////////////////////////////////////
//
// Initialises the first value for the backward sweep. The initialisation
// will always be an approximation (this is where the "infinite" in IIR
// breaks down) so the value will be calculated to a predefined precision.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::SplineColumn::init_bwd_sweep(double z, double lv, ExtrapolationType et, double prec) const
{
double iv = 0.0;
if (et == Mirror) {
iv = -z/(1.0-z*z) * (2.0*_col[_sz-1] - lv);
}
else {
unsigned int n = static_cast<unsigned int>((log(prec)/log(abs(z))) + 1.5);
n = (n > _sz) ? _sz : n;
iv = z * _col[_sz-1];
double z2i = z*z;
double *ptr=_col;
for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) {
iv += z2i * *ptr;
}
iv /= (z2i-1.0);
}
return(iv);
}
} // End namespace SPLINTERPOLATOR
#endif // End #ifndef splinterpolator.h
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment