jmespadero / pyDelaunay2D

A simple Delaunay 2D triangulation in python (with numpy)
GNU General Public License v3.0
186 stars 34 forks source link

From 2D to 3D #5

Open henryhansen23 opened 3 years ago

henryhansen23 commented 3 years ago

Hi I'm want to create Delaunay triangulation in 3D. I find your implementation interesting. Have you been giving any thought on what this line could look like in 3D: boundary.append((T[(edge+1) % 3], T[(edge-1) % 3], tri_op))?

jmespadero commented 3 years ago

First, I must remember that this implementation is done with didactic purposes. If you really need to use 3D triangulation (a better name is 3D tetrahedralization), consider using scipy.spatial.Delaunay. If you really want to learn how to implement 3D Delaunay, then go ahead.

First consideration is that in 3D you don't work with triangles, but with tetrahedrons (or simplex in 3D), so some code that in 2D is as simple as:

# Create two CCW triangles for the frame
        T1 = (0, 1, 3)
        T2 = (2, 3, 1)

must be rewritten as something quite more complex, like:

# Create five tetrahedrons for the cubic frame
        T1 = (1, 7, 4, 2)
        T2 = (1, 4, 7, 5)
        T3 = (1, 2, 4, 0)
        T4 = (1, 2, 7, 3)
        T5 = (7, 2, 4, 6)

Where the indexes come from the tetrahedralization of cubic initial frame. image

The line in 2D code:

boundary.append((T[(edge+1) % 3], T[(edge-1) % 3], tri_op))

Can be interpreted as: Insert a tuple formed by one edge and the index of its opposite triangle in the boundary. The T[(edge+1) % 3], T[(edge-1) % 3] formula means: Choose the next and previous vertex as extremes of this edge, which also could be interpreted as use the vertices of the triangle T that are not in the "edge" position and change the orientation (ie, new edge must be (next, previous), not (previous, next), because we want that edge to have the opposite orientation that has in the triangle T

In the 3D approach, you have to insert a new triangle (instead edge) using the three vertices that are not the "edge" position, sorted in the opposite orientation that in the current tetrahedron, and the index of the opposite tetrahedron (instead triangle) in the boundary. This can be coded as:

boundary.append((T[(edge+1) % 4], T[(edge-1) % 4], T[(edge+2) % 4], tetra_op))

Anyway, you are going to run into more problems just after this line, because there is not an easy way to «Check if boundary is a closed loop» as easy as done in 2D.

henryhansen23 commented 3 years ago

Terrific answer! Okay, so I was onto something then, but you are on top of stuff) I think I got your thorough explanations, thank you. "«Check if boundary is a closed loop» is not as easy as done in 2D." Right, but not impossible I guess. However, maybe I have to think differently since the complexity increased more than I thought when going from 2D to 3D. I had read that the Watson-algorithm could quite easily be extended to dimensions beyond 2D, so I was optimistic about this. I wish my life was so simple so that I could just use scipy.spatial.Delaunay.

jmespadero commented 3 years ago

I will let to others to judge if it easy or not to implement the multidimensional version of Bowyer-Watson algorithm because that depends of your experience in programming geometry algorithms and «understand dimensions beyond 3D».

The «Check if boundary is a closed loop» part can be avoided if you just iterate every triangle//tetrahedron in the bad_triangles list. The real problem will come later, because we have to «Retriangle the hole left by removal of bad_triangles» and then «Link the new triangles each another». This is trivial if our boundary is sorted in CCW order, because each triangle is linked with previous and next triangles in new_triangles list. If not, you have to do an extra search (find pairs of tetrahedrons that share 3 vertex)

This kind of complexity problems going to upper dimensions make the algorithm being not-so-didactic. Most people just can't imagine one «star-shaped hole in space being filled with tetrahedrons, while maintaining the coherence of orientation of the boundary».

Also... Could you share why you can't use scipy.spatial.Delaunay or QHull libraries? You can email me if you prefer to keep your reasons private.

henryhansen23 commented 3 years ago

That's very useful so thank you very much for you comments. I will take your thoughts into account and do some more research before would I go further.

"Most people just can't imagine one «star-shaped hole in space being filled with tetrahedrons, while maintaining the coherence of orientation of the boundary»".That's definitely true...