Skip to content
Snippets Groups Projects
Commit b5d7de25 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 dc101478
No related branches found
No related tags found
No related merge requests found
Pipeline #25441 skipped
Stage: fsl-ci-build
Hide whitespace changes
Inline Side-by-side
Showing
with 65 additions and 39 deletions
splinterpolator.h
+ 65
− 39
View file @ b5d7de25
...@@ -22,7 +22,14 @@ ...@@ -22,7 +22,14 @@
namespace SPLINTERPOLATOR { 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 class SplinterpolatorException: public std::exception
{ {
...@@ -657,27 +664,35 @@ double Splinterpolator<T>::value_at(const double *coord) const ...@@ -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 iwgt[8], jwgt[8], kwgt[8], lwgt[8], mwgt[8];
double *wgts[] = {iwgt, jwgt, kwgt, lwgt, mwgt}; double *wgts[] = {iwgt, jwgt, kwgt, lwgt, mwgt};
int start_inds[5];
int inds[5]; int inds[5];
unsigned int ni = 0; unsigned int ni = 0;
const T *cptr = coef_ptr(); const T *cptr = coef_ptr();
ni = get_start_indicies(coord,inds); ni = get_start_indicies(coord,start_inds);
get_wgts(coord,inds,wgts); get_wgts(coord,start_inds,wgts);
double val=0.0; double val=0.0;
for (unsigned int m=0, me=(_ndim>4)?ni:1; m<me; m++) { 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++) { 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++) { 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]; 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++) { 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]; 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++) { 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 ...@@ -696,9 +711,9 @@ double Splinterpolator<T>::value_at(const double *coord) const
template<class T> template<class T>
double Splinterpolator<T>::value_and_derivatives_at(const double *coord, double Splinterpolator<T>::value_and_derivatives_at(const double *coord,
const unsigned int *deriv, const unsigned int *deriv,
double *dval) double *dval)
const const
{ {
if (should_be_zero(coord)) { memset(dval,0,n_nonzero(deriv)*sizeof(double)); return(0.0); } if (should_be_zero(coord)) { memset(dval,0,n_nonzero(deriv)*sizeof(double)); return(0.0); }
...@@ -708,38 +723,47 @@ const ...@@ -708,38 +723,47 @@ const
double *dwgts[] = {diwgt, djwgt, dkwgt, dlwgt, dmwgt}; double *dwgts[] = {diwgt, djwgt, dkwgt, dlwgt, dmwgt};
double dwgt1[5]; double dwgt1[5];
double dwgt2[5]; double dwgt2[5];
int start_inds[5];
int inds[5]; int inds[5];
unsigned int dd[5]; unsigned int dd[5];
unsigned int nd = 0; unsigned int nd = 0;
unsigned int ni = 0; unsigned int ni = 0;
const T *cptr = coef_ptr(); const T *cptr = coef_ptr();
ni = get_start_indicies(coord,inds); ni = get_start_indicies(coord,start_inds);
get_wgts(coord,inds,wgts); get_wgts(coord,start_inds,wgts);
get_dwgts(coord,inds,deriv,dwgts); 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; } for (unsigned int i=0; i<_ndim; i++) if (deriv[i]) { dd[nd] = i; dval[nd++] = 0.0; }
double val=0.0; double val=0.0;
for (unsigned int m=0, me=(_ndim>4)?ni:1; m<me; m++) { 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++) { 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++) { 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]; double wgt1 = wgts[4][m]*wgts[3][l]*wgts[2][k];
get_dwgt1(wgts,dwgts,dd,nd,k,l,m,wgt1,dwgt1); 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++) { 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]; 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]; 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; double *iiwgt=iwgt;
for (unsigned int i=0; i<ni; i++, iiwgt++) { 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; val += c*(*iiwgt)*wgt2;
for (unsigned int d=0; d<nd; d++) { for (unsigned int d=0; d<nd; d++) {
double add = (dd[d]==0) ? c*diwgt[i]*dwgt2[d] : c*(*iiwgt)*dwgt2[d]; double add = (dd[d]==0) ? c*diwgt[i]*dwgt2[d] : c*(*iiwgt)*dwgt2[d];
dval[d] += add; dval[d] += add;
} }
} }
} }
} }
} }
} }
...@@ -1195,7 +1219,7 @@ inline unsigned int Splinterpolator<T>::indx2indx(int indx, unsigned int d) cons ...@@ -1195,7 +1219,7 @@ inline unsigned int Splinterpolator<T>::indx2indx(int indx, unsigned int d) cons
if (indx < 0) indx = 0; if (indx < 0) indx = 0;
else if (indx >= dim) indx = dim-1; 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 < 0) indx = 2*dim*((indx+1)/dim) - 1 - indx;
while (indx >= dim) indx = 2*dim*(indx/dim) - 1 - indx; while (indx >= dim) indx = 2*dim*(indx/dim) - 1 - indx;
} }
...@@ -1318,36 +1342,38 @@ T Splinterpolator<T>::coef(int *indx) const ...@@ -1318,36 +1342,38 @@ T Splinterpolator<T>::coef(int *indx) const
for (unsigned int i=0; i<_ndim; i++) { for (unsigned int i=0; i<_ndim; i++) {
if (indx[i] < 0) { if (indx[i] < 0) {
switch (_et[i]) { switch (_et[i]) {
case SoftZeros:
case Zeros: case Zeros:
return(static_cast<T>(0)); return(static_cast<T>(0));
case Constant: case Constant:
indx[i] = 0; indx[i] = 0;
break; break;
case Mirror: case Mirror:
indx[i] = 1-indx[i]; indx[i] = 1-indx[i];
break; break;
case Periodic: case Periodic:
indx[i] = _dim[i]+indx[i]; indx[i] = _dim[i]+indx[i];
break; break;
default: default:
break; break;
} }
} }
else if (indx[i] >= static_cast<int>(_dim[i])) { else if (indx[i] >= static_cast<int>(_dim[i])) {
switch (_et[i]) { switch (_et[i]) {
case SoftZeros:
case Zeros: case Zeros:
return(static_cast<T>(0)); return(static_cast<T>(0));
case Constant: case Constant:
indx[i] = _dim[i]-1; indx[i] = _dim[i]-1;
break; break;
case Mirror: case Mirror:
indx[i] = 2*_dim[i]-indx[i]-1; indx[i] = 2*_dim[i]-indx[i]-1;
break; break;
case Periodic: case Periodic:
indx[i] = indx[i]-_dim[i]; indx[i] = indx[i]-_dim[i];
break; break;
default: default:
break; break;
} }
} }
} }
...@@ -1643,7 +1669,7 @@ double Splinterpolator<T>::SplineColumn::init_fwd_sweep(double z, ExtrapolationT ...@@ -1643,7 +1669,7 @@ double Splinterpolator<T>::SplineColumn::init_fwd_sweep(double z, ExtrapolationT
double z2i=z; double z2i=z;
for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) iv += z2i * *ptr; 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 { else {
double *ptr=&_col[0]; double *ptr=&_col[0];
double z2i=z; double z2i=z;
...@@ -1678,7 +1704,7 @@ double Splinterpolator<T>::SplineColumn::init_bwd_sweep(double z, double lv, Ext ...@@ -1678,7 +1704,7 @@ double Splinterpolator<T>::SplineColumn::init_bwd_sweep(double z, double lv, Ext
else if (et == Mirror) { else if (et == Mirror) {
iv = -z/(1.0-z*z) * (2.0*_col[_sz-1] - lv); iv = -z/(1.0-z*z) * (2.0*_col[_sz-1] - lv);
} }
// et == Constant || et == Zeros // et == Constant || et == Zeros || et == SoftZeros
else { else {
iv = z / (z - 1) * _col[_sz-1]; 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