Implement ComputeTriangleAreas, GetNonManifoldEdges and RemoveNonMani…

…foldEdges in t::geometry::TriangleMesh (#6657)

Split the logic from ComputeSurfaceArea into a helper static function
and introduce a new method which computes triangle areas and writes the
resulting tensor into attributes of the mesh.

t::geometry::TriangleMesh::GetNonManifoldEdges mimics the logic of the
legacy method.

t::geometry::TriangleMesh::RemoveNonManifoldEdges follows the logic of
the legacy method but there are a few differences:
* the main difference is that the outer while-loop is removed. I don't
  see how after the first iteration any edge can have more than 2
  adjacent triangles, which makes the further iterations unnecessary.
* I count triangles with non-negative areas immediately and do not rely
  on the total number of adjacent triangles (which would also include
  triangles marked for removal).
* To choose a triangle with the minimal area out of the existing
  adjacent triangles I use a heap structure.

* Use unordered / sorted comparison for GetNonManifoldEdges() test, return empty areas property if there are no triangles.
---------
Co-authored-by: Sameer Sheorey <sameer.sheorey@intel.com>

cpp/open3d/t/geometry/TriangleMesh.cpp

@@ -283,6 +283,45 @@ TriangleMesh &TriangleMesh::ComputeVertexNormals(bool normalized) {
     return *this;
 }
 
+static core::Tensor ComputeTriangleAreasHelper(const TriangleMesh &mesh) {
+    const int64_t triangle_num = mesh.GetTriangleIndices().GetLength();
+    const core::Dtype dtype = mesh.GetVertexPositions().GetDtype();
+    core::Tensor triangle_areas({triangle_num}, dtype, mesh.GetDevice());
+    if (mesh.IsCPU()) {
+        kernel::trianglemesh::ComputeTriangleAreasCPU(
+                mesh.GetVertexPositions().Contiguous(),
+                mesh.GetTriangleIndices().Contiguous(), triangle_areas);
+    } else if (mesh.IsCUDA()) {
+        CUDA_CALL(kernel::trianglemesh::ComputeTriangleAreasCUDA,
+                  mesh.GetVertexPositions().Contiguous(),
+                  mesh.GetTriangleIndices().Contiguous(), triangle_areas);
+    } else {
+        utility::LogError("Unimplemented device");
+    }
+
+    return triangle_areas;
+}
+
+TriangleMesh &TriangleMesh::ComputeTriangleAreas() {
+    if (IsEmpty()) {
+        utility::LogWarning("TriangleMesh is empty.");
+        return *this;
+    }
+
+    if (!HasTriangleIndices()) {
+        SetTriangleAttr("areas", core::Tensor::Empty(
+                                         {0}, GetVertexPositions().GetDtype(),
+                                         GetDevice()));
+        utility::LogWarning("TriangleMesh has no triangle indices.");
+        return *this;
+    }
+
+    core::Tensor triangle_areas = ComputeTriangleAreasHelper(*this);
+    SetTriangleAttr("areas", triangle_areas);
+
+    return *this;
+}
+
 double TriangleMesh::GetSurfaceArea() const {
     double surface_area = 0;
     if (IsEmpty()) {
@@ -295,22 +334,7 @@ double TriangleMesh::GetSurfaceArea() const {
         return surface_area;
     }
 
-    const int64_t triangle_num = GetTriangleIndices().GetLength();
-    const core::Dtype dtype = GetVertexPositions().GetDtype();
-    core::Tensor triangle_areas({triangle_num}, dtype, GetDevice());
-
-    if (IsCPU()) {
-        kernel::trianglemesh::ComputeTriangleAreasCPU(
-                GetVertexPositions().Contiguous(),
-                GetTriangleIndices().Contiguous(), triangle_areas);
-    } else if (IsCUDA()) {
-        CUDA_CALL(kernel::trianglemesh::ComputeTriangleAreasCUDA,
-                  GetVertexPositions().Contiguous(),
-                  GetTriangleIndices().Contiguous(), triangle_areas);
-    } else {
-        utility::LogError("Unimplemented device");
-    }
-
+    core::Tensor triangle_areas = ComputeTriangleAreasHelper(*this);
     surface_area = triangle_areas.Sum({0}).To(core::Float64).Item<double>();
 
     return surface_area;
@@ -1326,6 +1350,163 @@ TriangleMesh TriangleMesh::RemoveUnreferencedVertices() {
     return *this;
 }
 
+template <typename T,
+          typename std::enable_if<std::is_integral<T>::value &&
+                                          !std::is_same<T, bool>::value,
+                                  T>::type * = nullptr>
+using Edge = std::tuple<T, T>;
+
+/// brief Helper function to get an edge with ordered vertex indices.
+template <typename T>
+static inline Edge<T> GetOrderedEdge(T vidx0, T vidx1) {
+    return (vidx0 < vidx1) ? Edge<T>{vidx0, vidx1} : Edge<T>{vidx1, vidx0};
+}
+
+/// brief Helper
+///
+template <typename T>
+static std::unordered_map<Edge<T>,
+                          std::vector<size_t>,
+                          utility::hash_tuple<Edge<T>>>
+GetEdgeToTrianglesMap(const core::Tensor &tris_cpu) {
+    std::unordered_map<Edge<T>, std::vector<size_t>,
+                       utility::hash_tuple<Edge<T>>>
+            tris_per_edge;
+    auto AddEdge = [&](T vidx0, T vidx1, int64_t tidx) {
+        tris_per_edge[GetOrderedEdge(vidx0, vidx1)].push_back(tidx);
+    };
+    const T *tris_ptr = tris_cpu.GetDataPtr<T>();
+    for (int64_t tidx = 0; tidx < tris_cpu.GetLength(); ++tidx) {
+        const T *triangle = &tris_ptr[3 * tidx];
+        AddEdge(triangle[0], triangle[1], tidx);
+        AddEdge(triangle[1], triangle[2], tidx);
+        AddEdge(triangle[2], triangle[0], tidx);
+    }
+    return tris_per_edge;
+}
+
+TriangleMesh TriangleMesh::RemoveNonManifoldEdges() {
+    if (!HasVertexPositions() || GetVertexPositions().GetLength() == 0) {
+        utility::LogWarning(
+                "[RemoveNonManifildEdges] TriangleMesh has no vertices.");
+        return *this;
+    }
+
+    if (!HasTriangleIndices() || GetTriangleIndices().GetLength() == 0) {
+        utility::LogWarning(
+                "[RemoveNonManifoldEdges] TriangleMesh has no triangles.");
+        return *this;
+    }
+
+    GetVertexAttr().AssertSizeSynchronized();
+    GetTriangleAttr().AssertSizeSynchronized();
+
+    core::Tensor tris_cpu =
+            GetTriangleIndices().To(core::Device()).Contiguous();
+
+    ComputeTriangleAreas();
+    core::Tensor tri_areas_cpu =
+            GetTriangleAttr("areas").To(core::Device()).Contiguous();
+
+    DISPATCH_FLOAT_INT_DTYPE_TO_TEMPLATE(
+            GetVertexPositions().GetDtype(), tris_cpu.GetDtype(), [&]() {
+                scalar_t *tri_areas_ptr = tri_areas_cpu.GetDataPtr<scalar_t>();
+                auto edges_to_tris = GetEdgeToTrianglesMap<int_t>(tris_cpu);
+
+                // lambda to compare triangles areas by index
+                auto area_greater_compare = [&tri_areas_ptr](size_t lhs,
+                                                             size_t rhs) {
+                    return tri_areas_ptr[lhs] > tri_areas_ptr[rhs];
+                };
+
+                // go through all edges and for those that have more than 2
+                // triangles attached, remove the triangles with the minimal
+                // area
+                for (auto &kv : edges_to_tris) {
+                    // remove all triangles which are already marked for removal
+                    // (area < 0) note, the erasing of triangles happens
+                    // afterwards
+                    auto tris_end = std::remove_if(
+                            kv.second.begin(), kv.second.end(),
+                            [=](size_t t) { return tri_areas_ptr[t] < 0; });
+                    // count non-removed triangles (with area > 0).
+                    int n_tris = std::distance(kv.second.begin(), tris_end);
+
+                    if (n_tris <= 2) {
+                        // nothing to do here as either:
+                        // - all triangles of the edge are already marked for
+                        // deletion
+                        // - the edge is manifold: it has 1 or 2 triangles with
+                        //   a non-negative area
+                        continue;
+                    }
+
+                    // now erase all triangle indices already marked for removal
+                    kv.second.erase(tris_end, kv.second.end());
+
+                    // find first to triangles with the maximal area
+                    std::nth_element(kv.second.begin(), kv.second.begin() + 1,
+                                     kv.second.end(), area_greater_compare);
+
+                    // mark others for deletion
+                    for (auto it = kv.second.begin() + 2; it < kv.second.end();
+                         ++it) {
+                        tri_areas_ptr[*it] = -1;
+                    }
+                }
+            });
+
+    // mask for triangles with positive area
+    core::Tensor tri_mask = tri_areas_cpu.Gt(0.0).To(GetDevice());
+
+    // pick up positive-area triangles (and their attributes)
+    for (auto item : GetTriangleAttr()) {
+        SetTriangleAttr(item.first, item.second.IndexGet({tri_mask}));
+    }
+
+    return *this;
+}
+
+core::Tensor TriangleMesh::GetNonManifoldEdges(
+        bool allow_boundary_edges /* = true */) const {
+    if (!HasVertexPositions()) {
+        utility::LogWarning(
+                "[GetNonManifoldEdges] TriangleMesh has no vertices.");
+        return {};
+    }
+
+    if (!HasTriangleIndices()) {
+        utility::LogWarning(
+                "[GetNonManifoldEdges] TriangleMesh has no triangles.");
+        return {};
+    }
+
+    core::Tensor result;
+    core::Tensor tris_cpu =
+            GetTriangleIndices().To(core::Device()).Contiguous();
+    core::Dtype tri_dtype = tris_cpu.GetDtype();
+
+    DISPATCH_INT_DTYPE_PREFIX_TO_TEMPLATE(tri_dtype, tris, [&]() {
+        auto edges = GetEdgeToTrianglesMap<scalar_tris_t>(tris_cpu);
+        std::vector<scalar_tris_t> non_manifold_edges;
+
+        for (auto &kv : edges) {
+            if ((allow_boundary_edges &&
+                 (kv.second.size() < 1 || kv.second.size() > 2)) ||
+                (!allow_boundary_edges && kv.second.size() != 2)) {
+                non_manifold_edges.push_back(std::get<0>(kv.first));
+                non_manifold_edges.push_back(std::get<1>(kv.first));
+            }
+        }
+
+        result = core::Tensor(non_manifold_edges,
+                              {(long int)non_manifold_edges.size() / 2, 2},
+                              tri_dtype, GetTriangleIndices().GetDevice());
+    });
+
+    return result;
+}
+
 }  // namespace geometry
 }  // namespace t
 }  // namespace open3d

