interpolation

In this article, we'll see how we can use cubic Bézier curves to create a smooth line that goes through a predefined set of points. If you don't know what Bézier curves are, you might want to check this post:

Bézier Curve Understand the mathematics of Bézier curves omaraflak.com

Cubic Bézier Curves

The goal is to fit [begin-latex-inline]n+1[end-latex-inline] given points [begin-latex-inline]{P_0, \ldots, P_n}[end-latex-inline]. In order to fit these points, we'll use one cubic Bézier curve (4 control points) between each consecutive pair of points.

cubic bezier curve

In this figure, [begin-latex-inline]G_0[end-latex-inline], [begin-latex-inline]G_1[end-latex-inline], and [begin-latex-inline]G_2[end-latex-inline] are three different cubic Bézier curves that start and end at [begin-latex-inline](P_0, P_1)[end-latex-inline], [begin-latex-inline](P_1, P_2)[end-latex-inline], and [begin-latex-inline](P_2, P_3)[end-latex-inline] respectively. Since any Bézier curve always starts and ends at the first and last control points, we are left to find 2 control points for each curve.

The general equation of the cubic Bézier curve is the following:

[begin-latex]\begin{align*} Z(t) &= \sum_{i=0}^3 B_i^3(t) \cdot K_i \\ &= \sum_{i=0}^3 {3 \choose i} t^i (1-t)^{3-i} \cdot K_i \\ & = (1 - t)^3 \cdot K_0 + 3t(1-t)^2 \cdot K_1 + 3t^2(1-t) \cdot K_2 + t^3 \cdot K_3 \end{align*}[end-latex]

Where [begin-latex-inline]K_i[end-latex-inline] are the 4 control points. In our case, [begin-latex-inline]K_0[end-latex-inline] and [begin-latex-inline]K_3[end-latex-inline] will be two consecutive points that we want to fit (e.g. [begin-latex-inline](P_0, P_1)[end-latex-inline], [begin-latex-inline](P_1, P_2)[end-latex-inline], etc.), and [begin-latex-inline]K_1[end-latex-inline] and [begin-latex-inline]K_2[end-latex-inline] are the 2 remaining control points we have to find. With [begin-latex-inline]t \in [0, 1][end-latex-inline].

Problem Setup

Given that we have [begin-latex-inline]n+1[end-latex-inline] points to fit, we will use a cubic Bézier curve to fit each consecutive pair of points. We denote [begin-latex-inline]\Gamma_i[end-latex-inline] the Bézier curve that fits [begin-latex-inline]P_i[end-latex-inline] to [begin-latex-inline]P_{i+1}[end-latex-inline]:

[begin-latex]\Gamma_i(t) = (1 - t)^3 \cdot P_i + 3t(1-t)^2 \cdot A_i + 3t^2(1-t) \cdot B_i + t^3 \cdot P_{i+1}, \quad i=0,\ldots,n-1[end-latex]

Where [begin-latex-inline]A_i[end-latex-inline] and [begin-latex-inline]B_i[end-latex-inline] are left to find. Notice that there are [begin-latex-inline]n[end-latex-inline] curves.

If we want the final curve to be smooth, we need to ensure that the transition between [begin-latex-inline]\Gamma_i[end-latex-inline] and [begin-latex-inline]\Gamma_{i+1}[end-latex-inline] is smooth around [begin-latex-inline]P_{i+1}[end-latex-inline]. In other words, that the curvature of [begin-latex-inline]\Gamma_i[end-latex-inline] matches the curvature of [begin-latex-inline]\Gamma_{i+1}[end-latex-inline] around [begin-latex-inline]P_{i+1}[end-latex-inline]. Mathematically, this means respecting the following conditions:

[begin-latex]\begin{align*} &\Gamma'_{i-1}(t=1) = \Gamma'_i(t=0), \quad &i=1,\ldots,n-1 \\[2ex] &\Gamma''_{i-1}(t=1) = \Gamma''_i(t=0), \quad &i=1,\ldots,n-1 \\ \end{align*}[end-latex]

