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)) {
            // 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)
                    if (*e1 == *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) {
        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...


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.