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 one-dimensional piecewise monotone cubic splines is available in the class splines.PiecewiseMonotoneCubic.

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

We use a few helper functions from helper.py for plotting:

from helper import plot_spline_1d, grid_lines

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

Let’s compare a piecewise monotone spline with a Catmull–Rom spline and a natural spline:

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

def plot_piecewise_monotone(*args, **kwargs):
    s = splines.PiecewiseMonotoneCubic(*args, **kwargs)
    plot_spline_1d(s)
    grid_lines(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.

—Fritsch and Carlson [FC80], section 1

Fritsch and Carlson [FC80] derive 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

For the first step, they suggest using the standard three-point difference, which we have already seen in the tangent vectors of non-uniform Catmull–Rom splines and which is implemented in the class splines.CatmullRom.

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}\).

—Fritsch and Carlson [FC80], section 4

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.

—Fritsch and Carlson [FC80], section 5

In the 2001 edition of the book by de Boor [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 Fritsch and Carlson [FC80], Dougherty et al. [DEH89] provide (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 algorithm is implemented in the class splines.PiecewiseMonotoneCubic.

PCHIP/PCHIM

A different approach for obtaining slopes that ensure monotonicity is described by Fritsch and Butland [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 by Fritsch [Fri82], who 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.

—Fritsch [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 online. 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 that they are in fact also equivalent to the aforementioned equation from Fritsch and Butland [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 by Moler [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

To illustrate the differences between the two approaches mentioned above, let’s plot a few examples. Both methods are piecewise monotone, but their exact shape is slightly different. Decide for yourself which one is more “visually pleasing”!

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', linestyle='--')
    plt.plot(
        plot_times,
        splines.PiecewiseMonotoneCubic(values, times).evaluate(plot_times),
        label='PiecewiseMonotoneCubic', linestyle='-.')
    plt.legend()
    grid_lines(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 entire 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(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

However, in some situations it might be useful to automatically infer the same slope at the beginning and end of the spline. This can be achieved with the cyclic flag.

s = splines.MonotoneCubic([0, 1, 5])

s_cyclic = splines.MonotoneCubic([0, 1, 5], cyclic=True)

plot_spline_1d(s, label='not cyclic')
plot_spline_1d(s_cyclic, label='cyclic')
grid_lines(s.grid)
plt.legend();

The cyclic flag is only allowed if the first and last slope is None:

try:
    splines.MonotoneCubic([0, 1, 5], slopes=[1, None, None], cyclic=True)
except Exception as e:
    print(e)
    assert 'cyclic' in str(e)
else:
    assert False

If "cyclic", the first and last slope must be None

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.

—Fritsch and Carlson [FC80], section 4

Fritsch and Carlson [FC80] recommend using the noncentered three point difference formula, however, they fail to mention what that actually is. Luckily, we can have a look at the code:

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.

—Moler [Mol04], section 3.4

In section 3.6, Moler [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(grid)

plot_pchip([0, 1, 0], [0, 1, 2])
plot_pchip([0, 1, 1], [0, 1, 2], linestyle='--')
grid_lines([0, 1, 2])

plot_pchip([0, 1, 0], [0, 1, 4])
plot_pchip([0, 1, 0], [0, 1, 1.5], linestyle='--')
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([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

The following example plots show different slopes at the beginning and end due to different end conditions.

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])