# Piecewise cubic Hermite interpolating polynomial (PCHIP)

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
static function SBL_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=SBL_setDerivativesPCHIP(w_nodesX, w_nodesY)
if(waveexists(xwave))
w_interp=SBL_InterpolatePCHIP(w_nodesX, w_nodesY, w_d, xwave[p])
else
w_interp=SBL_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 SBL_setDerivativesPCHIP(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]=SBL_EndDerivativePCHIP(w_delx[0], w_delx[1], w_m[0], w_m[1])
w_d[pmax+1]=SBL_EndDerivativePCHIP(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 SBL_EndDerivativePCHIP(delX0, delX1, m0, m1)
variable delX0, delX1, m0, m1

variable derivative

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

threadsafe static function SBL_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 SBL_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 static function SBL_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

Earlier today, WaveMetrics support was contacted about this post. Here is what we got. I'm posting his message here at his request.

----

I'd like to inform you of an error that I have found in the following code found on your page.

https://www.wavemetrics.com/code-snippet/piecewise-cubic-hermite-interpolating-polynomial-pchip

I have compared this code to the Octave / Matlab function that does the same and found it to be different. The Octave results seem more correct. The error is in the setEndDerivative function.

if(sign(derivative)!=sign(m1))

should actually be

if(sign(derivative)!=sign(m0))

Here is the Octave code for comparison. See line 186 and 244

http://octave.org/doxygen/4.0/d7/d35/pchim_8f_source.html

I don't know if this error has been fixed in your code already. I just thought I'd let you know so that you can check it.

Regards,

Peter Bone

Thanks, Adam, and thanks also to Peter Bone for pointing out the error.

I've edited the code snippet above. I'll post a revised version of the baseline spline package that uses this code some time soon. Since the code in that package is used to generate subjectively visually appealing baselines based on the user dragging nodes, the only important ramification is that reloading saved node positions may generate a slightly different baseline. My apologies for that.

Here are some comparisons of the splines generated by my code compared with those in figures from Fritsch and Butland (1984) A method for constructing local monotone piecewise cubic interpolants, SIAM Journal on Scientific and Statistical Computing 5: 300-304. For the set of nodes used in these figures the code change has no noticeable effect.

Forum

Support

Gallery