Introduce t::geometry::ComputeTriangleAreas method

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.

cpp/open3d/t/geometry/TriangleMesh.cpp

@@ -283,6 +283,42 @@ 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()) {
+        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 +331,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;

cpp/open3d/t/geometry/TriangleMesh.h

@@ -669,6 +669,10 @@ 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. 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

cpp/pybind/t/geometry/trianglemesh.cpp

@@ -962,6 +962,10 @@ or has a negative value, it is ignored.
         box = o3d.t.geometry.TriangleMesh.create_box()
         top_face = box.select_by_index([2, 3, 6, 7])
 )");
+
+    triangle_mesh.def("compute_triangle_areas",
+                      &TriangleMesh::ComputeTriangleAreas,
+                      "Compute triangle areas and save it as a mesh attribute");
 }
 
 }  // namespace geometry

cpp/tests/t/geometry/TriangleMesh.cpp

@@ -1212,5 +1212,20 @@ TEST_P(TriangleMeshPermuteDevices, SelectByIndex) {
             box_untouched.GetTriangleIndices()));
 }
 
+TEST_P(TriangleMeshPermuteDevices, ComputeTriangleAreas) {
+    core::Device device = GetParam();
+    t::geometry::TriangleMesh mesh_empty;
+    EXPECT_NO_THROW(mesh_empty.ComputeTriangleAreas());
+
+    std::shared_ptr<open3d::geometry::TriangleMesh> mesh =
+            open3d::geometry::TriangleMesh::CreateSphere(1.0, 3);
+    t::geometry::TriangleMesh t_mesh = t::geometry::TriangleMesh::FromLegacy(
+            *mesh, core::Float64, core::Int64, device);
+    std::vector<double> areas;
+    mesh->GetSurfaceArea(areas);
+    t_mesh.ComputeTriangleAreas();
+    EXPECT_TRUE(t_mesh.GetTriangleAttr("areas").AllClose(
+            core::Tensor(areas, {(int)areas.size()}, core::Float64)));
+}
 }  // namespace tests
 }  // namespace open3d

python/test/t/geometry/test_trianglemesh.py

@@ -660,3 +660,32 @@ def test_select_by_index_64(device):
         untouched_sphere.vertex.positions)
     assert sphere_custom.triangle.indices.allclose(
         untouched_sphere.triangle.indices)
+
+
+def check_compute_triangle_areas(device, int_t, float_t):
+    torus = o3d.t.geometry.TriangleMesh.create_torus(2, 1, 6, 3, float_t, int_t,
+                                                     device)
+
+    expected_areas = o3c.Tensor([
+        2.341874249399399, 1.1709371246996996, 1.299038105676658,
+        1.2990381056766576, 1.1709371246996996, 2.3418742493993996,
+        2.341874249399399, 1.1709371246996996, 1.299038105676658,
+        1.2990381056766573, 1.1709371246996996, 2.341874249399399,
+        2.341874249399399, 1.1709371246996998, 1.2990381056766582,
+        1.2990381056766576, 1.1709371246996993, 2.3418742493993996,
+        2.3418742493993987, 1.1709371246996996, 1.2990381056766578,
+        1.299038105676657, 1.1709371246996991, 2.3418742493993987,
+        2.3418742493993987, 1.1709371246996996, 1.299038105676658,
+        1.2990381056766573, 1.170937124699699, 2.341874249399399,
+        2.3418742493994, 1.1709371246997002, 1.299038105676659,
+        1.2990381056766582, 1.1709371246997, 2.3418742493994005
+    ], float_t, device)
+    assert torus.compute_triangle_areas().triangle.areas.allclose(
+        expected_areas)
+
+
+@pytest.mark.parametrize("device", list_devices())
+def test_compute_triangle_areas(device):
+    for int_t in [o3c.int32, o3c.int64]:
+        for float_t in [o3c.float32, o3c.float64]:
+            check_compute_triangle_areas(device, int_t, float_t)