Dealing with floating point precision
Last time I described the straight skeleton algorithm, and write about how precision error can impact calculation. Here, I’ll break down the methods and heuristics I’m using to try to calculate the skeleton for complex shapes.
Calculation windows
Instead of process edge event completely independently, let’s find a “window” of near simultaneous events to calculate each step. I found that without this, we could often get “missed motorcycles” – essentially a motorcycle crash should have happened, but never did, and the vertex proceeds along past it’s crash point. A missed motorcycle like this causes the edges of the wavefront to cross each other. This results in a shape that is unsolvable, and typically the entirely skeleton is ruined.
Near simultaneous events
Imagine that we have two events, A and B, which occur almost at the same calculated time, but B is some tiny amount of time later. A and B have both use some of the same vertices and edges in some way. The naive solution calculates A first, and then scans the resultant shape for the next event that should happen. Commonly, we’ll still find event B and proceed, and everything is fine. However, if the time difference between the events is small enough, the precision error related to processing A can be enough so that when we calculate B the second time, we now find that it happens slightly earlier than A (rather than after it).
Now we have a bit of a dilemma – do we assume that events that happen minutely before the last processed event are subject to precision error and still need to be processed? If so how big should that window big. How do we handle cases such as the 2 new vertices generated by a motorcycle crash event – these are generated on the same point and moving away from each other, so can be calculated as crashing into each other minutely in the past? (those vertices are generated on separate wavefronts, but as we’ll see that doesn’t solve this issue entirely).
This is one way in which the calculation window helps, as described a little below.
Chains of collapses
In the hex grid case, we often get multiple adjacent edges collapsing at the same time and point. Let’s imagine two collapse events, A and B, like this. If we calculate collapse A, we end up generate a vertex at the collapse point. That vertex will happen to be exactly where B is collapsing to (since both edges collapse at the same point). This means we’ve generated an infinitesimally short edge. If the length of an edge is less than or near the error of it’s vertices, that means we can’t accurately calculate the edge normal, which also means we can’t calculate the paths for the vertices nor the collapse point of the edge (meaning we might miss a collapse).
In the hex grid case, this gets even worse when we have a motorcycle crash happening at the same point and time as a the collapses. Here the motorcycle crash can be calculated to happen before the collapses, between the collapses or after the collapses; each of which can cause different problems.
Maintaining the calculation window
Every iteration we generate a list of the next few events that happen on a wavefront. We’ll accept all of the events, up until the gap between successive events exceeds some threshold (which will be very tiny). Now:
- we will process each and every event on this list, regardless of how one event impacts others
- we will process them mostly in order (but with some small reordering covered below)
- while processing the window, new vertices generated are considered to be frozen in place (ie, we don’t even bother generating vertex velocities until all events are processed)
Frequently, events that occur within this window are going to involve the same vertices and edges. So, when one event is processed that removed edges or vertices, we must adjust the other events in the window to take into account these cases. The simplest example of this is when a motorcycle crash is occurring on an edge that is simultaneously collapsing. If we process the collapse first, we must update the motorcycle crash event to now be crashing with the newly created collapse vertex (see the special case for motorcycle-vs-vertex below).
This gets a little more complex with the full range of events that can occur; but this is essential to ensure that every event within the window is actually calculated.
Processing collapse chains
As mentioned before, when adjacent edges are collapsing simultaneously we should process them together and generate a single collapse vertex from the result. We still have a separate event in the calculation window for each collapse, but when processing a collapse we will look for collapses on adjacent edges that occur at any point in the event window.
This is the only time we actually process events out-of-order. In the case where a motorcycle crash involving the same edges is calculated between the collapse events, this would splitting the collapse processing and end up generating a lot of overlapping points. Instead, once we start processing edge collapses, we just make sure we process all adjacent events that occur within the calculation window. Each chain then generates just a single new collapse vertex, rather than one new vertex per collapsed edge.
Special case for motorcycle-vs-vertex
Naively, motorcycle crashes are a vertex vs edge operation. However, as mentioned earlier, when motorcycle intersects another vertex very closely, there can be ambiguity about which side of the vertex it is really on, and therefore which edge to use. Furthermore, this risks generated infinitesimal short edges. So we need to introduce a special case for motorcycle-vs-vertex, wherein we imagine the vertices as small circles (that is, with the radius begin their error). This adds some complexity the motorcycle crash processing, but fortunately it doesn’t impact much else.
Generate velocity from collapse events
The calculations for vertex velocity, edge collapse points and motorcycle crash events are all fairly similar – they all basically boil down to edges moving along their tangents in different ways. However they are subject to precision error in different ways, meaning that we can end up with results that don’t agree with each other.
This can be an issue with vertices attached to very small edges, or edges with very acute angles between them. In these cases our accuracy is lowest and the likelihood of disagreements is highest.
In an ideal world, the collapse points for an edge should be exactly on the path traced by both it’s vertices. So, in order to minimize disagreement, we can use the edge collapse point and time directly as a way to calculate the vertex velocity. This raises a problem, though, in that all vertices are connected to two edges – so which collapse do we use? They should actually both be along the same vertex path, but one might have better precision than the other. I found that just using the soonest event was fine.
It also helps to use this calculated velocity in the calculation of motorcycle crashes. It is possible to calculate the motorcycle crash point directly from vertex points themselves (without using this secondarily derived velocity property). However, in order to ensure that collision point is as close to our calculated vertex path as possible, I found it was more reliable to look for crashes along the calculated vertex velocity. This also inherits the special cases for collinearity mentioned below.
Special cases for collinear edges
As mentioned earlier, the vertex velocity calculations do not hold up under the presence of collinear edges. Furthermore, we end up with very large precision error in near-collinear cases. We can do our best to try to avoid this type of geometry entering the calculations, but we still need some special cases to avoid extreme cases.
There are 2 problems to be wary of:
- edges that make a straight or nearly-straight line, with the vertex somewhere in the middle
- edges meeting in an extremely acute angle, like the vertex of a V squished horizontally
In case (1), there are infinite valid velocities. We need a special case to ensure we get the velocity that is the most intuitive (ie, just moving with the normal of the combined straight edge). In case (2) we approach an asymptote where the velocity is infinite.
To determine “near collinearity” we need some kind of epsilon/threshold value. This and the calculating window size are probably the main places where we use epsilon values. I found the best approach was to filter out case (1) before running the standard velocity calculation math, but to still calculate case (2) with the standard math. This will allow some very high velocities to enter the calculation occasionally. The calculation window actually helps resolve the worst cases of these, though – recall the edges with a collapse event at one end and a motorcycle crash event at the other mentioned in the previous page? These cases can result in extreme acute angles; but since both events will fall within the same calculation window, and since we don’t update velocities during a calculation window, we can side step the issue.
Collapses leading to collinearity
As mentioned above, sometimes non-adjacent collinear edges can become adjacent as a result of edge collapses between them. Adjacent collinear edges created in this fashion are at least as problematic as if they were part of the input data; so we have to do something about them.
Consider the hex grid skeleton from last time:
Here if we look at a part in detail, we can see an example of this
One option here is to just place the vertex at the collapse location as usual, and calculate it’s velocity as the normal of the straight edge that it is on. This adds some extra complexity, however, because as soon as we introduce the possibility of collinear adjacent edges, we need to handle that case in almost all of the calculations. That includes new special cases in both the edge collapse and motorcycle crash detection code.
Another option is to just not include that collapse vertex at all. The collapse vertex would follow the path in the dotted line shown, but – as will become clearer soon – this vertex doesn’t actually have any practical impact on the output result. In fact the output has all of the same properties if it’s there or if it’s not there.
So I elected to simply remove it from the wavefront, which means joining the 2 edges to either side of the collapse into a single edge. We use the same near collinearity detection as mentioned above. This also has a positive side-effect in that very short edges will be detected as collinear to many other edges, which can sometimes mean collapsing down problematic geometry before it causes any real problems.
That sounds great; but as it always seems to be with this algorithm, even this scenario introduces it’s own special cases. You may notice that removing that dotted line in the above diagram means that two of the grey input edges become connected by one shape in the diagram. This is actually significant, and this is the only case in which is can happen. So we can to record this so we can use it later.
Avoid changing vertices after they are constructed
Whenever we create a new vertex and connect it up via edges, we will need to calculate its velocity and resultant path. The velocity related to the collapse points of the connected edges, which is something we’ll need to calculate also. In the edge collapse case, this should be relatively trivial, because a collapse should not change the normal or collapse point of any other edge. But consider the motorcycle crash case, where we introduce edges that are entirely new, and so their collapse points and attached vertex velocities need to be calculated from scratch.
So, for each of the new edges, we’ll calculate its normal and use this to determine all of the properties we need to calculate. To get this normal, we need the position of the vertex at the same point in time. Since we don’t know the velocity of any newly created vertices just yet, we’ll set that point in time to be the time of the crash event itself. Then we can calculate the collapse points and from those collapse points we’ll find the velocity of the new vertices.
This introduces an interesting loop – we use vertex velocities to calculate edge collapses, and edge collapses to calculate vertex velocities. This runs into the compounding error problems mentioned previously. The newly generated vertices are very likely to be involved with future events, which will mean that their velocities will be required for generating tangents for future edges.
It may seem tempting to advance all vertices on the wavefront forward to the time of the event and then find the new positions and recalculate future collapses and motorcycle crashes based on this. In theory, this is valid because the loop created at every subsequent step should still be a valid input loop to the algorithm. However, this will most likely to compound precision error over time and might result in a final shape that is a bit distorted.
Instead I never recalculate the positions or velocities of vertices after they are first calculated. This means that even if a vertex is involved in many events, so long as it is never removed it will still have the same properties as when it was created (even if it was a vertex from the original input loop).
Also note that since we can calculate many events within a calculation window, and don’t update the velocities until the end of the window, we can create situations where both vertices of an edge were generated in the same window by different events. This means the vertices will have different initial times, but we won’t have a velocity to align them in time. In these cases we treat the vertices as if they were created at exactly the same time and assume that the different in calculated normal would be within the error threshold anyway.
Calculating epsilons
In some cases (such as the near collinearity checks) we need an epsilon value to represent the quantity of error we think is in the calculated results. This determines when we can be confident about a calculation and when it falls into the realm of uncertainty. This is mostly used for the collinearity checks and for calculation window. Given the way that I’m using “time” here, time and space have the same units, so the epsilon values in these cases also have the same units.
I’ve ended up just using a fixed value for this. That may not be perfectly ideal. In theory is should at least scale based on the size of the input shape, I guess. But this may be a case where the simple solution is actually the best.
It may be possible to do a little better – for example, we could calculate the amount of error that the underlying calculations (eg, collapse point determination) can deal with before they generate poor results. Or alternatively we could assign an epsilon value per vertex, and calculate this based on the amount of precision we get from the calculations used to generate that vertex. That might help with the compounding error problem; but it could also make things more complicated.
Summary
It took some time to figure out all of these little details and heuristics, I had to try out quite a few different possible approaches to some of the problems before settling these solutions; but these are what have gotten the best results for me. All of the difficult scenarios described in part one can be handled.
Often I found that what seemed like the simplest solution to precision errors – such as “magnetting” together close vertices – actually ended up causing more problems. So there’s often a kind of conflict between ignoring differences within error thresholds (such as collinearity checks, collapse event reordering, motor vs vertex) and trying to maintain as much precision as possible. It is kind of interesting how impactful precision error is on this algorithm in particular. It plays a big part in many computational geometry algorithms; but here we seem to have a combination of calculations converging on the same values while also having a lot of cases (such as a missed motorcycle or missed collapse) that will instantly ruin the output if the precision we need falls under the error thresholds.
Next time
There’s still a little to cover – next time cover calculating skeletons for different types of shapes (such as silhouettes created in art packages) and some applications in 3d.