cpp/open3d/t/geometry/TriangleMesh.h

@@ -671,6 +671,11 @@ class TriangleMesh : public Geometry, public DrawableGeometry {
     /// of the individual triangle surfaces.
     double GetSurfaceArea() const;
 
+    /// \brief Function to compute triangle areas and save it as a triangle
+    /// attribute "areas". Prints a warning, if mesh is empty or has no
+    /// triangles.
+    TriangleMesh &ComputeTriangleAreas();
+
     /// \brief Clip mesh with a plane.
     /// This method clips the triangle mesh with the specified plane.
     /// Parts of the mesh on the positive side of the plane will be kept and
@@ -974,6 +979,20 @@ class TriangleMesh : public Geometry, public DrawableGeometry {
     /// \return The reference to itself.
     TriangleMesh RemoveUnreferencedVertices();
 
+    /// Removes all non-manifold edges, by successively deleting triangles
+    /// with the smallest surface area adjacent to the
+    /// non-manifold edge until the number of adjacent triangles to the edge is
+    /// `<= 2`. If mesh is empty or has no triangles, prints a warning and
+    /// returns immediately. \return The reference to itself.
+    TriangleMesh RemoveNonManifoldEdges();
+
+    /// Returns the non-manifold edges of the triangle mesh.
+    /// If \param allow_boundary_edges is set to false, then also boundary
+    /// edges are returned.
+    /// \return 2d integer tensor with shape {n,2} encoding ordered edges.
+    /// If mesh is empty or has no triangles, returns an empty tensor.
+    core::Tensor GetNonManifoldEdges(bool allow_boundary_edges = true) const;
+
 protected:
     core::Device device_ = core::Device("CPU:0");
     TensorMap vertex_attr_;

cpp/pybind/t/geometry/trianglemesh.cpp

@@ -985,9 +985,42 @@ or has a negative value, it is ignored.
         top_face = box.select_by_index([2, 3, 6, 7])
 )");
 
-    triangle_mesh.def("remove_unreferenced_vertices",
-                      &TriangleMesh::RemoveUnreferencedVertices,
-                      "Removes unreferenced vertices from the mesh in-place.");
+    triangle_mesh.def(
+            "remove_unreferenced_vertices",
+            &TriangleMesh::RemoveUnreferencedVertices,
+            R"(Removes unreferenced vertices from the mesh in-place.)");
+
+    triangle_mesh.def(
+            "compute_triangle_areas", &TriangleMesh::ComputeTriangleAreas,
+            R"(Compute triangle areas and save it as \"areas\" triangle attribute.
+
+Returns:
+    The mesh.
+
+Example:
+
+    This code computes the overall surface area of a box:
+
+        import open3d as o3d
+        box = o3d.t.geometry.TriangleMesh.create_box()
+        surface_area = box.compute_triangle_areas().triangle.areas.sum()
+)");
+
+    triangle_mesh.def("remove_non_manifold_edges",
+                      &TriangleMesh::RemoveNonManifoldEdges,
+                      R"(Function that removes all non-manifold edges, by
+successively deleting  triangles with the smallest surface
+area adjacent to the non-manifold edge until the number of
+adjacent triangles to the edge is `<= 2`.
+
+Returns:
+    The mesh.
+)");
+
+    triangle_mesh.def("get_non_manifold_edges",
+                      &TriangleMesh::GetNonManifoldEdges,
+                      "allow_boundary_edges"_a = true,
+                      R"(Returns the list consisting of non-manifold edges.)");
 }
 
 }  // namespace geometry