This page was generated from doc/euclidean/piecewise-monotone.ipynb. Interactive online version: .

Piecewise Monotone Interpolation

When interpolating a sequence of one-dimensional data points, it is sometimes desirable to limit the interpolant between any two adjacent data points to a monotone function. This makes sure that there are no overshoots beyond the given data points. In other words, if the data points are within certain bounds, all interpolated data will also be within those same bounds. It follows that if all data points are non-negative, interpolated data will be non-negative as well. Furthermore, this makes sure that monotone data leads to a monotone interpolant (see also Monotone Interpolation below).

A Python implementation of piecewise monotone one-dimensional cubic splines is available in the splines.PiecewiseMonotoneCubic class.

The SciPy package provides a similar tool with the pchip_interpolate() function and the PchipInterpolator class (see below for more details).

The 3D animation software Blender provides an Auto Clamped property for creating piecewise monotone animation cuves.

Examples

import matplotlib.pyplot as plt
import numpy as np

import splines

helper.py

from helper import plot_spline_1d, grid_lines

values = 0, 3, 3, 7
times = 0, 3, 8, 10, 11

plot_spline_1d(
    splines.PiecewiseMonotoneCubic(values, times, closed=True),
    label='piecewise monotone')
plot_spline_1d(
    splines.CatmullRom(values, times, endconditions='closed'),
    label='Catmull–Rom')
plot_spline_1d(
    splines.Natural(values, times, endconditions='closed'),
    label='natural spline')
plt.legend()
grid_lines(times)

def plot_piecewise_monotone(*args, **kwargs):
    s = splines.PiecewiseMonotoneCubic(*args, **kwargs)
    plot_spline_1d(s)
    grid_lines(x=s.grid)

plot_piecewise_monotone([0, 1, 3, 2, 1])

Providing Slopes

By default, appropriate slopes are calculated automatically. However, those slopes can be overridden if desired. Specifying None falls back to the auto-generated default.

plot_piecewise_monotone([0, 1, 3, 2, 1], slopes=[None, 0, None, -3, -1.5])

Slopes that would lead to non-monotone segments are prohibited:

try:
    plot_piecewise_monotone([0, 1, 3, 2, 1], slopes=[None, 4, None, None, None])
except Exception as e:
    print(e)
    assert 'too steep' in str(e)
else:
    assert False

Slope too steep: 4

Generating and Modifying the Slopes at Segment Boundaries

In this paper we derive necessary and sufficient conditions for a cubic to be monotone in an interval. These conditions are then used to develop an algorithm which constructs a \(\mathscr{C}^1\) monotone piecewise cubic interpolant to monotone data. The curve produced contains no extraneous “bumps” or “wiggles”, which makes it more readily acceptable to scientists and engineers.

—[FC80], section 1: “Introduction”

[FC80] derives necessary and sufficient conditions for a cubic curve segment to be monotone, based on the slopes of the secant lines (i.e. the piecewise linear interpolant) and their endpoint derivatives. Furthermore, they provide a two-step algorithm to generate piecewise monotone cubics:

1. calculate initial tangents (with whatever method)

2. tweak the ones that don’t fulfill the monotonicity conditions

To implement Step 1 we have found the standard three-point difference formula to be satisfactory for \(d_2\), \(d_3\), \(\cdots\), \(d_{n-1}\).

—[FC80], section 4: “Monotone piecewise cubic interpolation algorithm”

This is what de Boor [[dB78], p. 53] calls cubic Bessel interpolation, in which the interior derivatives are set using the standard three point difference formula.

—[FC80], section 5: “Numerical examples”

In the 2001 edition of [dB78], piecewise cubic Bessel interpolation is defined on page 42.

For the following equations, we define the slope of the secant lines as

\begin{equation*} S_i = \frac{x_{i+1} - x_i}{t_{i+1} - t_i}. \end{equation*}

