Finding Catmull-Rom Spline and Line intersection. Part 2: Mathematical approach

As far as I can tell there’s no open-source solution out there for this problem to date. Maybe it’s trivial for those with significant mathetamical knowledge, but for those who are just starting to get familiar with the mathematical foundations of splines this problem may just seem too intimidating. As a result they might resort to a much worse solution they’re comfortable with. So did I at first.

So let this article be a guide for a solution. I will explain the math and also provide an implementation in Java using the libGDX library.

That said if you’re a least bit interested I’d recommend you trying to solve this for youself. All you need is some basic understanding of the Catmull-Rom Spline equation and the Line equation. As a hint I recommend reading this post. Though it actually contains a flaw it will put you on the right track.

catmull-rom formula

To determine a point on the Catmull-Rom Spline we need the 4 sorrounding control points (P0, P1, P2, P3 , predefined), a knot parameter (alpha, predefined) and a value in range of 0-1 (t) that will tell which point of the Spline we’ll get.

The formula can be written as:

q(t) = alpha * ((2*P1) + (-P0+P2) * t + (2*P0 – 5*P1 + 4*P2 – P3) * t^2 + (-P0 + 3*P1 – 3*P2 + P3) * t^3))

Which (if we disregard the predefined values) is a simple Cubic equation:

a*t^3 + b*t^2 + c*t + d = 0

Where

a = alpha*(-P0 + 3*P1 – 3*P2 + P3)

b = aplha*(2*P0 – 5*P1 + 4*P2 – P3)

c = alpha * (-P0+P2)

d = aplha * (2*P1)

Our points (P0-3) consist of two variables: the x and y coordinates. They’re independent so we can handle them seperately. Let’s say:

x = a1*t^3 + b1*t^2 + c1*t + d1

y = a2*t^3 + b2*t^2 + c2*t + d2

Where e.g. ‘a1’ is ‘a’ with the x coordinates substituted for each P point and ‘a2’ is ‘a’ with the y coordinates substituted for each P point and so on with b, c and d.

line FORMULA

Ax+By+C=0

Not much to see here, move along…

combined FORMULA

Substituting the values of x and y from the separated Spline formula into the Line formula we’ll get:

A*(a1*t^3 + b1*t^2 + c1*t + d1) + B*(a2*t^3 + b2*t^2 + c2*t + d2) + C = 0

Rearranged:

(A*a1 + B*a2)*t^3 + (A*b1 + B*b2)*t^2 + (A*c1 + B*c2)*t + (A*d1 + B*d2 + C) = 0

This is again a Cubic equation, combined from the Line’s and the Spline’s equation. If you solve this you can get up to 3 real roots for t.

You remember that t is ranged from 0 to 1, right? Choose those that fit, substitue them into the Spline formula we first talked about and you have the intersection point(s).

If none of them fits they don’t intersect. If more than one fits there’s more than one intersection (3 possible with a knot as seen on the uniform spline here). In case you only need one intersection point you should examine all possibilities and make a conscious decision (e.g. take first/last point on the Line/Spline) otherwise the outcome might be random.

multiple spans

The above solution is simplified to a Spline with 4 control points (1 span). What about Splines that has more spans?

You have to do the same calculation for all spans using their 4 sorrunding control points. Keep in mind that there may be intersections in every span.

implementation

The code (alongside with other Catmull-Rom utilities) is available on GitHub. The Span class takes care of the math, it takes values needed from the Line class’ interface. Finally, the Spline class’ intersect(Line) method will gather the results. There’s a desktop demo where you can try it for yourself.

Finding Catmull-Rom Spline and Line intersection. Part 1: Brute-force

This was an issue I tackled for a while and ended up with a pretty neat optimized solution – at least so I thought at the time. As it turns out it’s completely useless in light of the matchematical approach I later implemented. Still it was a pretty nice try and makes for a good story on how a programmer vs. a mathematician would approach a problem. At the end it seems like I’m both. Lucky me.

Anyway this line of thinking (pun not intended) may actually be useful in other cases where you’re left without a pure mathematical approach.

Brute-force 101

So the problem is given. Assuming you don’t want to take a deep-dive in equations you’re left with one choice. That one’s actually fairly straightforward.

Let’s take all points of your Line and cross-reference it with all points of your Spline and see if they have anything in common. Or at least within a certain margin of error. Being a smart developer you chose a library that has already a built-in way to iterate these objects. Your code might look something like this:

for (float l = 0f; l <= 1f; l+= precision ) {
    linePoint = line.valueAt(l);
    for (float s = 0f; s <= 1f; s += precision) {
        splinePoint = spline.valueAt(s);
        if (linePoint.dst(splinePoint) < delta) return splinePoint;
    }
}

Precision and delta should be chosen carefully. You’d probably want to chose delta first to set a maximum inaccuracy. Then you need to set the precision so that you can actually find the data you’re looking for. This may be much harder than it first seems especially considering how loose some libraries play with precision.

And there you have it: a very low-effort solution that can actually work just fine for you.

complexity cut

Well, it did not quite work for me. If you need to do this more than once, let alone continuously in real-time the above solution will really start to push it.

Time to get creative!

… but first of all you can’t really better anything if you can’t measure it. Assuming your  Line has L points and yor Spline has S the above solution yields a complexity of O(L*S). Simplifying it to a common precision (N) it’s O(N^2). Now even if you don’t know much about the order of algorithms you should know that anything with an exponent is bad.

Given that we’re brute-forcing (i.e. trying all possibilities) we can’t make our algorithm’s complexity independent from the object’s complexity. But maybe we can ignore some possibilities along the way.

X-flat spline

Here comes the idea: let’s construct our Spline in a way that it is continuous on one axis. My setting is vertical so I chose X.

