# Monotonic Cubic Spline interpolation (with some Rust)

Monotonic Cubic Spline interpolation (MCSI) is a popular and useful method which fits a smooth, continuous function through discrete data. MCSI has several applications in the field of computer vision and trajectory fitting. MCSI further guarantees monotonicity of the smoothed approximation, something which a cubic spline approximation alone cannot.
In this post I’ll show how to implement the method developed by F. N. Fritsch and R. E. Carlson^{1} in the Rust^{2} programming language.

## Rust

*Why* Rust? Definitely this is a type of solution so simple that it can be implemented in practically any programming language we can think of. However, I do find that the best way to get acquainted with a new language and its concepts is precisely to try to implement a simple and well-know solution. Although this post does not intend to be an introduction to the Rust language, some of the fundamentals will be presented as we go along.

Idiomatic Rust Object-Oriented Programming (OOP) has several characteristics which differ significantly from “traditional” OOP languages.
Rust achieves data and behaviour encapsulation by means of defining data structure blueprints (called `struct`

) and then defining their behaviour though a concrete implementation (through `impl`

). As an example, a simple “class” `Foo`

would consist of:

```
struct Foo {
}
impl Foo {
fn new() -> Foo {
return Foo {};
}
fn method(&mut self) {}
fn static_method() {}
}
pub fn main() {
let mut f = Foo::new();
f.method();
Foo::static_method();
}
```

The “constructor” is defined typically as `new()`

, but any “static” method which returns an initialised `struct`

*can* be a constructor and “object” methods include the passing of the `self`

instance not unlike languages such as Python. The `&mut self`

refers to the control or exclusive access to `self`

and it is not directly related to `mut`

mutability control. These concepts touch on Rust’s borrowing and ownership model which, unfortunately, are *way* beyond the scope of this blog post. A nice introduction is provided by the “*Rust programming book*” available here.
Our implementation aims at building a MCSI class `MonotonicCubicSpline`

by splitting the algorithm into the slope calculation at *construction* time, a *Hermite interpolation* function and a *partial application function generator*. This will follow the general structure

```
pub struct MonotonicCubicSpline {
m_x: Vec<f64>,
m_y: Vec<f64>,
m_m: Vec<f64>
}
impl MonotonicCubicSpline {
pub fn new(x : &Vec<f64>, y : &Vec<f64>) -> MonotonicCubicSpline {
// ...
}
pub fn hermite(point: f64, x : (f64, f64), y: (f64, f64), m: (f64, f64)) -> f64 {
// ...
}
pub fn interpolate(&mut self, point : f64) -> f64 {
// ...
}
fn partial(x: Vec<f64>, y: Vec<f64>) -> impl Fn(f64) -> f64 {
// ...
}
}
```

`Vec`

is a vector, a typed growable collection available in Rust’s standard library with documentation available here.

## Monotonic Cubic Splines

MCSI hinges on the concept of cubic Hermite interpolators. The Hermite interpolation for the unit interval for a generic interval $(x_k,x_{k+1})$ is

$$ p(x)=p_k h_{00}(t)+ h_{10}(t)(x_{k+1}-x_k)m_k + \\ h_{01}(t)p_{k+1} + h_{11}(t)(x_{k+1}-x_{k})m_{k+1}. $$

The $h_{\star}$ functions are usually called the *Hermite basis functions* in the literature and here we will use the factorised forms of:

$$ \begin{aligned} h_{00}(t) &= (1+2t)(1-t)^2 \\ h_{10}(t) &= t(1-t)^2 \\ h_{01}(t) &= t^2 (3-2t) \\ h_{11}(t) &= t^2 (t-1). \end{aligned} $$

This can be rewritten as

$$ \begin{aligned} p(x) = (p_k(1 + 2t) + \Delta x_k m_k t)(1-t)(1-t) + \\ (p_{k+1} (3 -2t) + \Delta x_k m_{k+1} (t-1))t^2 \end{aligned} $$

where

$$ \begin{aligned} \Delta x_k &= x_{k+1} - x_k \\ t &= \frac{x-x_k}{h}. \end{aligned} $$

This associated Rust method is the above mentioned “static” `MonotonicCubicSpline::hermite()`

:

```
pub fn hermite(point: f64, x : (f64, f64), y: (f64, f64), m: (f64, f64)) -> f64 {
let h = x.1 - x.0;
let t = (point - x.0) / h;
return (y.0 (1.0 + 2.0 t) + h m.0 t) (1.0 - t) (1.0 - t)
+ (y.1 (3.0 - 2.0 t) + h m.1 (t - 1.0)) t t;
}
```

