Obtaining the Cubic Hermite Spline Derivative on an Arbitrary Interval
Tags: gamedev, spline, derivative, math
Published: 2024-05-12 04:50:05

Recently, I had to work with the Hermite spline for the first time in a while. I realized that I had forgotten some details about how it works. I decided to wrap it up here while I was at it.
 
## 1. Finding Hermite Spline 
### 1.1. The 4 constraints

* `P1`: The initial position, `P`
* `m1`: The end point position, `Q`
* `P2`: The slope/velocity at the initial position, `P`
* `m2`: The slope/velocity at the end point, `Q`

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/u91e0hngn906g1yogamz.png)

(The red points represent `P1`, `P2`, while the blue slopes represent `m1`, and `m2`)

Since it is expressed as a cubic function that has 4 coefficients, it can have 4 constraints. The constraints of Hermite spline are the positions and the velocities (=slope) of the start and end points.

### 1.2. Finding the blending functions 

Let's consider the Hermite spline to be obtained, denoted as `P(t)`.

This, a nondescript cubic function, can be represented in a multiplication of 2 matrices. 

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/9lhda3b964nqexewrqko.png)

A multiplication of a <u>1-by-4</u> matrix and a <u>4-by-1</u> matrix yields a <u>1-by-1</u> matrix. The one and only element in the matrix would be the evaluation of the spline given the `t` value.

To be able to evaluate the Hermite spline function, we need to figure out what the `A` matrix consist of. The `A` matrix comprises the 4 coefficients of a cubic polynomial function, the typical `a`, `b`, `c`, and `d`, as we all are familiar with.

The 4 constraints kick in here to determine the coefficients. 

By definition, `P1` represents the value `P(t)` where `t = 0`, and `P2` does, `P(1)`. `m1` and `m2` represents the slope of `P(t)` where `t = 0, 1` respectively.

Substituting each `t` values in the `T` matrix with 0 and 1 accordingly, the constraints can be grouped as below. 

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/jggdoka4jrcjc1esz0dg.png)

To rearrange only the lateral sides of the equation using matrices, We could represent the matrix C (as in *Constraints*) as a multiplication of the coefficient matrix above and the `A` matrix. 

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/5zcgymq03ygn4hbt712u.png)

By rearranging the equation around it, We could represent `A` in the form above. 

> Note: When finding an inverse matrix, https://matrix.reshish.com/inverse.php comes in quite handy.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/cntf0arai05qqc1da1pk.png)

By substituting the newly found `A` matrix into the existing expression `P(t)`, we get the expression above. 

We can take one step further and arrange the equation in the form of `ax + by + cz + d = 0`, but it's not gonna be useful. We will want to have it parameterized by the constraints, As your code will likely provide the constrains, not the coefficients. We will want to plug in them as-is.

The polynomial multiplied by each constraint is called a **basis function** or a **blending function**.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/4lptuh1d53g8yls9gbfq.png)

Then, in conclusion, the spline, and its derivative could be represented in the form below, where `0 <= t <= 1`.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/d5uoeqhxdncungrs9ylv.png)

## 2. On an arbitrary interval `[x_{k}, x_{k+1}]`
The Hermite spline obtained above is valid only in the interval where `t` is in the canonical `[0, 1]` interval. 

* When `t = 0`, the evaluation and the slope would be `P1`, and `m1` .
* When `t = 1`, they would be `P2` and `m2`.

In the ranges before and after 0 and 1, the evaluated values and the slopes will be completely meaningless to the given constraints setting, as it only works in the canonical [0, 1] range. 

However, there are times when we want to move `t` along an arbitrary interval, but the canonical one. For example, what if we want to interpolate the values ​​with the desired position and velocity when `t` is in the interval `[3, 5]`? 

You might think: *'I can just normalize `t` to the length of the interval, to map `[3, 5]` to `[0, 1]`. Simple 1:1 interpolation. It would work, right?'*

Nope, it didn't work.
 
For instance, let's say you want to find a Hermite spline with the constraints setting:
* Starting point: (position: 1, velocity: 1)
    * `P1` = 1, `m1` = 1
* End point: (position: 0.7, velocity: 0)
    * `P2` = 0.7, `m2` = 1

This will yield a spline valid only in the `[0, 1]` range. We could map the parameter `t` from it to the `[3, 5]` range by normalizing it as below.

{%katex%}
t_{normalized} = (t - 3) / (5 - 3)
{%endkatex%}

