## Piecewise cubic Hermite interpolating polynomial (PCHIP)

Average rating

The Baseline Spline package now includes an option for PCHIP splines. Here's the PCHIP code:

```// populate w_interp with interpolated values from piecewise cubic
// Hermite spline
function PCHIP(w_nodesX, w_nodesY, w_interp, [xwave])
wave w_nodesX, w_nodesY, w_interp
wave /Z xwave

if(numpnts(w_nodesX)==2) // linear interpolation
if(waveexists(xwave))
w_interp=w_nodesY[0]+(xwave[p]-w_nodesX[0])/(w_nodesX[1]-w_nodesX[0])*(w_nodesY[1]-w_nodesY[0])
else
w_interp=w_nodesY[0]+(x-w_nodesX[0])/(w_nodesX[1]-w_nodesX[0])*(w_nodesY[1]-w_nodesY[0])
endif
return 1
endif

wave w_d=setDerivatives(w_nodesX, w_nodesY)
if(waveexists(xwave))
w_interp=InterpolatePCHIP(w_nodesX, w_nodesY, w_d, xwave[p])
else
w_interp=InterpolatePCHIP(w_nodesX, w_nodesY, w_d, x)
endif
end

// The trick is to choose good values for the slope at each node. See:
// FN Fritsch and RE Carlson (1980) Monotone piecewise cubic
// interpolation. SIAM J. Numer. Anal. 17: 238. DOI:10.1137/0717021
// KW Brodlie (1980) A review of methods for curve and function drawing.
// In Mathematical Methods in Computer Graphics and Design, KW Brodlie,
// ed., Academic Press, London, pp. 1-37.
// FN Fritsch and J Butland (1984) A method for constructing local
// monotone piecewise cubic interpolants, SIAM Journal on Scientific and
// Statistical Computing 5: 300-304.
static function /WAVE setDerivatives(w_nodesX, w_nodesY)
wave w_nodesX, w_nodesY

duplicate /free w_nodesX, w_d, w_a
// w_d will be set to desired derivatives at node positions
// w_a will be used as weighting for harmonic means
w_d=0

make /free/n=(numpnts(w_nodesX)-1), w_m, w_delx
// w_m[i] is gradient between node i and i+1
// w_delx[i] and w_delx[i-1] are x offsets of nodes i+1 and i-1

w_m = (w_nodesY[p+1]-w_nodesY)/(w_nodesX[p+1]-w_nodesX)
w_delx = w_nodesX[p+1]-w_nodesX
variable pmax=numpnts(w_d)-2  // pmax is numpnts(w_m)-1 = numpnts(w_d)-2

w_a[1,pmax] = (1+(w_delx[p-1])/(w_delx[p-1]+w_delx))/3

// Brodlie modification of Butland formula (for same sign slopes)
w_d [1,pmax] = ( w_m!=0 && w_m[p-1]!=0 && (sign(w_m)==sign(w_m[p-1])) ) ? w_m[p-1]*w_m/(w_a*w_m[p-1]+(1-w_a)*w_m) : 0

// deal with ends
w_d[0]=setEndDerivative(w_delx[0], w_delx[1], w_m[0], w_m[1])
w_d[pmax+1]=setEndDerivative(w_delx[pmax], w_delx[pmax-1], w_m[pmax], w_m[pmax-1])

// clean up any division by zero errors
w_d = numtype(w_d)==0 ? w_d : 0

return w_d
end

// one-sided three-point estimate for the derivative adjusted to be shape
// preserving
static function setEndDerivative(delX0, delX1, m0, m1)
variable delX0, delX1, m0, m1

variable derivative

derivative = ( m0*(2*delX0+delX1) - delX0*m1 ) / (delX0+delX1)
if(sign(derivative)!=sign(m1))
derivative=0
elseif( (sign(m0)!=sign(m1)) && (abs(derivative)>3*abs(m0)) )
derivative=3*m0
endif
return derivative
end

threadsafe function InterpolatePCHIP(w_nodesX, w_nodesY, w_nodesM, x)
wave w_nodesX, w_nodesY, w_nodesM
variable x

variable p0, p1
if(x<wavemin(w_nodesX))
p0=0
elseif(x>=wavemax(w_nodesX))
p0=numpnts(w_nodesX)-2
else
findlevel /P/EDGE=1/Q w_nodesX, x
return nan
endif
p0=floor(v_levelX)
endif

p1=p0+1
return InterpolateHermite(x, w_nodesX[p0], w_nodesY[p0], w_nodesX[p1], w_nodesY[p1], w_nodesM[p0], w_nodesM[p1])
end

// https: en.wikipedia.org/wiki/Cubic_Hermite_spline
// calculate cubic function f(x) such that f(x1)=y1, f(x2)=y2, f'(x1)=m1 and f'(x2)=m2
threadsafe function InterpolateHermite(x, x0, y0, x1, y1, m0, m1)
variable x, x0, y0, x1, y1, m0, m1

variable t=(x-x0)/(x1-x0)
m0*=(x1-x0)
m1*=(x1-x0)
return y0*(2*t^3-3*t^2+1) + m0*(t^3-2*t^2+t) + y1*(-2*t^3+3*t^2) + m1*(t^3-t^2)
end```