Feature/Geometry type for tetrahedral meshes (#1109)

* Delaunay tetramesh (#2)

Added new geometry type TetraMesh.

* added missing headers

* removed superfluous check

* typo and style fixes

* added references/description
-qhull for tetrahedralization
-primal contouring for surface extraction

* renamed ComputeDelaunayTriangulation3D ...
to ComputeDelaunayTetrahedralization

* style fix

src/Open3D/Geometry/Geometry.h

@@ -41,6 +41,7 @@ class Geometry {
         HalfEdgeTriangleMesh = 6,
         Image = 7,
         RGBDImage = 8,
+        TetraMesh = 9,
     };
 
 public:

src/Open3D/Geometry/LineSet.h

@@ -37,6 +37,7 @@ namespace geometry {
 
 class PointCloud;
 class TriangleMesh;
+class TetraMesh;
 
 class LineSet : public Geometry3D {
 public:
@@ -94,6 +95,10 @@ class LineSet : public Geometry3D {
     static std::shared_ptr<LineSet> CreateFromTriangleMesh(
             const TriangleMesh &mesh);
 
+    /// Factory function to create a LineSet from edges of a tetra mesh
+    /// \param mesh.
+    static std::shared_ptr<LineSet> CreateFromTetraMesh(const TetraMesh &mesh);
+
 public:
     std::vector<Eigen::Vector3d> points_;
     std::vector<Eigen::Vector2i> lines_;

src/Open3D/Geometry/LineSetFactory.cpp

@@ -28,6 +28,7 @@
 
 #include "Open3D/Geometry/LineSet.h"
 #include "Open3D/Geometry/PointCloud.h"
+#include "Open3D/Geometry/TetraMesh.h"
 #include "Open3D/Geometry/TriangleMesh.h"
 
 namespace open3d {
@@ -79,5 +80,31 @@ std::shared_ptr<LineSet> LineSet::CreateFromTriangleMesh(
     return line_set;
 }
 
+std::shared_ptr<LineSet> LineSet::CreateFromTetraMesh(const TetraMesh &mesh) {
+    auto line_set = std::make_shared<LineSet>();
+    line_set->points_ = mesh.vertices_;
+
+    std::unordered_set<Eigen::Vector2i,
+                       utility::hash_eigen::hash<Eigen::Vector2i>>
+            inserted_edges;
+    auto InsertEdge = [&](int vidx0, int vidx1) {
+        Eigen::Vector2i edge(std::min(vidx0, vidx1), std::max(vidx0, vidx1));
+        if (inserted_edges.count(edge) == 0) {
+            inserted_edges.insert(edge);
+            line_set->lines_.push_back(Eigen::Vector2i(vidx0, vidx1));
+        }
+    };
+    for (const auto &tetra : mesh.tetras_) {
+        InsertEdge(tetra(0), tetra(1));
+        InsertEdge(tetra(1), tetra(2));
+        InsertEdge(tetra(2), tetra(0));
+        InsertEdge(tetra(3), tetra(0));
+        InsertEdge(tetra(3), tetra(1));
+        InsertEdge(tetra(3), tetra(2));
+    }
+
+    return line_set;
+}
+
 }  // namespace geometry
 }  // namespace open3d

src/Open3D/Geometry/Qhull.cpp

@@ -25,7 +25,9 @@
 // ----------------------------------------------------------------------------
 
 #include "Open3D/Geometry/Qhull.h"
+#include "Open3D/Geometry/TetraMesh.h"
 #include "Open3D/Geometry/TriangleMesh.h"
+#include "Open3D/Utility/Console.h"
 
 #include "libqhullcpp/PointCoordinates.h"
 #include "libqhullcpp/Qhull.h"
@@ -97,5 +99,83 @@ std::shared_ptr<TriangleMesh> Qhull::ComputeConvexHull(
     return convex_hull;
 }
 
+std::shared_ptr<TetraMesh> Qhull::ComputeDelaunayTetrahedralization(
+        const std::vector<Eigen::Vector3d>& points) {
+    typedef decltype(TetraMesh::tetras_)::value_type Vector4i;
+    auto delaunay_triangulation = std::make_shared<TetraMesh>();
+
+    if (points.size() < 4) {
+        utility::LogWarning(
+                "[ComputeDelaunayTriangulation3D] not enough points to create "
+                "a tetrahedral mesh.\n");
+        return delaunay_triangulation;
+    }
+
+    // qhull cannot deal with this case
+    if (points.size() == 4) {
+        delaunay_triangulation->vertices_ = points;
+        delaunay_triangulation->tetras_.push_back(Vector4i(0, 1, 2, 3));
+        return delaunay_triangulation;
+    }
+
+    std::vector<double> qhull_points_data(points.size() * 3);
+    for (size_t pidx = 0; pidx < points.size(); ++pidx) {
+        const auto& pt = points[pidx];
+        qhull_points_data[pidx * 3 + 0] = pt(0);
+        qhull_points_data[pidx * 3 + 1] = pt(1);
+        qhull_points_data[pidx * 3 + 2] = pt(2);
+    }
+
+    orgQhull::PointCoordinates qhull_points(3, "");
+    qhull_points.append(qhull_points_data);
+
+    orgQhull::Qhull qhull;
+    qhull.runQhull(qhull_points.comment().c_str(), qhull_points.dimension(),
+                   qhull_points.count(), qhull_points.coordinates(),
+                   "d Qbb Qt");
+
+    orgQhull::QhullFacetList facets = qhull.facetList();
+    delaunay_triangulation->tetras_.resize(facets.count());
+    std::unordered_map<int, int> vert_map;
+    std::unordered_set<int> inserted_vertices;
+    int tidx = 0;
+    for (orgQhull::QhullFacetList::iterator it = facets.begin();
+         it != facets.end(); ++it) {
+        if (!(*it).isGood()) continue;
+
+        orgQhull::QhullFacet f = *it;
+        orgQhull::QhullVertexSet vSet = f.vertices();
+        int tetra_subscript = 0;
+        for (orgQhull::QhullVertexSet::iterator vIt = vSet.begin();
+             vIt != vSet.end(); ++vIt) {
+            orgQhull::QhullVertex v = *vIt;
+            orgQhull::QhullPoint p = v.point();
+
+            int vidx = p.id();
+            delaunay_triangulation->tetras_[tidx](tetra_subscript) = vidx;
+            tetra_subscript++;
+
+            if (inserted_vertices.count(vidx) == 0) {
+                inserted_vertices.insert(vidx);
+                vert_map[vidx] = int(delaunay_triangulation->vertices_.size());
+                double* coords = p.coordinates();
+                delaunay_triangulation->vertices_.push_back(
+                        Eigen::Vector3d(coords[0], coords[1], coords[2]));
+            }
+        }
+
+        tidx++;
+    }
+
+    for (auto& tetra : delaunay_triangulation->tetras_) {
+        tetra(0) = vert_map[tetra(0)];
+        tetra(1) = vert_map[tetra(1)];
+        tetra(2) = vert_map[tetra(2)];
+        tetra(3) = vert_map[tetra(3)];
+    }
+
+    return delaunay_triangulation;
+}
+
 }  // namespace geometry
 }  // namespace open3d

src/Open3D/Geometry/Qhull.h

@@ -34,11 +34,15 @@ namespace open3d {
 namespace geometry {
 
 class TriangleMesh;
+class TetraMesh;
 
 class Qhull {
 public:
     static std::shared_ptr<TriangleMesh> ComputeConvexHull(
             const std::vector<Eigen::Vector3d>& points);
+
+    static std::shared_ptr<TetraMesh> ComputeDelaunayTetrahedralization(
+            const std::vector<Eigen::Vector3d>& points);
 };
 
 }  // namespace geometry