We use \(x_i\) to represent the given data points and and \(t_i\) to represent the corresponding parameter values. The slope at those values is represented by \(\dot{x}_i\).

Note

In the literature, the parameter values are often represented by \(x_i\), so try not to be confused!

Based on [FC80], [DEH89] provides (in equation 4.2) an algorithm for modifying the initial slopes to ensure monotonicity. Adapted to our notation, it looks like this:

\begin{equation*} \dot{x}_i \leftarrow \begin{cases} \min(\max(0, \dot{x}_i), 3 \min(|S_{i-1}|, |S_i|)), & \sigma_i > 0,\\ \max(\min(0, \dot{x}_i), -3 \min(|S_{i-1}|, |S_i|)), & \sigma_i < 0,\\ 0, & \sigma_i=0, \end{cases} \end{equation*}

where \(\sigma_i = \operatorname{sgn}(S_i)\) if \(S_i S_{i-1} > 0\) and \(\sigma_i=0\) otherwise.

This algoritm is implemented in the splines.PiecewiseMonotoneCubic class.

PCHIP/PCHIM

A different approach for obtaining slopes that ensure monotonicity is described in [FB84], equation (5):

\begin{equation*} G(S_1, S_2, h_1, h_2) = \begin{cases} \frac{S_1 S_2}{\alpha S_2 + (1 - \alpha) S_1} \quad & \text{if } S_1 S_2 > 0,\\ 0 & \text{otherwise,} \end{cases} \end{equation*}

where

\begin{equation*} \alpha = \frac{1}{3} \left(1 + \frac{h_2}{h_1 + h_2}\right) = \frac{h_1 + 2 h_2}{3 (h_1 + h_2)}. \end{equation*}

The function \(G\) can be used to calculate the slopes at segment boundaries, given the slopes \(S_i\) of the neighboring secant lines and the neighboring parameter intervals \(h_i = t_{i+1} - t_i\).

Let’s define this using SymPy for later reference:

import sympy as sp

h1, h2 = sp.symbols('h1:3')
S1, S2 = sp.symbols('S1:3')

alpha = (h1 + 2 * h2) / (3 * (h1 + h2))
G1 = (S1 * S2) / (alpha * S2 + (1 - alpha) * S1)

This has been implemented in a Fortran package described in [Fri82], which has coined the acronym PCHIP, originally meaning Piecewise Cubic Hermite Interpolation Package.

It features software to produce a monotone and “visually pleasing” interpolant to monotone data.

—[Fri82]

The package contains many Fortran subroutines, but the one that’s relevant here is PCHIM, which is short for Piecewise Cubic Hermite Interpolation to Monotone data.

The source code (including some later modifications) is available at https://people.sc.fsu.edu/~jburkardt/f77_src/pchip/pchip.html. This is the code snippet responsible for calculating the slopes:

C
C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
C
   45    CONTINUE
         HSUMT3 = HSUM+HSUM+HSUM
         W1 = (HSUM + H1)/HSUMT3
         W2 = (HSUM + H2)/HSUMT3
         DMAX = MAX( ABS(DEL1), ABS(DEL2) )
         DMIN = MIN( ABS(DEL1), ABS(DEL2) )
         DRAT1 = DEL1/DMAX
         DRAT2 = DEL2/DMAX
         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)

This looks different from the function \(G\) defined above, but if we transform the Fortran code into math …

HSUM = h1 + h2

W1 = (HSUM + h1) / (3 * HSUM)
W2 = (HSUM + h2) / (3 * HSUM)

… and use separate expressions depending on which of the neighboring secant slopes is larger …

G2 = S1 / (W1 * S1 / S2 + W2 * S2 / S2)
G3 = S2 / (W1 * S1 / S1 + W2 * S2 / S1)

… we see that the two cases are mathematically equivalent …

assert sp.simplify(G2 - G3) == 0

… and they are in fact also equivalent to the aforementioned equation from [FB84]:

assert sp.simplify(G1 - G2) == 0

