Skip to content
Snippets Groups Projects
Commit ba42d0b6 authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

ENH: NEW "SoftZeros" extrapolation method, which sets out-of-bounds locations 

to zero, but still interpolates out-of-bounds, to allow a smooth transition
from in-to-out of bounds
parent 433217f8
No related branches found
No related tags found
1 merge request!17ENH: Allow `Splinterpolator` instances to be created from an existing set of spline coefficients, and extend extrapolation options
This commit is part of merge request !17. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
Showing
with 65 additions and 39 deletions
splinterpolator.h
+ 65
− 39
Edit View file @ ba42d0b6
......@@ -22,7 +22,14 @@
namespace SPLINTERPOLATOR {
enum ExtrapolationType {Zeros, Constant, Mirror, Periodic};
// Controls how the field behaves outside of the FOV:
// - Zeros: Out-of-bounds coordinates are set to zero
// - Constant: Out-of-bounds coordinates are set to the boundary voxels
// - Mirror: The field is mirrored in all directions
// - Periodic: The field is repeated in all directions
// - SoftZeros: The field boundary is set to zero, but out-of-bounds voxels
// are calculated via interpolation
enum ExtrapolationType {Zeros, Constant, Mirror, Periodic, SoftZeros};
class SplinterpolatorException: public std::exception
{
......@@ -657,27 +664,35 @@ 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 start_inds[5];
int inds[5];
unsigned int ni = 0;
const T *cptr = coef_ptr();
ni = get_start_indicies(coord,inds);
get_wgts(coord,inds,wgts);
ni = get_start_indicies(coord,start_inds);
get_wgts(coord,start_inds,wgts);
double val=0.0;
for (unsigned int m=0, me=(_ndim>4)?ni:1; m<me; m++) {
inds[4] = start_inds[4] + m;
for (unsigned int l=0, le=(_ndim>3)?ni:1; l<le; l++) {
inds[3] = start_inds[3] + l;
for (unsigned int k=0, ke=(_ndim>2)?ni:1; k<ke; k++) {
inds[2] = start_inds[2] + k;
double wgt1 = wgts[4][m]*wgts[3][l]*wgts[2][k];
unsigned int linear1 = indx2linear(inds[2]+k,inds[3]+l,inds[4]+m);
for (unsigned int j=0, je=(_ndim>1)?ni:1; j<je; j++) {
inds[1] = start_inds[1] + j;
double wgt2 = wgt1*wgts[1][j];
int linear2 = add2linear(linear1,inds[1]+j);
double *iiwgt=iwgt;
double *iiwgt = iwgt;
for (unsigned int i=0; i<ni; i++, iiwgt++) {
val += cptr[linear2+indx2indx(inds[0]+i,0)]*(*iiwgt)*wgt2;
inds[0] = start_inds[0] + i;
val += coef(inds) * (*iiwgt) * wgt2;
}
}
}
}
}
}
......@@ -696,9 +711,9 @@ double Splinterpolator<T>::value_at(const double *coord) const
template<class T>
double Splinterpolator<T>::value_and_derivatives_at(const double *coord,
const unsigned int *deriv,
double *dval)
const
const unsigned int *deriv,
double *dval)
const
{
if (should_be_zero(coord)) { memset(dval,0,n_nonzero(deriv)*sizeof(double)); return(0.0); }
......@@ -708,38 +723,47 @@ const
double *dwgts[] = {diwgt, djwgt, dkwgt, dlwgt, dmwgt};
double dwgt1[5];
double dwgt2[5];
int start_inds[5];
int inds[5];
unsigned int dd[5];
unsigned int nd = 0;
unsigned int ni = 0;
const T *cptr = coef_ptr();
ni = get_start_indicies(coord,inds);
get_wgts(coord,inds,wgts);
get_dwgts(coord,inds,deriv,dwgts);
ni = get_start_indicies(coord,start_inds);
get_wgts(coord,start_inds,wgts);
get_dwgts(coord,start_inds,deriv,dwgts);
for (unsigned int i=0; i<_ndim; i++) if (deriv[i]) { dd[nd] = i; dval[nd++] = 0.0; }
double val=0.0;
for (unsigned int m=0, me=(_ndim>4)?ni:1; m<me; m++) {
inds[4] = start_inds[4] + m;
for (unsigned int l=0, le=(_ndim>3)?ni:1; l<le; l++) {
inds[3] = start_inds[3] + l;
for (unsigned int k=0, ke=(_ndim>2)?ni:1; k<ke; k++) {
inds[2] = start_inds[2] + k;
double wgt1 = wgts[4][m]*wgts[3][l]*wgts[2][k];
get_dwgt1(wgts,dwgts,dd,nd,k,l,m,wgt1,dwgt1);
unsigned int linear1 = indx2linear(inds[2]+k,inds[3]+l,inds[4]+m);
for (unsigned int j=0, je=(_ndim>1)?ni:1; j<je; j++) {
inds[1] = start_inds[1] + j;
double wgt2 = wgt1*wgts[1][j];
for (unsigned int d=0; d<nd; d++) dwgt2[d] = (dd[d]==1) ? dwgt1[d]*dwgts[1][j] : dwgt1[d]*wgts[1][j];
int linear2 = add2linear(linear1,inds[1]+j);
double *iiwgt=iwgt;
for (unsigned int i=0; i<ni; i++, iiwgt++) {
double c = cptr[linear2+indx2indx(inds[0]+i,0)];
inds[0] = start_inds[0] + i;
double c = coef(inds);
val += c*(*iiwgt)*wgt2;
for (unsigned int d=0; d<nd; d++) {
double add = (dd[d]==0) ? c*diwgt[i]*dwgt2[d] : c*(*iiwgt)*dwgt2[d];
dval[d] += add;
}
}
}
}
}
}
}
}
}
......@@ -1200,7 +1224,7 @@ inline unsigned int Splinterpolator<T>::indx2indx(int indx, unsigned int d) cons
if (indx < 0) indx = 0;
else if (indx >= dim) indx = dim-1;
}
else if (_et[d] == Zeros || _et[d] == Mirror) {
else if (_et[d] == Zeros || _et[d] == SoftZeros || _et[d] == Mirror) {
while (indx < 0) indx = 2*dim*((indx+1)/dim) - 1 - indx;
while (indx >= dim) indx = 2*dim*(indx/dim) - 1 - indx;
}
......@@ -1323,36 +1347,38 @@ T Splinterpolator<T>::coef(int *indx) const
for (unsigned int i=0; i<_ndim; i++) {
if (indx[i] < 0) {
switch (_et[i]) {
case SoftZeros:
case Zeros:
return(static_cast<T>(0));
return(static_cast<T>(0));
case Constant:
indx[i] = 0;
break;
indx[i] = 0;
break;
case Mirror:
indx[i] = 1-indx[i];
break;
indx[i] = 1-indx[i];
break;
case Periodic:
indx[i] = _dim[i]+indx[i];
break;
indx[i] = _dim[i]+indx[i];
break;
default:
break;
break;
}
}
else if (indx[i] >= static_cast<int>(_dim[i])) {
switch (_et[i]) {
case SoftZeros:
case Zeros:
return(static_cast<T>(0));
return(static_cast<T>(0));
case Constant:
indx[i] = _dim[i]-1;
break;
indx[i] = _dim[i]-1;
break;
case Mirror:
indx[i] = 2*_dim[i]-indx[i]-1;
break;
indx[i] = 2*_dim[i]-indx[i]-1;
break;
case Periodic:
indx[i] = indx[i]-_dim[i];
break;
indx[i] = indx[i]-_dim[i];
break;
default:
break;
break;
}
}
}
......@@ -1648,7 +1674,7 @@ double Splinterpolator<T>::SplineColumn::init_fwd_sweep(double z, ExtrapolationT
double z2i=z;
for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) iv += z2i * *ptr;
}
// et == Constant || et == Zeros
// et == Constant || et == Zeros || et == SoftZeros
else {
double *ptr=&_col[0];
double z2i=z;
......@@ -1683,7 +1709,7 @@ double Splinterpolator<T>::SplineColumn::init_bwd_sweep(double z, double lv, Ext
else if (et == Mirror) {
iv = -z/(1.0-z*z) * (2.0*_col[_sz-1] - lv);
}
// et == Constant || et == Zeros
// et == Constant || et == Zeros || et == SoftZeros
else {
iv = z / (z - 1) * _col[_sz-1];
}
......
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