where the tuples correspond to $x \to (x_k, x_{k+1})$, $t \to (y_k, y_{k+1})$ and $m \to (m_k, m_{k+1})$

For a series of data points $(x_k, y_k)$ with $k=1,\dots,n$ we then calculate the slopes of the secant lines between consecutive points, that is:

$$ \Delta_k = \frac{\Delta y_{k}}{\Delta x_k},\qquad \text{for}\ k=1,\dots,n-1 $$

with $Delta y_k = y_{k+1}-y_k$ and $\Delta x_k$ as defined previously.

Since the data is represented by the vectors `x : Vec<f64>`

and `y : Vec<f64>`

we implement this in the “constructor”:

```
let mut secants = vec![0.0 ; n - 1];
let mut slopes = vec![0.0 ; n];
for i in 0..(n-1) {
let dx = x[i + 1] - x[i];
let dy = y[i + 1] - y[i];
secants[i] = dy / dx;
}
```

The next step is to average the secants in order to get the tangents, such that

$$ m_k = \frac{\Delta_{k-1}+\Delta_k}{2},\qquad \text{for}\ k=2,\dots,n-1. $$

This is achieved by the code:

```
slopes[0] = secants[0];
for i in 1..(n-1) {
slopes[i] = (secants[i - 1] + secants[i]) * 0.5;
}
slopes[n - 1] = secants[n - 2];
```

By definition, we want to ensure monotonicity of the interpolated points, but to guarantee this we must avoid the interpolation spline to go too far from a certain radius of the control points. If we define $\alpha_k$ and $\beta_k$ as

$$ \begin{aligned} \alpha_k &= \frac{m_k}{\Delta_k} \\ \beta_k &= \frac{m_{k+1}}{\Delta_k}, \end{aligned} $$

to ensure the monotonicity of the interpolation we can impose the following constraint on the above quantities:

$$ \phi(\alpha, \beta) = \alpha - \frac{(2\alpha+\beta-3)^2}{3(\alpha+\beta-2)}\geq 0, $$

that is

$$ \alpha + 2\beta - 3 \leq 0, \text{or}\ 2\alpha+\beta-3 \leq 0 $$

Typically the vector $(\alpha_k, \beta_k)$ is restricted to a circle of radius 3, that is

$$ \alpha^2_l + \beta_k^2>9, $$

and then setting

$$ m_{k+1} = t\beta_k\Delta_k, $$

where

$$ \begin{aligned} h &= \sqrt{\alpha^2_k + \beta^2_k} \\ t &= \frac{3}{h}. \end{aligned} $$

One of the ways in which Rust implements polymorphism is through method dispatch. The `f64`

primitive provides a shorthand for the quantity $\sqrt{\alpha^2_k + \beta^2_k}$ as $\alpha.\text{hypot}(\beta)$. The relevant Rust code will then be:

```
for i in 0..(n-1) {
if secants[i] == 0.0 {
slopes[i] = 0.0;
slopes[i + 1] = 0.0;
} else {
let alpha = slopes[i] / secants[i];
let beta = slopes[i + 1] / secants[i];
let h = alpha.hypot(beta);
if h > 3.0 {
let t = 3.0 / h;
slopes[i] = t * alpha * secants[i];
slopes[i + 1] = t * beta * secants[i];
}
}
}
```

We are now able to define a “smooth function” generator using MCSI. We generate a smooth function $g(.)$ given a set of $(x_k, y_k)$ points, such that

$$ f(x_k, y_k, p) \to g(p). $$

## Partial application

Before anything, it is important to recall the difference between partial application and *currying*, since the two are (incorrectly) used interchangeably quite often.
*Function currying* allows to factor functions with multiple arguments into a chain of single-argument functions, that is

$$ f(x, y, z) = h(x)(y)(z) $$

The concept is prevalent in functional programming, since its initial formalisation^{3}. Partial application, however, generally aims at using an existing function conditioned on some argument as a basis to build functions with a reduced arity. In this case this would be useful since ultimately we want to create a smooth, continuous function based on the control points $(x_k, y_k)$. The partial application implementation is done in Rust as

```
pub fn partial(x: Vec<f64>, y: Vec<f64>) -> impl Fn(f64) -> f64 {
move |p| {
let mut spline = MonotonicCubicSpline::new(&x, &y);
spline.interpolate(p)
}
}
```

An example of how to generate a concrete smoothed continuous function from a set of control points can be:

```
let x = vec![0.0, 2.0, 3.0, 10.0];
let y = vec![1.0, 4.0, 8.0, 10.5];
let g = partial(x, y);
// calculate an interpolated point
let point = g(0.39);
```

The full code can be found here.