Presumably, the Fortran code uses the larger one of the pair of secant slopes in the denominator in order to reduce numerical errors if one of the slopes is very close to zero.

Yet another variation of this theme is shown in [Mol04], section 3.4, which defines the slope \(d_k\) as a weighted harmonic mean of the two neighboring secant slopes:

\begin{equation*} \frac{w_1 + w_2}{d_k} = \frac{w_1}{\delta_{k-1}} + \frac{w_2}{\delta_k}, \end{equation*}

with \(w_1 = 2 h_k + h_{k-1}\) and \(w_2 = h_k + 2 h_{k-1}\). Using the notation from above, \(d_k = \dot{x}_k\) and \(\delta_k = S_k\).

Again, when defining this using SymPy …

w1 = 2 * h2 + h1
w2 = h2 + 2 * h1

G4 = (w1 + w2) / (w1 / S1 + w2 / S2)

… we can see that it is actually equivalent to the previous equations:

assert sp.simplify(G4 - G1) == 0

The PCHIM algorithm, which is nowadays known by the less self-explanatory name PCHIP, is available in the SciPy package in form of the pchip_interpolate() function and the PchipInterpolator class.

from scipy.interpolate import PchipInterpolator

More Examples

def compare_pchip(values, times):
    plot_times = np.linspace(times[0], times[-1], 100)
    plt.plot(
        plot_times,
        PchipInterpolator(times, values)(plot_times),
        label='PCHIP')
    plt.plot(
        plot_times,
        splines.PiecewiseMonotoneCubic(values, times).evaluate(plot_times),
        label='PiecewiseMonotoneCubic')
    plt.legend()
    grid_lines(x=times)

compare_pchip([0, 0, 1.5, 4, 4], [-1, 0, 1, 8, 9])

compare_pchip([0, 0, 1.5, 4, 4], [-1, 0, 6, 8, 9])

There is even a slight difference in the uniform case:

compare_pchip([0, 0, 3.3, 4, 4], [-1, 0, 1, 2, 3])

compare_pchip([0, 0, 0.7, 4, 4], [-1, 0, 1, 2, 3])

For differences at the beginning and the end of the curve, see the section about end conditions.

Monotone Interpolation

When using the aforementioned piecewise monotone algorithms with monotone data, the whole interpolant will be monotone.

The class splines.MonotoneCubic works very much the same as splines.PiecewiseMonotoneCubic, except that it only allows monotone data values.

Since the resulting interpolation function is monotone, it can be inverted. Given a function value, the method .get_time() can be used to find the associated parameter value.

s = splines.MonotoneCubic([0, 2, 2, 6, 6], grid=[0, 2, 3, 6, 8])

probes = 1, 3, 5

fig, ax = plt.subplots()
plot_spline_1d(s)
ax.scatter(s.get_time(probes), probes)
grid_lines(x=s.grid)

If the solution is not unique (i.e. on plateaus), the return value is None:

assert s.get_time(2) is None

Closed curves are obviously not possible:

try:
    splines.MonotoneCubic([0, 2, 2, 6, 6], closed=True)
except Exception as e:
    print(e)
    assert 'closed' in str(e)
else:
    assert False

The "closed" argument is not allowed

End Conditions

The usual end conditions don’t necessarily lead to a monotone interpolant, therefore we need to come up with custom end conditions that preserve monotonicity.

For the end derivatives, the noncentered three point difference formula may be used, although it is sometimes necessary to modify \(d_1\) and/or \(d_n\) if the signs are not appropriate. In these cases we have obtained better results setting \(d_1\) or \(d_n\) equal to zero, rather than equal to the slope of the secant line.

—[FC80], section 4: “Monotone piecewise cubic interpolation algorithm”

[FC80] recommends using the noncentered three point difference formula, however, it fails to mention what that actually is. Luckily, we can have a look at the code at https://people.sc.fsu.edu/~jburkardt/f77_src/pchip/pchip.html:

