Triangulating a 3D mesh from a set of data points

Question

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:

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...