Delaunay3D can produce overlapping tetrahedra

[permelin]

Tested versions

4.3.dev.custom_build.29b3d9e9e

System information

Godot v4.3.dev (248e58a4e) - Debian GNU/Linux trixie/sid trixie - X11 - Vulkan (Forward+) - dedicated NVIDIA GeForce RTX 2080 (nvidia) - Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz (8 Threads)

Issue description

Delaunay3D::tetrahedralize() when called directly from C++ or from GDScript via Geometry3D.tetrahedralize_delaunay() can produce overlapping tetrahedra (with points inside their circumsphere). All examples of two overlapping tetrahedra that I've seen have had two or three vertices in common and intersecting volumes.

It doesn't happen all the time and I don't know exactly what circumstances trigger it, but it seems that both of these must be true:

• there are some points that are several hundreds distance units apart, and
• there are some points that are only a couple of units apart

I noticed this problem while working on issues with LightmapGI probe triangulation and partitioning.

It's difficult to create an MRP that doesn't involve thousands of points and doesn't require extra tooling to find and visualize the few invalid cases. The attached MRP is the smallest and simplest thing I can do.

The MRP is just one scene with one node and one script. I'll dump the script here so you can check my code for obvious bugs without running the project.

@tool
extends MeshInstance3D

const points: PackedVector3Array = [
	Vector3(204.537, 2.64082, 0),
	Vector3(208.53, 2.64082, 0),
	Vector3(204.537, 2.64082, 2),
	Vector3(206.53, 2.64082, 2),
	Vector3(-501.4156, -2.415572, -501.4156),
	Vector3(501.4156, -2.415572, -501.4156),
	Vector3(-501.4156, 248.2922, -501.4156),
	Vector3(501.4156, 248.2922, -501.4156),
	Vector3(-501.4156, -2.415572, 501.4156),
	Vector3(501.4156, -2.415572, 501.4156),
	Vector3(-501.4156, 248.2922, 501.4156),
	Vector3(501.4156, 248.2922, 501.4156),
	Vector3(-250.7078, -2.415572, -501.4156),
	Vector3(-250.7078, 248.2922, -501.4156),
	Vector3(-501.4156, -2.415572, -250.7078),
	Vector3(-250.7078, -2.415572, -250.7078),
	Vector3(-501.4156, 248.2922, -250.7078),
	Vector3(-250.7078, 248.2922, -250.7078),
	Vector3(0, -2.415572, -501.4156),
	Vector3(0, 248.2922, -501.4156),
	Vector3(0, -2.415572, -250.7078),
	Vector3(0, 248.2922, -250.7078),
	Vector3(-501.4156, -2.415572, 0),
	Vector3(-250.7078, -2.415572, 0),
	Vector3(-501.4156, 248.2922, 0),
	Vector3(-250.7078, 248.2922, 0),
	Vector3(0, -2.415572, 0),
	Vector3(0, 248.2922, 0),
	Vector3(250.7078, -2.415572, -501.4156),
	Vector3(250.7078, 248.2922, -501.4156),
	Vector3(250.7078, -2.415572, -250.7078),
	Vector3(250.7078, 248.2922, -250.7078),
	Vector3(501.4156, -2.415572, -250.7078),
	Vector3(501.4156, 248.2922, -250.7078),
	Vector3(250.7078, -2.415572, 0),
	Vector3(250.7078, 248.2922, 0),
	Vector3(501.4156, -2.415572, 0),
	Vector3(501.4156, 248.2922, 0),
	Vector3(-501.4156, -2.415572, 250.7078),
	Vector3(-250.7078, -2.415572, 250.7078),
	Vector3(-501.4156, 248.2922, 250.7078),
	Vector3(-250.7078, 248.2922, 250.7078),
	Vector3(0, -2.415572, 250.7078),
	Vector3(0, 248.2922, 250.7078),
	Vector3(-250.7078, -2.415572, 501.4156),
	Vector3(-250.7078, 248.2922, 501.4156),
	Vector3(0, -2.415572, 501.4156),
	Vector3(0, 248.2922, 501.4156),
	Vector3(250.7078, -2.415572, 250.7078),
	Vector3(250.7078, 248.2922, 250.7078),
	Vector3(501.4156, -2.415572, 250.7078),
	Vector3(501.4156, 248.2922, 250.7078),
	Vector3(250.7078, -2.415572, 501.4156),
	Vector3(250.7078, 248.2922, 501.4156)]

var imesh: ImmediateMesh

func _ready() -> void:
	imesh = ImmediateMesh.new()
	mesh = imesh

	var tetrahedra: PackedInt32Array = Geometry3D.tetrahedralize_delaunay(points)
	assert(tetrahedra.size() % 4 == 0)

	# These two are just one example of overlapping tetrahedra
	draw_tetrahedron(tetrahedra, 129, Color.RED)
	draw_tetrahedron(tetrahedra, 130, Color.BLUE)

func draw_tetrahedron(tetrahedra: PackedInt32Array, index: int, color: Color) -> void:
	var vertices: Array[Vector3] = [
		points[tetrahedra[index * 4]],
		points[tetrahedra[index * 4 + 1]],
		points[tetrahedra[index * 4 + 2]],
		points[tetrahedra[index * 4 + 3]]
	]
	print(vertices)

	imesh.surface_begin(Mesh.PRIMITIVE_TRIANGLES)
	for i in 4:
		imesh.surface_set_color(color)
		imesh.surface_add_vertex(vertices[i])
		imesh.surface_add_vertex(vertices[(i + 1) % 4])
		imesh.surface_add_vertex(vertices[(i + 2) % 4])
	imesh.surface_end()

The two tetrahedra chosen here, and in the screenshot above, are:
[(250.7078, 248.2922, 0), (501.4156, 248.2922, 0), (250.7078, -2.415572, 250.7078), (250.7078, 248.2922, 250.7078)]
and:
[(250.7078, 248.2922, 0), (501.4156, -2.415572, 0), (250.7078, -2.415572, 250.7078), (250.7078, 248.2922, 250.7078)]

I just noticed that in this particular example, the two have exactly the same circumsphere, but that is likely a red herring since it is not true for many other cases (e.g. 117 together with 167).

Steps to reproduce

Load the MRP. It's a tool script and will show the issue directly in the editor.

Minimal reproduction project (MRP)

invalid-delaunay3d-mrp.zip

[permelin]

The problem is the constant used here:

godot/core/math/delaunay_3d.h

Line 176 in f47f4a0

return radius2 < (p_simplex.circum_r2 - R128(0.00001));

When distances between points span a large range, from 2 to 350 units (meters) in my experiments, this constant is too big. Making it smaller by two orders of magnitude improves the situation a lot, but at the price of introducing some instability for smaller point ranges. I don't know where this instability comes from. Going to 1e-7 is pushing beyond the significant digits of Vector3, but the relevant calculations are all done with R128.