Delaunay2D fails when trying to triangulate certain set of points

[Kemikal1]

Godot version

3.5.1

System information

Windows10, Ryzen 5600

Issue description

Self expalnatory, calling Delaunay2D on a set of points fails to return a triangulation. This triangulation can exist in matlab, so most likely godot bug.

Steps to reproduce

Make an array of the points ill show below, calling the Geometry.Delaunay2D(list) returns an empty int array.

Minimal reproduction project

List<Vector2> testList = new List<Vector2>();
testList.Add(new Vector2(0,0));
testList.Add(new Vector2(-13.3f,0));
testList.Add(new Vector2(-6.5f,-0.15f));
testList.Add(new Vector2(-6.8f,-0.15f));
GD.Print(Geometry.TriangulateDelaunay2d(testList.ToArray()).Count());

Use this or the equivalent in GDScript and you get 0 printed

[AThousandShips]

Also happens in 4.x

[Kemikal1]

I can do that triangulation with Matlab tho, doesnt this mean its a godot implementation problem?

[AThousandShips]

Then that is a bug likely, sorry didn't know you had verified it was correct

[Zireael07]

@Kemikal1 The fact you can triangulate this in another Delaunay implementation is relevant and should be up there in the OP

[Kemikal1]

Making the triangle sufficiently large partially solves it, but it now becomes order dependent, will investigate some further

Maybe it has some angle calculations in its implementation idk for sure just spitballin

[AThousandShips]

That was a false trail, that's why I deleted the comment (which I did an hour or so ago)

[Kemikal1]

This isx how it looks in matlab

[AThousandShips]

Note that this doesn't happen if the y-values are positive for this specific case, indicating this is a strange edge case

Be aware that in Godot y values go the other way, so the plot is upside down, and also auto scaled meaning it ignores the significantly stretched shape, reason I'm pointing this out is that this plot doesn't give a proper idea of how the shape actually looks

[Kemikal1]

What do you mean the plot is upside down? Also my usage of the function is for generating meshes in 3D and for that you need point projections and i cant avoid some coordinates being 0 there.

[AThousandShips]

The Y axis in Godot is positive down, not negative down, it's not necessarily relevant but it does affect things like triangle ordering

What I meant with positive is that this case stops breaking if the coordinates are instead [(0,0),(-13.3,0),(-6.5,0.15),(-6.8,0.15)], I only mention this for identifying the bug, not as a workaround

[Kemikal1]

Ah, ok i see.

[fire]

@smix8 I remember you had a lot of comments and maybe some fixes for this.

[smix8]

I don't have the time to investigate this but my blind guess would be that this is either an unsolved bug in Godot 3 cause the Delaunay code in Godot 4 is slightly different, or this is just another precision issues as there are a lot of float ops in the Delaunay function. Also a lot of length_squared() and distance_squared_to() checks which are fast but not anywhere as precise as e.g. MathLab.

[AThousandShips]

Note that this happens on double precision builds as well

[kleonc]

So the issue seems to be the algorithm itself, it's not guaranteed to returna a proper Delaunay triangulation for the given set of N points. Currently used algorithm starts with an arbitrarily choosen bounding triangle and keeps subdividing until a proper triangulation is found for such set of N + 3 points (original N points and 3 bounding vertices). At the end all triangles including at least 1 bounding vertex are discarded from the final result. The issue is: the remaining triangles are not guaranteed to form a valid triangulation of the original set of N points. At worst no triangles will be left, like in the MRP (before the removal algorithm ends with 9 triangles but all of them are discarded).

Possibly it fails only for cocircular points (all lying on the same circle), like the 4 points from the MRP. But not sure about it.

Anyway, the algorithm needs tweaking/rewrite. I'll look into it.

[MatthewLJensen]

I just created a pull request for a bug I initially assumed to be the same as the one discussed here. But after trying the problematic example provided in this issue, I realized that my bug was actually different. In some further testing, it looks like the Bowyer Watson algorithm has problems with nearly colinear points, and thus the bounding triangle needs to potentially be infinite to handle these scenarios. This is because the circumscribed circle grows rapidly as the points approach a line. Thus, the bounding triangle technically needs to be large enough to avoid having one of its vertices enclosed by that circle.

