# Triangulating a 3D mesh from a set of data points

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);
triangles.push_back(tetrahedron->faces);
triangles.push_back(tetrahedron->faces);
triangles.push_back(tetrahedron->faces);
}
}

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