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.

JPA – the uniqueness dilemma & solution

The dilemma

JPA does not support non-primary (often referred to as business) keys. In simple words if you want your entity to have an Id (as you normally do) and an additional composite key – you’re on your own. There are numerous articles and discussions on the subject worth checking out (yep, each word is a different link).

Let’s see if we can find a solution…

Solution #1 – Drop the Id and use a composite key

To have an Id is standard for many reasons. Even if you think you don’t need the Id (yes you do) you’re being forced to deviate from the standard just because of the lack of a better solution – no way.

Solution #2 – Keep the Id as your primary key and define unique constraints

It will be hard to read (as you’ll have to annotate the entity class instead of the fields) but it doesn’t sound that bad. Until.. welcome to EDD, exception driven development. – Not as good as TDD. The JPA has no intelligent way of letting you know whether the object already exists (often referred to as mergeOrCreate) instead it just occasionally throws a ConstraintViolationException in your face.

Meaning your implementation will look something like:

  • Try persisting
  • Got an exception? It’s already in the DB, go and fetch it..
  • No exception? Good.

Don’t you dare calling this a solution user453265436 on StackOverflow!

Solution #3 – Keep the Id, keep the constraints and make your own mergeOrCreate

The idea is to check if the object with the given unique fields already exists before trying to insert it into the DB. Luckily for us the good people of the internet offer some pseudo and actual implementations. Now all you have to do is write hashCode and equals functions for all your entities to be able to compare with tarnsient objects. And write equal criteria as well to be able to compare with persisted ones.

Do I seriously have to code this?! I fill sick already. Best part of my day.

Oh and make sure you can test is somehow. Don’t forget to keep hashCode and equals in sync. Also never misuse hashCode as key. Keep in mind what happens with your hash code if you update a collection.

Best part of my week. And we’re nowhere near multithreading yet – oh man, we’re gonna have such a great time!

Do yourself a favour and do it in a genral way as much as posibble. You might want to think about an Interface for your entities to provide their uniqueness criteria and a common EntityFactory or superclass. For some more detailed insight I recommend reading this post in particular from the above ones.

The solution I came up with is available on GitHub.

Below are some additional issues worth addressing. In some cases they’re part of the implementation, in others I’ll explain why they’re not.

Entity relationships

Pssst.. hey! Listen, don’t even bother just put everything in one entity.

If somehow by the grace of the ORM god you didn’t face the issue of uniquness yet, chances are you’re about to when normalizing your schema. Think about building up your entity in-memory that has a ManyToOne connection (Person and City in my example). When persisting it (besides checking the uniqueness of the entity itself) you need to check whether any connected entities already exist.

Solution

You’ll have to implement DAOs for every entity with ManyToOne connection. You’ll need to query that One for each of your Many.

Thought I was exaggerating? Best week ever!

On a sidenote here’s one reason why I prefer Hibernate’s Session API over JPA’s EntityManager API: it updates your object to the persisted state (i.e. assings Id when persisted). It makes synchronizing objects easier.

Those who say the both APIs are the same are dirty liars. They also tend to prefer EM. Don’t listen to them!

PARALLELISM

Aka. just when you thought it shouldn’t get any worse. I don’t even say “it can’t get worse” anymore.

We have to stop here and clarify how the Persistence Context cache works before we make more harm than good. Despite the many articles and examples that indicate otherwise JPA’s EntityManager and Hibernate’s Session offer an extended Persistence Context which means the cache is not bound to transaction but lives across transactions made from the same EntityManager/Session.

This is where Session’s getCurrentSession() comes into play. Another huge advantage of the Session API over EntityManager is that this method allows the threadsafe use of the API.

All you need to do is set ‘hibernate.current_session_context_class’ in your persistence config to ‘thread’. If you forget about it though you’ll get a nice exception during runtime which is, well, not so very nice…

Aaaaand that’s all for the good news. As for the bad news…

As mentioned before when using constraints you will get an ConstraintViolationException when trying to insert an entry that already exists in the DB. This may very well happen when multithreading even if you check whether the objects exists before inserting via mergeOrCreate because another thread may have created it in the meantime.

This might sound unlikely but when working with big amount of data it will happen – and this is exactly the case when you’re relying on your cache the most.

Why is this a problem? The Session / EntityManager object will be invalidated upon an expection and therefore the cache is lost. Meaning parallel execution will likely make your application slower.

Solution

You can avoid the Session invalidation by using nested transactions, however netiher of the mentioned APIs support it. Some others do though (e.g. Spring). If single threading is not an option (e.g. you’re developing a distributed app) you can still use locks.

Batch processing

In some cases this is something that could give you the performance boost you need. It’s not that simple here.

Persisting bigger chunks of data in one Transaction will leave you with probably even more ConstraintViolationExceptions – invalidating the whole Transaction. So instead of persisting multiple entities at once you’ll lose even those that we already persisted within the Transaction and have to start over.

In a single threaded environment that is only the case if the same entries are present within the same chunk. We can actually eliminate this possibility so this might be a possible improvement.

THE final solution?

It’s simple. We kill the batman.

We need support for business keys across application and DB layers. We need a persistence provider that can handle if a business key already exists in the persistence context or the DB.

It’s that simple.

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