#math #geometry #spline #interpolation #number


Pythaogrean hodograph splines

10 releases (4 breaking)

0.5.3 Dec 31, 2022
0.5.2 Dec 28, 2022
0.4.0 Dec 27, 2022
0.3.3 Dec 27, 2022
0.1.0 Dec 10, 2022

#228 in Math

Download history 50/week @ 2022-12-09 28/week @ 2022-12-16 122/week @ 2022-12-23 72/week @ 2022-12-30 2/week @ 2023-01-06 21/week @ 2023-01-13 8/week @ 2023-01-20

108 downloads per month


294 lines


Build Status License Crates.io API

Implementations of Pythagorean hodograph (PH) splines in Rust, based on the work of R. Farouki and this paper by Zbyněk Šír & Bert Jüttler.

What are Pythagorean hodographs?

Pythaogrean hodographs are interpolating polynomial curves that satisfy an analogue of the Pythagorean identity. To illustrate, let $\mathbf{p}(t)$ be a polynomial parametric curve in $\mathbb{R}^2$. If there exist polynomials $u(t), v(t), w(t)$ such that the following holds:

$$\begin{equation} \frac{d\mathbf{p}(t)}{dt} = \begin{bmatrix} u(t)^2 - v(t)^2 \ 2u(t)v(t) \ \end{bmatrix}\cdot w(t) \end{equation}$$

then $\mathbf{p}(t)$ is a Pythaogrean hodograph, because $(u^2-v^2)^2+(2uv)^2=(u^2+v^2)^2,$ a perfect square of a polynomial. It then follows that $|\mathbf{p}'(t)|=(u(t)^2+v(t)^2)\cdot w(t).$ $\mathbf{p}(t)$ is then a planar PH curve.

There are many useful properties of Pythagorean hodographs as a result of this:

  1. The tangent vector, $\frac{\mathbf{p}'(t)}{|\mathbf{p}'(t)|},$ has rational functions as components. This makes computing the tangent fast, although square roots nowadays have dedicated CPU instructions.

  2. The arc-length function $s(t)=\int_0^t|\mathbf{p}'(\tau)|\ d\tau$ is a polynomial function.

  3. As a corollary of #2, computing the arc-length parameterization is trivial using Newton's method.

  4. In three dimensions, Pythagorean hodographs have a naturally associated moving frame called the Euler-Rodrigues frame.

Spatial PH curves

Well-studied readers may have noticed a connection to complex numbers. In particular, if $x(t)=u(t)+iv(t),$ then $p(t)=x(t)^2$ (where $\mathbf{p}(t)$ and $p(t)$ are equivalent via the map from $\mathbb{R}^2\rightarrow\mathbb{C}$).

Complex numbers are isomorphic to the plane, but we need something else if we want to model three dimensions -- quaternions. One has the following presentation: $\mathbf{r}(t)$ is a spatial PH curve if there exists a quaternionic polynomial curve $\mathcal{A}(t)$ such that

$$\begin{equation} \begin{gathered} \begin{aligned} \mathcal{A}(t)&=u(t)+v(t)\cdot\mathbf{i}+p(t)\cdot\mathbf{j}+q(t)\cdot\mathbf{k} \ r'(t)&=-\mathcal{A}(t)\ \mathbf{k}\ \bar{\mathcal{A}}(t)\ &=-2[u(t)p(t)+v(t)q(t)]\cdot\mathbf{i}+2[u(t)v(t)-p(t)q(t)]\cdot\mathbf{j}+2[u(t)^2-v(t)^2-p(t)^2+q(t)^2]\cdot\mathbf{k} \end{aligned} \end{gathered} \end{equation}$$

It then holds that $|\mathbf{r}'(t)|=|\mathcal{A}(t)|^2.$

The Euler-Rodrigues frame

With spatial PH curves, there is a naturally associated moving frame called the Euler-Rodrigues frame. One has the following formula for the tangent, normal, and binormal:

$$\begin{equation} \begin{bmatrix} \mathbf{T}(t) \ \mathbf{N}(t) \ \mathbf{B}(t) \end{bmatrix} = \begin{bmatrix} \mathcal{A}(t)\ \mathbf{i}\ \bar{\mathcal{A}}(t)\ \mathcal{A}(t)\ \mathbf{j}\ \bar{\mathcal{A}}(t) \ \mathcal{A}(t)\ \mathbf{k}\ \bar{\mathcal{A}}(t) \end{bmatrix} \end{equation}$$

or, more succinctly, the rotation matrix $[\mathbf{T\ N\ B}]$ is equivalent to the 3-dimensional rotation associated with the quaternion $\mathcal{A}(t).$


~65K SLoC