The triangulations output are not Delaunay triangulations

[Treer]

Edit: Pull Request #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:

• p0 = (-6189.595f, 8209.541f)
• p1 = (-5733.924f, 8823.252f)
• p2 = (-5748.702f, 8231.538f)
• p3 = (-5687.935f, 7984.462f)
• p4 = (-6189.595f, 7709.541f)
• p5 = (-6893.049f, 7682.856f)
• p6 = (-6759.087f, 8353.894f)
• p7 = (-6763.313f, 8504.048f)
• p8 = (-6284.938f, 8771.748f)

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:

[Triangle p6 p5 p0] passes Delauney constraint
[Triangle p2 p1 p0] FAILED Delauney constraint:
    p8 is inside the circumcircle of [Triangle p2 p1 p0], circumcircle (-5984.762, 8533.476) with radiusSquared 146890.2 (r 383.262536487719)
[Triangle p3 p2 p0] passes Delauney constraint
[Triangle p3 p1 p2] passes Delauney constraint
[Triangle p4 p3 p0] passes Delauney constraint
[Triangle p6 p8 p7] FAILED Delauney constraint:
    p0 is inside the circumcircle of [Triangle p6 p8 p7], circumcircle (-6412.698, 8438.779) with radiusSquared 127190.8 (r 356.638253346441)
[Triangle p5 p4 p0] passes Delauney constraint
[Triangle p6 p0 p1] FAILED Delauney constraint:
    p7 is inside the circumcircle of [Triangle p6 p0 p1], circumcircle (-6343.088, 8799.527) with radiusSquared 371644.4 (r 609.626473752576)
    p8 is inside the circumcircle of [Triangle p6 p0 p1], circumcircle (-6343.088, 8799.527) with radiusSquared 371644.4 (r 609.626473752576)
[Triangle p7 p5 p6] passes Delauney constraint
[Triangle p8 p6 p1] FAILED Delauney constraint:
    p0 is inside the circumcircle of [Triangle p8 p6 p1], circumcircle (-5924.076, 7884.329) with radiusSquared 917734.9 (r 957.9848315605)
    p2 is inside the circumcircle of [Triangle p8 p6 p1], circumcircle (-5924.076, 7884.329) with radiusSquared 917734.9 (r 957.9848315605)
    p3 is inside the circumcircle of [Triangle p8 p6 p1], circumcircle (-5924.076, 7884.329) with radiusSquared 917734.9 (r 957.9848315605)
    p4 is inside the circumcircle of [Triangle p8 p6 p1], circumcircle (-5924.076, 7884.329) with radiusSquared 917734.9 (r 957.9848315605)

[Test]
public void DelaunayTriangulation() {

	var points = new Vector2[] {
		new Vector2(-6189.595f, 8209.541f), // p0
		new Vector2(-5733.924f, 8823.252f), // p1
		new Vector2(-5748.702f, 8231.538f), // p2
		new Vector2(-5687.935f, 7984.462f), // p3
		new Vector2(-6189.595f, 7709.541f), // p4
		new Vector2(-6893.049f, 7682.856f), // p5
		new Vector2(-6759.087f, 8353.894f), // p6
		new Vector2(-6763.313f, 8504.048f), // p7
		new Vector2(-6284.938f, 8771.748f), // p8
	}.ToList();

	List<Triangle2D> allTriangles = new List<Triangle2D>();
	List<Triangle2D> actualTriangles = new List<Triangle2D>();

	DelaunayTriangulation graph = new DelaunayTriangulation();
	graph.Triangulate(points);
	graph.GetAllTriangles(allTriangles);

	// exclude triangles generated by the supertriangle
	List<int> trianglesToRemove = graph.DiscardedTriangles;
	for (int i = 0; i < allTriangles.Count; i++) {
		if (!trianglesToRemove.Contains(i)) {
			actualTriangles.Add(allTriangles[i]);
		}
	}
	Assert.AreEqual(10, actualTriangles.Count);

	// check the triangulation is valid
	// In a Delaunay Triangulation, the circumcircle of each triangle contains no other points
	bool allPassed = true;
	foreach (var triangle in actualTriangles) {

		List<string> constraintFailures = new List<string>();
		foreach (var point in points.Where(p => !TriangleContainsVertex(triangle, p))) {

			if (MathUtils.IsPointInsideCircumcircle(triangle.p0, triangle.p1, triangle.p2, point)) {
				CalculateCircumcircle(triangle.p0, triangle.p1, triangle.p2, out Vector2 circumCenter, out float radiusSquared);
				constraintFailures.Add($"{PointToString(point, points)} is inside the circumcircle of {TriangleToString(triangle, points)}, circumcircle {circumCenter} with radiusSquared {radiusSquared} (r {Math.Sqrt(radiusSquared)})");
			}
		}
		if (constraintFailures.Any()) {
			allPassed = false;
			Console.WriteLine($"{TriangleToString(triangle, points)} FAILED Delauney constraint:\r\n    " + string.Join("\r\n    ", constraintFailures));
		} else { 
			Console.WriteLine($"{TriangleToString(triangle, points)} passes Delauney constraint"); 
		}
	}
	Assert.True(allPassed);

}



// =====================
#region Helper functions
// =====================

private bool TriangleContainsVertex(Triangle2D triangle, Vector2 vertex) {
	return triangle.p0 == vertex || triangle.p1 == vertex || triangle.p2 == vertex;
}

private string TriangleToString(Triangle2D trangle, IList<Vector2> points) {

	string v0 = PointToString(trangle.p0, points);
	string v1 = PointToString(trangle.p1, points);
	string v2 = PointToString(trangle.p2, points);
	return $"[Triangle {v0} {v1} {v2}]";
}

private string PointToString(Vector2 point, IList<Vector2> points) {
	int pointIndex = points.IndexOf(point);
	return $"p{pointIndex}";
}

private void CalculateCircumcircle(Vector2 p0, Vector2 p1, Vector2 p2, out Vector2 circumcenter, out float radiusSquared) {
	// https://codefound.wordpress.com/2013/02/21/how-to-compute-a-circumcircle/#more-58
	// https://en.wikipedia.org/wiki/Circumscribed_circle
	var dA = (double)p0.x * p0.x + (double)p0.y * p0.y;
	var dB = (double)p1.x * p1.x + (double)p1.y * p1.y;
	var dC = (double)p2.x * p2.x + (double)p2.y * p2.y;

	var aux1 = dA * ((double)p2.y - p1.y) + dB * ((double)p0.y - p2.y) + dC * ((double)p1.y - p0.y);
	var aux2 = -(dA * ((double)p2.x - p1.x) + dB * ((double)p0.x - p2.x) + dC * ((double)p1.x - p0.x));
	var div = 2 * (p0.x * ((double)p2.y - p1.y) + p1.x * ((double)p0.y - p2.y) + p2.x * ((double)p1.y - p0.y));

	if (div == 0) {
		throw new Exception($"Div is 0: {p0}, {p1}, {p2}");
	}

	var centerX = aux1 / div;
	var centerY = aux2 / div;
	circumcenter = new Vector2((float)centerX, (float)centerY);
	radiusSquared = (float)((centerX - p0.x) * (centerX - p0.x) + (centerY - p0.y) * (centerY - p0.y));
}
#endregion Helper functions

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.

[QThund]

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.

[Treer]

Sloan's algorithm looks damn hard (but so fast), so thanks for publishing and diagramming your work - gives some shoulders to stand on.

[Treer]

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.