2
\$\begingroup\$

I'm trying to implement the Bowyer-Watson point insertion version of the Delaunay triangulation algorithm, but in 3D. I previously implemented the 2D version without problems, but when transitioning to 3D using Tetrahedrons instead of Triangles (and triangles instead of edges), I don't get proper triangulation.

For instance, given 8 points defining the corners of a unit cube (-1,-1,-1 through to 1,1,1), I get the following result:

broken cube

A summary of the algorithm:

  • Create a super tetrahedron that encapsulates all of my points and add it to a tetrahedron list.
  • One at a time, check each point against the circumsphere of each tetrahedron in the list, copying its triangles.
    • Delete all tetrahedrons who were intersected
    • Delete all duplicate triangles in the list
    • Using each remaining triangle, form new tetrahedron from triangle's 3 points + the new point
  • After iterating through all points, delete all tetrahedrons that share a vertex with the starting super tetrahedron

    vector<vec3> Geometry_Tools::Delaunay_Triangulate(const vector<vec3>& points)
    {
        vector<Tetrahedron> tetrahedrons;
        vec3 p1, p2, p3, p4;
        Calculate_Supertetrahedron(points, p1, p2, p3, p4);
        tetrahedrons.push_back(Tetrahedron(p1, p2, p3, p4));
    
        // Add the points one at a time, branch tetrahedrons off of it.
        for (auto point = begin(points); point != end(points); point++) {
            // Mark all the tetrahedrons which contain this point for deletion, grab their triangles too.
            vector<Tetrahedron> badTetrahedrons;
            vector<Triangle> triangles;
            for (auto tetrahedron = begin(tetrahedrons); tetrahedron != end(tetrahedrons); tetrahedron++) {
                if (tetrahedron->CS_Contains_Point(*point)) {
                    badTetrahedrons.push_back(*tetrahedron);
                    triangles.push_back(tetrahedron->faces[0]);
                    triangles.push_back(tetrahedron->faces[1]);
                    triangles.push_back(tetrahedron->faces[2]);
                    triangles.push_back(tetrahedron->faces[3]);
                }
            }
    
            // Delete all the invalid tetrahedrons in 1 pass
            tetrahedrons.erase(std::remove_if(begin(tetrahedrons), end(tetrahedrons), [badTetrahedrons](Tetrahedron &tetrahedron) {
                for (auto invalidTetrahedron = begin(badTetrahedrons); invalidTetrahedron != end(badTetrahedrons); invalidTetrahedron++)
                    if (*invalidTetrahedron == tetrahedron)
                        return true;
                return false;
            }), end(tetrahedrons));
    
            // All the duplicate triangles cancel each other out and aren't needed. Find and mark them for deletion. 
            vector<Triangle> badTriangles;
            for (auto e1 = begin(triangles); e1 != end(triangles); e1++) {
                for (auto e2 = begin(triangles); e2 != end(triangles); e2++) {
                    // We are reading the same list twice, so skip if the itterators are the same
                    if (e1 == e2)
                        continue;
                    if (*e1 == *e2) {
                        badTriangles.push_back(*e1);
                        badTriangles.push_back(*e2);
                    }
                }
            }
    
            // Delete all the invalid triangles in 1 pass. Creates the hole we need to insert the point.
            triangles.erase(std::remove_if(begin(triangles), end(triangles), [badTriangles](Triangle &triangle) {
                for (auto invalidTriangle = begin(badTriangles); invalidTriangle != end(badTriangles); invalidTriangle++) 
                    if (*invalidTriangle == triangle)
                        return true;            
                return false;
            }), end(triangles));
    
            // What remains is an triangle which when used with point == new tetrahedron
            for each (const auto &triangle in triangles)
                tetrahedrons.push_back(Tetrahedron(triangle.p0, triangle.p1, triangle.p2, *point));
        }
    
        // Delete the contribution of the original super tetrahedron.
        // If a tetrahedron shares a point with the original super tetrahedron, then it should be deleted
        tetrahedrons.erase(std::remove_if(begin(tetrahedrons), end(tetrahedrons), [p1, p2, p3, p4](Tetrahedron &tetrahedron) {
            return tetrahedron.Contains_Vertex(p1) || tetrahedron.Contains_Vertex(p2) || tetrahedron.Contains_Vertex(p3) || tetrahedron.Contains_Vertex(p4);
        }), end(tetrahedrons));
    
        vector<vec3> vertices;
        vertices.reserve(tetrahedrons.size() * 3);
        for each (const auto &tetrahedron in tetrahedrons) {
            for each (const auto &triangle in tetrahedron.faces) {
                vertices.push_back(triangle.p0);
                vertices.push_back(triangle.p1);
                vertices.push_back(triangle.p2);
            }
        }   
        return vertices;
    }
    

Can anyone help me solve my problem?

I've verified that I'm calculating the circumsphere correctly; I double checked the formula and my implementation & visualizing the sphere against the tetrahedron.

It seems like most of the final tetrahedrons share 2-3 points with the boundary super tetrahedron. I don't know why so few new tetrahedra between the data points exist...

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.