This will enable us to simplify things a lot. Take a random point. Take the first and last point of your Spline. Just by looking at the X coordinates you can tell if that point has no chance at all of being on the Spline (i.e. it’s not between the Spline’s ends). If it does you still have check though but remember, we’re all about eliminating possibilities here. If your iteration allows it you can actually do the same thing with the Spline’s spans (parts between adjacent control points). At this point if your Spline has P control points you essentially saved (P-2)/(P-1)th of the work (that is, the most) not even counting points that are outside of the Spline’s reach ((e.g. if you have 6 control points which result in 5 spans and you identify the span you need to look at you saved looking at the other 4)).

By this point you might have realied that it’s possible to do a binary search on the Spline (or only  on 1 span) for the X coordinate. Once you found the Spline’s point that has the same X coordinate as your random point you check the Y as well for the final verdict. And that actually only takes O(log(S)) time which is a very significant improvement. Your code might look something like this:

var spline, delta;

binaryFindByX (possibleIntersection, firstX, lastX) {
    middleX = firstX + (lastX - firstX) / 2;
    middlePoint = spline.pointAt(middleX);

    if (Math.abs(possibleIntersection.x - middlePoint.x) <= delta)
        return possibleIntersection.dst(middlePoint) <= delta ? possibleIntersection : null;
    else if (possibleIntersection.x < middlePoint.x)
        return binaryFindByX(possibleIntersection, firstX, middleX);
    else
        return binaryFindByX(possibleIntersection, middleX, lastX);
 }

binaryFindByX(pointOfLine, spline.getFirstX, spline.getLastX)

That said you still have to iterate through the Line’s points but do note that I was trying to shorten the processing of the Spline. In my case a Spline was longer and by definition more complex than the Line so it made more sense, but you can very much apply this logic the other way around.

The result is O(L*log(S)) complexity. That is a great improvement even by itself but in practice when we apply these modifications on the more complex object this really makes a huge difference in performance.

implementation

The XFlatSpline implementation (alongside with other Catmull-Rom utilities) is available on GitHub. The ‘calculateIntersection’ function will yield all intersections and uses ‘getPointOnSpline’ which implements the binary search based on the X coordinate. There’s a desktop demo where you can try it for yourself.

libGDX – Iterating a Catmull-Rom Spline by fixed units

If you want a well defined curve Catmull-Rom is your go-to spline. It’s being used very often in graphics so I was glad to see that libGDX (a cross-platform game development library) has built-in support for it. But if you want to do some more advanced calculations you can easily run into problems.

out of the box iteration with libgdx

The most basic operation you’ll ever do with a spline is iterating over its points. Let’s look at rendering as an example:

path = new CatmullRomSpline<>(controlPoints, false); //false means non-continuous; bad design right there

for (float t = 0f; t <= 1f; t+= precision ) {
 Vector2 point = new Vector2();
 path.valueAt(point, t);
 shapeRenderer.circle(point.x, point.y, height);
}

In the above code we’re rendering small circles on the path of the spline (a dotted line) using libGDX’s CatmullRomSpline API.

Questions is: how do we choose the precision? Essentially it will determine the number of points rendered, therefore the smoothness of the spline but it also highly affects performace. It should obviously correlate with how long or pointed the spline is.

As length is fairly trivial and supported by the API we set up a small demo using something along the lines of:

int precisionOverLength = 5;
float length = path.approxLength(numberOfSamples);
precision = precisionOverLength / length;

Where precision over length determines your relative precision. Assuming length is in pixels: 1 would mean 1 point per pixel (very precise), 5 is 1 point per 5 pixels (less precise).

Let’s see the result:

Capture3

The big dots – in case you’re not familiar with Catmull-Rom splines – are the control points allowing you to shape the curve as it will always pass over those points. The actual curve is calculated from these points using a mathematical equation. Another term I must define is a span. A span is the section of the spline between two control points – this one has seven.

Did you notice the density in spans 3 and 4? Despite our efforts to achieve universal precision in these spans the precision is way higher than expected. The reason is that the number of points we’ll calculate on the spline is divided equally between spans. A shorter span therefore will have a higher density of points, a higher precision.

the upside

In many cases (and especially when rendering) this is helpful as sharp changes in direction in a small area could cause our spline to be fragmented. Having lots of points in these areas will make them smoother.

the downside

It will take as much time to compute a small area as to compute a large one even if that area is very simple and would not warrant a higher precision. Having no control over this may lead to performance issues.

You can’t guarantee a fixed or minimal density of points during iteration. This lack of guaranteed precision is most painful if you’re working with approximations and margins of error.

Imagine you’re looking for a specific point (e.g. one third) or an intersection of your spline. It has one very long and many short spans. If you want to be precise on the long span you’ll end up spending the unnecessarily big efforts on the shorter spans. And there’s still no guarantee that you’ll find what you’re looking for within the margin of error.

iteration by fixed units

Luckily for us there’s another method provided by the API where you can iterate only points of a specific span. This way we can move our precision calculation to the level of spans and bypass the equal point division between spans.

This will require quite some additional logic to keep track of where we are in the spline and what precision we should use in each span.

the result

SinglePointSplineIterator iterator = new SinglePointSplineIterator(path);
SplinePoint splinePoint;

while ((splinePoint = iterator.getNext()) != null) {
    shapeRenderer.circle(splinePoint.point.x, splinePoint.point.y, height);
}

Capture4

The code (alongside with other Catmull-Rom utilities) is available on GitHub. There’s an XFlatSpline class using the iterator (here’s why) and a desktop demo where you can try it for yourself. As a QoL improvement there’re separate single– and double point iterators available as in many cases you’ll need the current and previous points as well during iteration.

Finally a side-by-side comparison:

capture