Assume we are given the set $S$ of $n$ points on the real plane and want to draw a parametrized cubic curve (actually a segment of Bézier spline) with fixed startpoint in such a way, that the resulting curve is closest to $S$. $S$ may be assumed to contain at least three different points. By closest I mean such a curve, that the sum of distances between it and every point of $S$ is minimal (however I think the sum of squares of distances will be ok too). Drawing a cubic curve with fixed startpoint amounts to specifying 6 coefficients in the parametric presentation and the problem is to find way of expressing these coefficients through the given coordinates of points in $S$. Is there any solution to the problem?
$\begingroup$
$\endgroup$
4
-
$\begingroup$ You probably want to work with the sum of squares of distances. Then standard least-squares methods should work. But I haven't actually done the algebra. $\endgroup$– Michael LugoCommented Apr 7, 2011 at 15:58
-
$\begingroup$ LMC will work, but it minimizes slightly different thing, namely the sum of squares of distances between points of S and some specific knots (each for every point of S) on the curve. For programming purposes this is enough, of course. $\endgroup$– isnmrCommented Apr 7, 2011 at 16:18
-
2$\begingroup$ Honestly, I do not see how the LSM works here because, unlike the distance to the line, the (signed) distance to the curve is not linear in the parameters (actually, it has no nice explicit formula whatsoever). I wonder if there is some convexity in the problem that would ensure at least that the local minimizer is unique. Then the gradient descend method has a chance to work :) $\endgroup$– fedjaCommented Apr 7, 2011 at 17:57
-
$\begingroup$ fedja, that's a good point. Like I said, I didn't actually do the algebra. $\endgroup$– Michael LugoCommented Apr 8, 2011 at 13:58
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
The method I'll present below, is described (briefly) in the paper "A Revisit to Least Squares Orthogonal Distance Fitting of Parametric Curves and Surfaces" by Liu and Wang. Other methods for solving this problem are also presented in that paper and its references.
We wish to minimize the following sum-of-squares error of the cubic Bézier curve $C(u)$:
\begin{equation} \min \sum_{k=0}^{n-1}(C(u_k) - S_k)^2 = \min \sum_{k=0}^{n-1}(\sum_{i=0}^{3}P_i B_i(u_k) - S_k)^2 \end{equation}
where $P_i$ are the Bézier control points to be computed (the output), $B_i(u)$ are the Bézier basis functions, $S_k$ are the input points and $u_k$ are the parameters of the curve corresponding to the input points.
If $u_k$ are given as input, then this can be solved as a linear least square problem. The equation can be interpreted as finding the minimal-norm solution (in the least square sense) of an over-determined set of linear equations of the form:
\begin{equation} \left[ \begin{array}{cccc} B_0(u_0) & B_1(u_0) & B_2(u_0) & B_3(u_0) \\ & \vdots & & \\ & \vdots & & \\ B_0(u_k) & B_1(u_k) & B_2(u_k) & B_3(u_k) \\ & \vdots & & \\ & \vdots & & \\ B_0(u_{n-1}) & B_1(u_{n-1}) & B_2(u_{n-1}) & B_3(u_{n-1}) \end{array} \right] % \left( \begin{array}{c} P_0 \\ P_1 \\ P_2 \\ P_3 \end{array} \right) = \left( \begin{array}{c} S_0 \\ \vdots \\ \vdots \\ S_k \\ \vdots \\ \vdots \\ S_{n-1} \end{array} \right). \end{equation}
Denoting $B$ as the left-hand matrix, $P$ as the column vector $(P_0,P_1,P_2,P_3)^T$, and $S$ as the right hand column vector (of size $n$), the standard least squares solution to this over-determined linear system is the solution of: \begin{equation} (B^T B) P = B^T S. \end{equation}
However, in the context of the question, $u_k$ are not given and are part of the unknowns. They are the parameters where the minimal distance should be attained. Thus, the problem is non-linear (as @fedja noted in the comments).
A standard heuristic of fitting a curve is to assign (normalized) chord-length values to the $u_k$. While a practical solution, it only gives a coarse approximation to the real minimum attaining curve. However, it can be used as an initial approximation for an iterative process that converges to the minimum. In each iteration $j$, one alternates between solving the linear least square problem to attain the current Bézier curve $C_j(u)$, and then projecting the points $S_k$ onto $C_j$ to attain improved $u_k$ parameters. By executing these two steps iteratively, improved parameters are obtained.
A note on projecting a point onto a Bézier curve (sometimes called point-inversion, see for example Chapter 6.1 of The NURBS Book). Given a point $S_k$ and a curve $C(u)$, the parameter $u_k$ can be attained by finding the roots of the equation $f(u) = C'(u) \cdot (C(u) - S_k)$. In the case of a cubic Bézier curve $f(u)$ is a degree-5 polynomial and numerical methods such as Newton–Raphson (or any other polynomial solving method) can be used to find the root $u_k$ corresponding to the closest projected point.
The following figure demonstrates using this method. The blue crosses are the input points, the orange curve is the initial guess (achieved with chord-length parameter assignment) and the green curve is the result after ten iterations (the start point is fixed as requested in the original question, which means $P_0$ is set and therefore some small modifications to the formulations above were needed).
An interesting thing to notice is that for some configurations there is no minimal solution. The following figure shows a "corner" configuration where except for the first point, all other points are on a line. The curves in the figure correspond to results after 100, 200 and 300 iterations.
It can be seen that, since there is no constraint on the length of the curve, as the curve extends more and more to the left, it will achieve better and better approximation.
