Bezier Fitting

The goal is to fit a sequence of cubic Bezier curve over a sequence of ordered points. We use a classic Bezier fitting technique developed by Philip Schneider. A Bezier curve $Q$ of order $n$ is defined by its control points $v_i$: $Q(u) = \sum_{i=0}^n v_i B_i^n(u)$, where $B_i^n(u)$ are the Bernstein polynomials of degree $n$, $B_i^n(u) = {n \choose i} u^i (1-u)^{n-i}$. In our case, we choose $n=3$, i.e. the cubic Bezier curves with four control points.

We need to minimize the distance from a Bezier curve to a set $S$ of $n$ points. Formally, we define an objective function $F = \sum_{j=1}^n ||p_i - Q(u_i)||^2$, where $u_i$ is the parameter of the closest point on cubic Bezier curve $Q$ to the point $p_i$.

Let $t_1$ and $t_2$ be the tangents at the start and end of the curve. Then write the second and the third control points in terms of the tangents, $v_1 = \alpha_1 t_1 + v_0$ and $v_2 = \alpha_2 t_2 + v_3$ where $\alpha_1$ and $\alpha_2$ are the distances to the nearest end points in the tangent direction. Therefore, the function $F$ can be minimized by finding the right parameters $\alpha_1$ and $\alpha_2$. In particular, we set the derivatives of $F$ with respect to $\alpha_1$ and $\alpha_2$, and set them to zero:

$\frac{\partial F}{\partial \alpha_1} = \sum_{i=1}^n 2\left(p_i - Q(u_i)\right) \cdot t_1 B_1^3(u_i) = 0$

$\frac{\partial F}{\partial \alpha_2} = \sum_{i=1}^n 2\left(p_i - Q(u_i)\right) \cdot t_2 B_2^3(u_i) = 0$

This linear system of equations can be solved via matrix inversion. Let $A_{i,j} = t_j B_{j}^3(u_1)$. Then,

$\left(\sum_{i=1}^n A_{i,1}^2\right) \alpha_1 + \left(\sum_{i=1}^n A_{i,1} \cdot A_{i, 2}\right) \alpha_2 = \sum_{i=1}^n \left(p_i - \left( v_0 B_0^3(u_i) + v_0 B_1^3(u_1) + v_3 B_2^3(u_i) + v_3 B_3^3(u_i) \right)\right) \cdot A_{i,1}$

$\left(\sum_{i=1}^n A_{i,2}^2\right) \alpha_2 + \left(\sum_{i=1}^n A_{i,1} \cdot A_{i, 2}\right) \alpha_1 = \sum_{i=1}^n \left(p_i - \left( v_0 B_0^3(u_i) + v_0 B_1^3(u_1) + v_3 B_2^3(u_i) + v_3 B_3^3(u_i) \right)\right) \cdot A_{i,2}$

Condensing the equations into a matrix,

$\begin{bmatrix} c_{1,1} & c_{1,2} \\ c_{2,1} & c_{2,2} \end{bmatrix} \begin{bmatrix}\alpha_1 \\ \alpha_2 \end{bmatrix} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$

Then solving these two equations allow us to obtain a good approximation.

There is still the problem of finding the parameter of the closest point for each point $p_i$. We can use Newton-Raphson iteration to get a better $u$: $u \leftarrow u - \frac{f(u)}{f'(u)}$ where $f = (p_i - Q(u)) \cdot Q'(u)$ and $f' = Q'(u) \cdot Q'(u) + (p_i - Q(u)) \cdot Q''(u)$.