We need to find all [begin-latex-inline]A_i[end-latex-inline] and [begin-latex-inline]B_i[end-latex-inline]. Since we have one pair of them in each Bézier curve, and since we have [begin-latex-inline]n[end-latex-inline] curves, we need to find [begin-latex-inline]2n[end-latex-inline] variables. However, here we have [begin-latex-inline]2(n-1)[end-latex-inline] equations. We are missing 2 equations to solve the system. Therefore, we impose the following (arbitrary) boundary conditions:

[begin-latex]\Gamma''_0(t=0) = 0 \\[2ex] \Gamma''_{n-1}(t=1) = 0[end-latex]

Write The System

Before solving the system, we need to calculate the first and second derivatives of [begin-latex-inline]\Gamma_i[end-latex-inline] and write the system down.

[begin-latex]\begin{align*} &\Gamma'_i(t) = 3 \left[-(1-t)^2 \cdot P_i + (1-3t)(1-t) \cdot A_i + t(2-3t) \cdot B_i + t^2 \cdot P_{i+1} \right] \\[2ex] &\Gamma''_i(t) = 6 \left[ (1-t) \cdot P_i + (3t-2) \cdot A_i + (1-3t) \cdot B_i + t \cdot P_{i+1} \right] \end{align*}[end-latex]

Injecting Equations

The boundary condition for the first derivative becomes:

[begin-latex]\begin{align*} & \Gamma'_{i-1}(t=1) = \Gamma'_i(t=0), \quad i=1,\ldots,n-1 \\ \iff & A_i + B_{i-1} = 2P_i \end{align*}[end-latex]

The boundary condition for the second derivative becomes:

[begin-latex]\begin{align*} & \Gamma''_{i-1}(t=1) = \Gamma''_i(t=0), \quad i=1,\ldots,n-1 \\ \iff & A_{i-1} + 2A_i = 2B_{i-1} + B_i \end{align*}[end-latex]

The first arbitrary boundary condition becomes:

[begin-latex]\begin{align*} & \Gamma''_0(t=0) = 0 \\ \iff & P_0 - 2A_0 + B_0 = 0 \end{align*}[end-latex]

The second arbitrary boundary condition becomes:

[begin-latex]\begin{align*} & \Gamma''_{n-1}(t=1) = 0 \\ \iff & A_{n-1} - 2B_{n-1} + P_n = 0 \end{align*}[end-latex]

Solve The System

To sum up, we have the following [begin-latex-inline]2n[end-latex-inline] equations where we need to find [begin-latex-inline]A_i[end-latex-inline] and [begin-latex-inline]B_i[end-latex-inline].