C
C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      HSUM = H1 + H2
      W1 = (H1 + HSUM)/HSUM
      W2 = -H1/HSUM
      D(1,1) = W1*DEL1 + W2*DEL2
      IF ( PCHST(D(1,1),DEL1) .LE. ZERO)  THEN
         D(1,1) = ZERO
      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         DMAX = THREE*DEL1
         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
      ENDIF

The function PCHST is a simple sign test:

PCHST = SIGN(ONE,ARG1) * SIGN(ONE,ARG2)
IF ((ARG1.EQ.ZERO) .OR. (ARG2.EQ.ZERO))  PCHST = ZERO

This implementation seems to be used by “modern” PCHIP/PCHIM implementations as well.

This defines the pchip slopes at interior breakpoints, but the slopes \(d_1\) and \(d_n\) at either end of the data interval are determined by a slightly different, one-sided analysis. The details are in pchiptx.m.

—[Mol04], section 3.4

Section 3.6 of [Mol04] shows the implementation of pchiptx.m:

function d = pchipend(h1,h2,del1,del2)
%  Noncentered, shape-preserving, three-point formula.
    d = ((2*h1+h2)*del1 - h1*del2)/(h1+h2);
    if sign(d) ~= sign(del1)
        d = 0;
    elseif (sign(del1)~=sign(del2))&(abs(d)>abs(3*del1))
        d = 3*del1;
    end

Apparently, this is the same as the above Fortran implementation.

The class scipy.interpolate.PchipInterpolator uses the same implementation (ported to Python).

This implementation ensures monotonicity, but it might seem a bit strange that for calculating the first slope, the second slope is not directly taken into account.

Another awkward property is that for calculating the inner slopes, only the immediately neighboring secant slopes and time intervals are considered, while for calculating the initial and final slopes, both the neighboring segment and the one next to it are considered. This makes the curve less locally controlled at the ends compared to the middle.

def plot_pchip(values, grid, **kwargs):
    pchip = PchipInterpolator(grid, values)
    times = np.linspace(grid[0], grid[-1], 100)
    plt.plot(times, pchip(times), **kwargs)
    plt.scatter(grid, pchip(grid))
    grid_lines(x=grid)

plot_pchip([0, 1, 0], [0, 1, 2])
plot_pchip([0, 1, 1], [0, 1, 2])
grid_lines([0, 1, 2])

plot_pchip([0, 1, 0], [0, 1, 4])
plot_pchip([0, 1, 0], [0, 1, 1.5])
grid_lines([0, 1, 1.5, 4])

In both of the above examples, the very left slope depends on properties of the very right segment.

The slope at \(t = 1\) is clearly zero in both cases and apart from that fact, the shape of the curve at \(t > 1\) should, arguably, not have any influence on the slope at \(t = 0\).

To provide an alternative to this behavior, the class splines.PiecewiseMonotoneCubic uses end conditions that depend on the slope at \(t = 1\), but not explicitly on the shape of the curve at \(t > 1\):

plot_piecewise_monotone([0, 1, 0], grid=[0, 1, 1.5])
plot_piecewise_monotone([0, 1, 0], grid=[0, 1, 4])
grid_lines(x=[0, 1, 1.5, 4])

The initial and final slopes of splines.PiecewiseMonotoneCubic are implemented like this:

def monotone_end_condition(inner_slope, secant_slope):
    if secant_slope < 0:
        return -monotone_end_condition(-inner_slope, -secant_slope)
    assert 0 <= inner_slope <= 3 * secant_slope
    if inner_slope <= secant_slope:
        return 3 * secant_slope - 2 * inner_slope
    else:
        return (3 * secant_slope - inner_slope) / 2

Even More Examples

compare_pchip([1, 2, 1], [1, 3.5, 5])

compare_pchip([1, 2, 3.5, 4, 3], [1, 1.5, 4, 5, 6])

compare_pchip([1, 2, 1.9, 1], [1, 3, 4, 6])