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
        if(v_flag==1) // not found         
            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

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More