[begin-latex]\left\{ \begin{align} & A_i + B_{i-1} = 2P_i, \quad &i=1,\ldots,n-1 \quad \quad \\ & A_{i-1} + 2A_i = 2B_{i-1} + B_i, \quad &i=1,\ldots,n-1 \quad \quad \\ & P_0 - 2A_0 + B_0 = 0 \quad \quad \\ & A_{n-1} - 2B_{n-1} + P_n = 0 \quad \quad \end{align} \right.[end-latex]

To solve the system, we'll eliminate all the [begin-latex-inline]B_i[end-latex-inline] by injecting [begin-latex-inline](1)[end-latex-inline] into [begin-latex-inline](2), (3), (4)[end-latex-inline]. To start off:

[begin-latex]\begin{align*} & A_i + B_{i-1} = 2P_i, \quad &i=1,\ldots,n-1 \\ \iff & B_i = 2P_{i+1} - A_{i+1}, \quad &i=0,\ldots,n-2 \end{align*}[end-latex]

Injecting into [begin-latex-inline](2)[end-latex-inline]:

[begin-latex]\begin{align*} & A_{i-1} + 2A_i = 2B_{i-1} + B_i, \quad i=1,\ldots,n-1 \\[3ex] \iff & \left\{ \begin{align*} & A_{i-1} + 2A_i = 2(2P_i - A_i) + 2P_{i+1} - A_{i+1}, & \quad i&=1,\ldots,n-2 \\ & A_{n-2} + 2A_{n-1} = 2B_{n-2} + B_{n-1}, & \quad i&=n-1 \end{align*} \right. \\[3ex] \iff & \left\{ \begin{align*} & A_{i-1} + 4A_i + A_{i+1} = 4P_i + 2P_{i+1}, & \quad i&=1,\ldots,n-2 \\ & A_{n-2} + 2A_{n-1} = 2B_{n-2} + B_{n-1}, & \quad i&=n-1 \end{align*} \right. \end{align*}[end-latex]

Injecting into [begin-latex-inline](3)[end-latex-inline]:

[begin-latex]\begin{align*} & P_0 - 2A_0 + B_0 = 0 \\ \iff & P_0 - 2A_0 + 2P_1 - A_1 = 0 \\ \iff & 2A_0 + A_1 = P_0 + 2P_1 \end{align*}[end-latex]

Injecting into [begin-latex-inline](4)[end-latex-inline]:

[begin-latex]\begin{align*} & A_{n-1} - 2B_{n-1} + P_n = 0 \\ \iff & A_{n-1} - 2(A_{n-2} + 2A_{n-1} - 2B_{n-2}) + P_n = 0 \\ \iff & A_{n-1} - 2(A_{n-2} + 2A_{n-1} - 2(2P_{n-1} - A_{n-1})) + P_n = 0 \\ \iff & 2A_{n-2} + 7A_{n-1} = 8P_{n-1} + P_n \end{align*}[end-latex]

Alright! To sump up, after getting rid of [begin-latex-inline]B_i[end-latex-inline], we now have:

[begin-latex]\left\{ \begin{align*} & 2A_0 + A_1 = P_0 + 2P_1 \\ & A_{i-1} + 4A_i + A_{i+1} = 4P_i + 2P_{i+1}, \quad i=1,\ldots,n-2 \\ & 2A_{n-2} + 7A_{n-1} = 8P_{n-1} + P_n \end{align*} \right.[end-latex]

We can write this system as a matrix multiplication and solve it!

[begin-latex]\begin{align*} & \left\{ \begin{align*} & 2A_0 + A_1 = P_0 + 2P_1 \\ & A_{i-1} + 4A_i + A_{i+1} = 4P_i + 2P_{i+1}, \quad i=1,\ldots,n-2 \\ & 2A_{n-2} + 7A_{n-1} = 8P_{n-1} + P_n \end{align*} \right. \\ \\ \iff & \left\{ \begin{align*} & 2A_0 + A_1 = P_0 + 2P_1 \\ & A_0 + 4A_1 + A_2 = 4P_1 + 2P_2 \\ & A_1 + 4A_2 + A_3 = 4P_2 + 2P_3 \\ & \cdots \\ & A_{n-3} + 4A_{n-2} + A_{n-1} = 4P_{n-2} + 2P_{n-1} \\ & 2A_{n-2} + 7A_{n-1} = 8P_{n-1} + P_n \end{align*} \right. \\ \\ \iff & \begin{bmatrix} 2 & 1 & 0 & 0 & 0 & \cdots & 0 \\ 1 & 4 & 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 4 & 1 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 1 & 4 & 1 & 0 \\ 0 & \cdots & 0 & 0 & 1 & 4 & 1 \\ 0 & \cdots & 0 & 0 & 0 & 2 & 7 \end{bmatrix} \begin{bmatrix} A_0 \\ A_1 \\ A_2 \\ \vdots \\ A_{n-3} \\ A_{n-2} \\ A_{n-1} \end{bmatrix} = \begin{bmatrix} P_0 + 2P_1 \\ 4P_1 + 2P_2 \\ 4P_2 + 2P_3 \\ \vdots \\ 4P_{n-3} + 2P_{n-2} \\ 4P_{n-2} + 2P_{n-1} \\ 8P_{n-1} + P_n \end{bmatrix} \end{align*}[end-latex]

As you can see, the first matrix is mostly zeros except for the main 3 diagonals. This kind of matrix is called, reasonably enough, a tridiagonal matrix. Algorithms exist to solve this type of systems efficiently, such as Thomas Algorithm which runs linearly in time.

We don't really care about optimizing here. We'll use the built-in function of Numpy in Python to solve the system.

However, we are still missing the [begin-latex-inline]B_i[end-latex-inline] points. To find those, we use equation [begin-latex-inline](1)[end-latex-inline] which works out all the [begin-latex-inline]B_i[end-latex-inline] up until [begin-latex-inline]B_{n-2}[end-latex-inline] and then equation [begin-latex-inline](4)[end-latex-inline] which gives us the last term, [begin-latex-inline]B_{n-1}[end-latex-inline].

[begin-latex]\left\{ \begin{align*} & B_i = 2P_{i+1} - A_{i+1}, & \quad i&=0,\ldots,n-2 \\ & B_{n-1} = \dfrac{A_{n-1} + P_n}{2}, & \quad i&=n-1 \end{align*} \right.[end-latex]

We are done at last! Let's see how we can program this using Python.

Python Code

import numpy as np
import matplotlib.pyplot as plt

# find the A & B points
def get_bezier_coef(points):
    # since the formula assumes n+1 points, then n must be:
    n = len(points) - 1

    # build coefficents matrix
    C = 4 * np.identity(n)
    np.fill_diagonal(C[1:], 1)
    np.fill_diagonal(C[:, 1:], 1)
    C[0, 0] = 2
    C[n - 1, n - 1] = 7
    C[n - 1, n - 2] = 2

    # build points vector
    P = [4 * points[i] + 2 * points[i + 1] for i in range(n)]
    P[0] = points[0] + 2 * points[1]
    P[n - 1] = 8 * points[n - 1] + points[n]

    # solve system, find a & b
    A = np.linalg.solve(C, P)
    B = [0] * n
    for i in range(n - 1):
        B[i] = 2 * points[i + 1] - A[i + 1]
    B[n - 1] = (A[n - 1] + points[n]) / 2

    return A, B

# returns the general cubic Bezier formula given 4 control points
def get_bezier_cubic(a, b, c, d):
    return lambda t: np.power(1 - t, 3) * a + 3 * np.power(1 - t, 2) * t * b + 3 * (1 - t) * np.power(t, 2) * c + np.power(t, 3) * d

# returns one cubic Bezier curve for each consecutive pair of points
def get_bezier_interpolation(points):
    A, B = get_bezier_coef(points)
    return [
        get_bezier_cubic(points[i], A[i], B[i], points[i + 1])
        for i in range(len(points) - 1)
    ]

# evalute each cubic Bezier curve on the range [0, 1] sliced in `total` points
def evaluate_bezier_interpolation(points, total):
    curves = get_bezier_interpolation(points)
    return np.array([fun(t) for fun in curves for t in np.linspace(0, 1, total)])


# generate 5 (or any number that you want) random points that we want to fit (or set them youreself)
points = np.random.rand(5, 2)

# fit the points with Bezier interpolation
# use 50 points between each consecutive points to draw the curve
path = evaluate_bezier_interpolation(points, 50)

# extract x & y coordinates of points
x, y = points[:,0], points[:,1]
px, py = path[:,0], path[:,1]

# plot
plt.figure(figsize=(11, 8))
plt.plot(px, py, 'b-')
plt.plot(x, y, 'ro')
plt.show()

Here's a demo running in the browser. Drag the red dots with the mouse to move them around!