I tried increasing the multiplier when creating the bounding triangle from 16 to 100, and I believe I was able to reduce the non-zero y-coordinates in the example in this issue down to -0.045 while still getting valid results. I'm running out of time to work on this right now, but I'm planning on creating a graphic in desmos later tonight to show the problem in more detail.

[MatthewLJensen]

Further evidence to indicate that the problem is with a bounding triangle that is too small. The circle at the top is the circumcircle of three of the points provided by @Kemikal1 [(0,0), (-6.5,0.15), (-13.3,0)]. The blue triangle is the bounding triangle created by the algorithm. Notice that the vertex of the bounding triangle is found within the circumcircle of the points. I believe this is why the algorithm is breaking down.

As noted by @AThousandShips , if the y-values are positive, the Delaunay triangulation does work, but I think this is because it effectively flips the circumcircle (circle in red) which then no longer contains any bounding triangle vertices.

See the Desmos graph.

I'm still not sure the best way to fix this algorithm. In the stack exchange forum post I mentioned above, it looks like you can theoretically set the bounding triangle to an infinite size, but I'm pretty sure it would involve some more maths to handle situations when some circumcircles are consequently infinite. I'm surprised I don't see more details about this problem with the Bowyer Watson algorithm other places online.

I'm also wondering if we could add a check prior to the circumcircle check on line 116, such that if all three points of the considered triangle are non-bounding-triangle vertices, and if the current point i is, then continue. I haven't had a chance to test this hypothesis yet.

[AThousandShips]

That's right! I was pondering why, though it isn't necessary relevant that it does contain them, but I suspect the issue is that flipping with nearly collinear points creates a new degenerate case

One method could be to remap the points, this would theoretically improve numerical stability as well, instead of scaling the super-triangle to large numbers, instead scale down the vertices to the range (0,1] (might be useful to exclude 0, unsure), and then it should be possible to pick super-vertices such that they prevent this issue, tested some when investigating this but didn't complete it

Been using and working with some of this, using CGAL, which isn't compatible licence wise, but maybe their bibliography can help, and I'll do some testing to see what I can come up with and share if I do

[Kemikal1]

, but I'm pretty sure it would involve some more maths to handle situations when some circumcircles are consequently infinite.

The only time a circumcircle is infinite is if the points are colinear. So they dont need to be considered. My two cents is to make the bounding riangle as big as the error allows it, if you dont want to consider an infinite bounding triangle. There must be some formula that makes a bounding triangle from the smallest error considered(im talking about angle error).

[kleonc]

Assigned myself to this issue to not forget about it again (I remember getting into it after https://github.com/godotengine/godot/issues/77279#issuecomment-1560147702 but I guess I focused on something else later and forgot about this issue). I'm not getting straight into it right now so if anyone else wants to do so - feel free to do so. Setting milestone to 4.2 as a kinda deadline/reminder for myself.

[YuriSizov]

I'll set it to 4.3 as a reminder for you, @kleonc :P But if this is not a priority, we should probably put it in 4.x, or remove it entirely until the fix becomes available.

[mr11stein]

I think the correct way is to treat the vertices of the super triangle with infinites and then handle the circumcircle with those vertices as an edge case, which is hinted at in this answer: https://math.stackexchange.com/a/4001825

[christianzski]

I think the correct way is to treat the vertices of the super triangle with infinites and then handle the circumcircle with those vertices as an edge case, which is hinted at in this answer: https://math.stackexchange.com/a/4001825

Creating an infinite super triangle will most definitely fix this issue. I encountered this same issue while implementing Bowyer-Watson in C++. You might fix this specific case by scaling the super triangle more, but you'll just find more cases where this breaks down as 3 points approach collinearity.

Not sure if it is frowned upon to link my own code, but for reference you can see how I dealt with the edge cases introduced by an infinite super triangle here. Mostly it involves detecting whether a new point is inside a half-plane or not, but this will incur some performance hit due to more branching.

[fire]

I would help review and support a pull request to godot engine for detecting whether a new point is inside a half-plane or not, but this will incur some performance hit due to more branching. or similar solutions.