# How to calculate a spline for points in general position?

I want to find a curve passing through (or near) $n$ points in the plane. The catch is that the curve need not be a function. That is, a vertical line might pass through the curve in more than one place. I want to treat the $x$ and $y$ coordinates identically.

So, I have a set of points, $(x_1, y_1), ..., (x_n, y_n)$ where it is NOT necessarily the case that $x_1 < x_2 < ... < x_n$ as is the usual condition for a spline, but I'd like to use a spline-like approximation.

I'm sure this must be a problem that has been thoroughly addressed in the literature, but I can't figure out what to search for. If someone could at least tell me some keywords to search for, I would appreciate it.

Edit: After doing some more thinking and hunting around, I believe the key is going to be to use separate splines to represent the $x$ and $y$ coordinates in terms of some other parameter, $t$, i.e., as a set of parametric equations. So, for example, if I had four points, I could represent $x$ and $y$ as cubic polynomials:

$x = at^3 + bt^2 + ct + d$ and $y = et^3 + ft^2 + gt + h$. Since $t$ is an arbitrary parameter, I could let it take on integer values, e.g., $t = 0,1,2,3$ in order to make solving the system easier.

The problem is that $t$ isn't really a natural parameter. I would prefer something like arc length so that I could calculate and presumably minimize the bending energy. Somehow, I will have to incorporate that into my calculation at a later stage.

Darrell

Your thinking is correct -- you need to use a parametric curve (either a spline or a single polynomial) that gives $x$ and $y$ as functions of some parameter $t$. To compute the curve, you must first assign a parameter value $t_i$ to each point $P_i$. As you say, these parameter values are entirely fabricated -- they are not an intrinsic part of the original interpolation problem. Furthermore, different choices for these parameter values will produce significantly different curves. The usual approach is to choose them so that $t_{i+1}-t_i$ is equal to the chord-length between the points $P_i$ and $P_{i+1}$.

Using arclength seems like the natural approach, as you say, but this doesn't work because you don't know the arclengths until you've constructed the curve -- a chicken-and-egg problem.

A good reference is this web site.

See Raph Levien's phd thesis, which pretty much does what you want, and produces an arc length parameterization:

http://www.levien.com/phd/phd.html

Here's a super-simplified overview of Levien's approach:

1. Write a polynomial; this will be the curvature function. Give it enough degrees of freedom to satisfy all constraints (e.g. for G2 curves this is 2).
2. Integrate that with respect to s to get the angle function $th(s)$.
3. The x/y equations are now:

a. $x = \int_{0}^{_s}$ sin(th(s))

b. $y = \int_{0}^{_s}$ cos(th(s))

4. For each desired segment in the spline. . .

1. Set the coefficients of the curvature polynomial to zero.
2. To evaluate a point on the curve at position $s$. . .
1. Find the start and end points of the curve by integrating x and y numerically.
1. Calculate the translation, rotation and scale needed for the curve endpoints to coincide with the two input points that define the input segment.
2. Evaluate the curve at position $s$ using numerical integration. Apply transformation to result.
3. Using Newton-Raphson, solve for the coefficients that satisfy the constraints you want (e.g. tangent continuity, G2 continuity, etc). Note that for curvature continuity you will have to reverse the scaling you did earlier to get world-space curvature.

You do end up solving a large system of nonlinear equations, but the results are awesome. It produces much better curvature behavior than you get with simple polynomials.

Note that Levien has code on his site that performs these steps for the polynomial he came up with in his thesis.