Open Treer opened 2 years ago
Hi Treer. Thank you so much for testing it and sharing your results. I'm so sorry for not having enough time to review my code and being able to fix it soon. Since I implemented it as part of the development of my current game I could not spend too much time to make it work perfectly as there are many other components that require my attention. I encourage whoever wants to fork the project and dig into the code to make it work properly. Really, I would love to do it myself, but currently it's physically impossible for me.
Sloan's algorithm looks damn hard (but so fast), so thanks for publishing and diagramming your work - gives some shoulders to stand on.
rough update for people from the future, since I may be switching to another library for unrelated reasons.
Edit: I just noticed the part of the algorithm in question has changed between Sloan's two papers. My explanation below assumed the former paper but this code implements a mix of both systems (e.g. pushing triangles 6 & 7 instead of 0 & 3 is the latter system). The wrong "P" is being checked against when 6 & 7 get popped off the stack. I will have to read Sloan's '92 paper more closely to figure out how the correct P is known - it sounds like [in the 1992 system] you keep the same P until the stack is empty, rather than determining it by the adjacent[OPPOSITE_TRIANGLE_INDEX]
of the most recently popped triangle [the 1987 system].
Conceptually I think the problem is in adjacentTrianglesToProcess
: the stack of triangles needing their Delaunay constraint checked, it's really a stack of edges - with the quadrilateral of each edge checked. Using triangle indexes is a memory optimization but we need to be thinking edges (or quadrilaterals).
Everything works when points are being added, as the quadrilateral is always on the triangle's second edge (the one opposite p[0], given by adjacent[OPPOSITE_TRIANGLE_INDEX]
) and the code knows this. However, when a triangle is added to the stack because it happened to be next to one which just had its edges swapped (i.e. triangles 6 and 7 in the diagram below) then its second edge is often the wrong triangle for it to check against.
e.g. instead of checking triangle 6 against triangle 0. it checks triangle 6 against one of 6's other neighbours - a quadrilateral which hasn't changed. (pretend triangle 6 has other neighbours in that diagram)
There's a comment in the Fortran code on page 49, "PUT EDGES L-A and R-B ON STACK". Below that comment I think Sloan's code is NOT adding triangle 6 and 7 to the stack, but rather adding triangles 0 and 3, i.e. adding the edge between 0-6 and the edge between 3-7, with the edge-swap routine which just rearranged triangles 0 and 3 having been performed in such a way that triangle 0's second edge gives triangle 6, and triangle 3's opposite is triangle 7.
the text on page 38 is similiar (triangles L & R in the PDF would be 3 & 0 in the diagram above, likewise A & B → 7 & 6):
If a swap is necessary, as shown in Fig. 7, the vertex and adjacent arrays for triangles L and R are updated (again with P as the first vertex in the vertex arrays), as are the adjacency arrays for triangles A and C. Provided that there is a triangle opposite P which is adjacent to L, triangle L is placed on the stack. Similarly for triangle R.
Stepping through the triangulation of my test points to watch when it goes astray is how I've come to believe it's testing the correct-triangle-but-wrong-quadrilateral after an edge swap.
Edit: Pull Request https://github.com/QThund/ConstrainedDelaunayTriangulation/pull/3 looks like it was a fix for this, if so then ticket can be closed
I put this code into a project which unit-tests its Delaunay triangulation, and it suggests that the triangulations this code produced are wrong. The set of points it tested are not close together, or inline with each other, or zero on any of their axes, and there are no holes:
Debugging an implementation of Sloan's algorithm looks difficult,
so before attempting that I'd like to be certain it's actually failing and I haven't missed anything - does anybody get the correct triangulation with these test points? Is anybody else checking that their triangulation output conforms to the Delaunay constraint (i.e. no other points in a circumcircle)?(diagram can be inspected at https://www.desmos.com/calculator/j2fipckoxg)
I've noticed the order of the points passed to
DelaunayTriangulation.Triangulate()
changes the triangulation this code produces, so there might be an order where the test passes, but that would be luck.The code listing below gives me the following output:
FWIW, the test-points come from a randomly generated set which (by chance) had p1 fall so close to the circumcircle of p0-p2-p8 that it tripped up a single-point precision circumcircle calculation, however that is not the issue here, as moving p1 doesn't fix anything, there are more triangles wrong than just that pair, and I'd converted
IsPointInsideCircumcircle()
to double precision anyway.