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

Added methods to get continous derivatives on voxel centers

parent d4686fb3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 308 additions and 0 deletions
splinterpolator.h
+ 308
− 0
View file @ 2511fbb0
......@@ -102,6 +102,25 @@ public:
}
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.
......@@ -195,11 +214,17 @@ private:
unsigned int add2linear(unsigned int lin, int j) const;
double value_at(const double *coord) const;
double value_and_derivatives_at(const double *coord, const unsigned int *deriv, double *dval) const;
void derivatives_at_i(const unsigned int *indx, const unsigned int *deriv, double *dval) const;
unsigned int get_start_indicies(const double *coord, int *sinds) const;
unsigned int get_start_indicies_at_i(const unsigned int *indx, int *sinds) const;
unsigned int get_wgts(const double *coord, const int *sinds, double **wgts) const;
unsigned int get_wgts_at_i(const unsigned int *indx, const int *sinds, double **wgts) const;
unsigned int get_dwgts(const double *coord, const int *sinds, const unsigned int *deriv, double **dwgts) const;
unsigned int get_dwgts_at_i(const unsigned int *indx, const int *sinds, const unsigned int *deriv, double **dwgts) const;
double get_wgt(double x) const;
double get_wgt_at_i(int i) const;
double get_dwgt(double x) const;
double get_dwgt_at_i(int i) const;
void get_dwgt1(const double * const *wgts, const double * const *dwgts, const unsigned int *dd, unsigned int nd,
unsigned int k, unsigned int l, unsigned int m, double wgt1, double *dwgt1) const;
......@@ -339,6 +364,75 @@ T Splinterpolator<T>::ValAndDerivs(double x, double y, double z, std::vector<T>&
return(rval);
}
/////////////////////////////////////////////////////////////////////
//
// Routine that returns a 3D gradient at an integer location.
//
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
//
// Routine that returns a single derivative at an integer location.
//
/////////////////////////////////////////////////////////////////////
template<class T>
T Splinterpolator<T>::Deriv(const std::vector<unsigned int>& indx, unsigned int dd) const
{
if (!_valid) throw SplinterpolatorException("Deriv: Cannot take derivative of un-initialized object");
if (indx.size() != _ndim) SplinterpolatorException("Deriv: Input indx of wrong dimension");
if (dd > (_ndim-1)) throw SplinterpolatorException("Deriv: derivative specified for invalid direction");
double dval;
unsigned int lindx[5] = {0,0,0,0,0};
unsigned int deriv[5] = {0,0,0,0,0};
for (unsigned int i=0; i<_ndim; i++) lindx[i]=indx[i];
deriv[dd] = 1;
derivatives_at_i(lindx,deriv,&dval);
return(static_cast<T>(dval));
}
template<class T>
T Splinterpolator<T>::DerivXYZ(unsigned int i, unsigned int j, unsigned int k, unsigned int dd) const
{
if (!_valid) throw SplinterpolatorException("DerivXYZ: Cannot take derivative of un-initialized object");
if (_ndim!=3 || dd>2) throw SplinterpolatorException("DerivXYZ: Input has wrong dimensionality");
double dval;
unsigned int lindx[5] = {i,j,k,0,0};
unsigned int deriv[5] = {0,0,0,0,0};
deriv[dd] = 1;
derivatives_at_i(lindx,deriv,&dval);
return(static_cast<T>(dval));
}
template<class T>
void Splinterpolator<T>::Grad3D(unsigned int i, unsigned int j, unsigned int k, T *xg, T *yg, T *zg) const
{
if (!_valid) throw SplinterpolatorException("Grad3D: Cannot take derivative of un-initialized object");
if (_ndim != 3) SplinterpolatorException("Grad3D: Input of wrong dimension");
unsigned int lindx[5] = {i,j,k,0,0};
unsigned int deriv[5] = {1,1,1,0,0};
double dval[5] = {0.0,0.0,0.0,0.0,0.0};
derivatives_at_i(lindx,deriv,dval);
*xg=static_cast<T>(dval[0]); *yg=static_cast<T>(dval[1]); *zg=static_cast<T>(dval[2]);
return;
}
template<class T>
void Splinterpolator<T>::Grad(const std::vector<unsigned int>& indx, std::vector<T>& grad) const
{
if (!_valid) throw SplinterpolatorException("Grad: Cannot take derivative of un-initialized object");
if (indx.size() != _ndim || grad.size() != _ndim) SplinterpolatorException("Grad: Input indx or grad of wrong dimension");
unsigned int lindx[5] = {0,0,0,0,0};
unsigned int deriv[5] = {0,0,0,0,0};
double dval[5] = {0.0,0.0,0.0,0.0,0.0};
for (unsigned int i=0; i<_ndim; i++) { lindx[i]=indx[i]; deriv[i]=1; }
derivatives_at_i(lindx,deriv,dval);
for (unsigned int i=0; i<_ndim; i++) grad[i] = static_cast<T>(dval[i]);
return;
}
/////////////////////////////////////////////////////////////////////
//
// Returns the value of the coefficient given by indx (zero-offset)
......@@ -596,6 +690,57 @@ const
return(val);
}
template <class T>
void Splinterpolator<T>::derivatives_at_i(const unsigned int *indx,
const unsigned int *deriv,
double *dval)
const
{
double iwgt[8], jwgt[8], kwgt[8], lwgt[8], mwgt[8];
double *wgts[] = {iwgt, jwgt, kwgt, lwgt, mwgt};
double diwgt[8], djwgt[8], dkwgt[8], dlwgt[8], dmwgt[8];
double *dwgts[] = {diwgt, djwgt, dkwgt, dlwgt, dmwgt};
double dwgt1[5];
double dwgt2[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_at_i(indx,inds);
get_wgts_at_i(indx,inds,wgts);
get_dwgts_at_i(indx,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++) {
for (unsigned int l=0, le=(_ndim>3)?ni:1; l<le; l++) {
for (unsigned int k=0, ke=(_ndim>2)?ni:1; k<ke; 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++) {
// 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)];
// 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;
}
}
}
}
}
}
// return(val);
return;
}
/////////////////////////////////////////////////////////////////////
//
// Returns (in sinds) the indicies of the first coefficient in all
......@@ -628,6 +773,19 @@ unsigned int Splinterpolator<T>::get_start_indicies(const double *coord, int *si
return(ni);
}
// Does the same thing, but for integer (spot on voxel centre) index
template<class T>
unsigned int Splinterpolator<T>::get_start_indicies_at_i(const unsigned int *indx, int *sinds) const
{
unsigned int ni = (odd(_order)) ? _order : _order+1;
for (unsigned int i=0; i<_ndim; i++) {
sinds[i] = indx[i] - (_order/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
......@@ -650,6 +808,22 @@ unsigned int Splinterpolator<T>::get_wgts(const double *coord, const int *sinds,
return(ni);
}
// Same for integer (spot on voxel centre) index
template<class T>
unsigned int Splinterpolator<T>::get_wgts_at_i(const unsigned int *indx, const int *sinds, double **wgts) const
{
unsigned int ni = (odd(_order)) ? _order : _order+1;
for (unsigned int dim=0; dim<_ndim; dim++) {
for (unsigned int i=0; i<ni; i++) {
wgts[dim][i] = get_wgt_at_i(indx[dim]-(sinds[dim]+i));
}
}
for (unsigned int dim=_ndim; dim<5; dim++) wgts[dim][0] = 1.0;
return(ni);
}
template<class T>
unsigned int Splinterpolator<T>::get_dwgts(const double *coord, const int *sinds, const unsigned int *deriv, double **dwgts) const
{
......@@ -678,6 +852,140 @@ unsigned int Splinterpolator<T>::get_dwgts(const double *coord, const int *sinds
return(ni);
}
// Same for integer (spot on voxel centre) index
template<class T>
unsigned int Splinterpolator<T>::get_dwgts_at_i(const unsigned int *indx, const int *sinds, const unsigned int *deriv, double **dwgts) const
{
unsigned int ni = (odd(_order)) ? _order : _order+1;
for (unsigned int dim=0; dim<_ndim; dim++) {
if (deriv[dim]) {
switch (_order) {
case 0: case 1:
throw SplinterpolatorException("get_dwgts_at_i: invalid order spline");
break;
case 2: case 3: case 4: case 5: case 6: case 7:
for (unsigned int i=0; i<ni; i++) {
dwgts[dim][i] = get_dwgt_at_i(indx[dim]-(sinds[dim]+i));
}
break;
default:
throw SplinterpolatorException("get_dwgts_at_i: invalid order spline");
}
}
}
return(ni);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the weight for a spline at integer index i, where i is
// relative to the centre index of the spline.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::get_wgt_at_i(int i) const
{
double val = 0.0;
int ai = std::abs(i);
switch (_order) {
case 0: case 1:
val = (ai) ? 1.0 : 0.0;
break;
case 2:
if (!ai) val = 0.75;
else if (ai==1) val = 0.125;
break;
case 3:
if (!ai) val = 0.666666666666667;
else if (ai==1) val = 0.166666666666667;
break;
case 4:
if (!ai) val = 0.598958333333333;
else if (ai==1) val = 0.197916666666667;
else if (ai==2) val = 0.002604166666667;
break;
case 5:
if (!ai) val = 0.55;
else if (ai==1) val = 0.216666666666667;
else if (ai==2) val = 0.008333333333333;
break;
case 6:
if (!ai) val = 0.511024305555556;
else if (ai==1) val = 0.228797743055556;
else if (ai==2) val = 0.015668402777779;
else if (ai==3) val = 8.680555555555556e-05;
break;
case 7:
if (!ai) val = 0.479365079365079;
else if (ai==1) val = 0.236309523809524;
else if (ai==2) val = 0.023809523809524;
else if (ai==3) val = 1.984126984126984e-04;
break;
default:
throw SplinterpolatorException("get_wgt_at_i: invalid order spline");
break;
}
return(val);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the weight for the first derivative of a spline at integer
// index i, where i is relative to the centre index of the spline.
//
/////////////////////////////////////////////////////////////////////
template<class T>
double Splinterpolator<T>::get_dwgt_at_i(int i) const
{
double val = 0.0;
int ai = std::abs(i);
int sign = (ai) ? i/ai : 1;
switch (_order) {
case 0: case 1:
throw SplinterpolatorException("get_dwgt: invalid order spline");
break;
case 2:
if (!ai) val = 0.0;
else if (ai==1) val = sign * (-0.5);
break;
case 3:
if (!ai) val = 0.0;
else if (ai==1) val = sign * (-0.5);
break;
case 4:
if (!ai) val = 0.0;
else if (ai==1) val = sign * (-0.458333333333333);
else if (ai==2) val = sign * (-0.020833333333333);
break;
case 5:
if (!ai) val = 0.0;
else if (ai==1) val = sign * (-0.416666666666667);
else if (ai==2) val = sign * (-0.041666666666667);
break;
case 6:
if (!ai) val = 0.0;
else if (ai==1) val = sign * (-0.376302083333333);
else if (ai==2) val = sign * (-0.061458333333334);
else if (ai==3) val = sign * (-2.604166666666667e-04);
break;
case 7:
if (!ai) val = 0.0;
else if (ai==1) val = sign * (-0.340277777777778);
else if (ai==2) val = sign * (-0.077777777777778);
else if (ai==3) val = sign * (-0.001388888888889);
break;
default:
throw SplinterpolatorException("get_dwgt_at_i: invalid order spline");
break;
}
return(val);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the weight for a spline at coordinate x, where x is relative
......
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