However, after the normalization, what you're getting is the curve in black color, whereas the actual curve we intended to have is in green color.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/1suj827xbbxeilzz1z5q.png)

(The orange curve is the spline drawn in the normalized interval.)

When `t` is normalized as above, the original orange graph on the left (`[0, 1]` section) gets stretched left and right by the length of the target interval, then moved to the side by the beginning value of the target interval, yielding the black graph.

Stretching the graph left and right like this, **the slope (=velocity) of the graph will be stretched as well**, resulting in a stretched velocity different from the original. As this creates a spline that does not align with the intended constraints, **this is an incorrect way to map the intervals.**

What we want is the <u>green graph</u>, where the slope remains consistent across the range aligning with the initial constraints. To achieve this, we need to **compensate** for the reduced velocity that results from stretching the graph horizontally. To do this, simply **multiply** the velocity **by the length of the interval** to account for the stretch.

The spline curve ` P(x)`, the green one, on an arbitrary section below is as follows. 

{%katex%}
interval = [x_k, x_{k+1}]
{%endkatex%}

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/morjayixuf5mxoin5j10.png)

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/44zd4ljkwno152f0lxbo.png)

`p_k`, `p_{k+1}` are the positions, `m_k`, `m_{k+1}` are the velocity. 

You can see that the velocity constraints are multiplied by the length of the interval for stretch compensation.

This slope factor compensation is also mentioned in the [Wikipedia article](https://en.wikipedia.org/wiki/Cubic_Hermite_spline#Interpolation_on_an_arbitrary_interval 
).

## 3. The Derivative of the Cubic Hermit Spline

The code I was working on required me to provide **the derivative of the cubic Hermite spline**.

The derivative of spline `P(x)` on an arbitrary interval is as follows. 

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/cysstznn4b0ykeriikd7.png)

(Since the constraints are mere scalars, they could be placed out of the differentiation brackets)

However, the blending functions are functions of `t`, not of `x`. Currently, they are differentiated with respect to `x`, so they need to be converted to differentiation with respect to `t`.

In fact, conveniently, the derivatives of each function with respect to `x` and `t` have the following relationship.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/nza1mtcist4wh59t0m7p.png)

Let's take `h_00(t)` as an example. 
1. First, substitute `t` with `x` and differentiate
2. Then, rearrange it for `t`,

as shown below.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/otkd5w62f5b10ge9zw5x.png)

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ea7kcle6imcbag5wiqu2.png)

The numerator of the final result is **identical** to <u>the derivative of `h00(t)` with respect to `t`</u>. Substituting yields the form on the right side above.

This pattern unfolds similarly for the other blending functions, leading to the following final result.

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/rqpax6endpz57wuwx0i2.png)

(As for the slope constraints, since `(x_{k+1} - x_{k})` is already in the numerator, it cancels itself out.)

That is, in fact, the green graph shown above. 

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/lgueu9bqxwadkdix9186.png)

* The black graph is the incorrectly stretched one.
* The green graph is the correctly stretched, and compensated one.

Instead of code, I have parameterized the derivative of a cubic Hermite spline in the Desmos link below so that you can examine it. 

Try changing the value of `t_{slider}` and observe how the slope graph touches and moves along the green curve.

* https://www.desmos.com/calculator/h8taue3l9q

> Note: When calculating derivatives, this website comes in handy: https://www.derivative-calculator.net/ 

## Conclusion

1. Understanding the Hermite Spline Constraints:
    * The cubic Hermite spline is defined by four constraints: the initial and final positions (`P1` and `P2`) and their corresponding velocities (`m1` and `m2`). 
    * Understanding these constraints and their relationship to the cubic polynomial's coefficients helps you draw the blending functions of the cubic Hermite spline. 
    * This step is crucial for mapping the spline's range and performing further differentiation.

2. Mapping from `[0, 1]` to another Arbitrary Range:
    * Mapping the Hermite spline from the canonical `[0, 1]` interval to an arbitrary range `[x, y]` requires compensation for the interval's lateral scaling effect. 
    * The velocity constraints, need to be adjusted to maintain the initially desired spline across different ranges. 
    * This makes sure that the spline accurately interpolates the points and maintains the correct curvature.

3. Deriving the Spline's Derivative:
    * It is necessary to find the derivative with respect to `x` and then represent it as a derivative with respect to `t`.
    * Together with the arbitrary range mapping above, be aware that velocity compensation must also be included. Pay close attention to which terms are canceled out in this process.