Add GetSurfaceArea() for t::TriangleMesh (#6248)

cpp/open3d/core/linalg/kernel/Matrix.h

@@ -72,6 +72,18 @@ OPEN3D_HOST_DEVICE OPEN3D_FORCE_INLINE void cross_3x1(
             A_3x1_input[0] * B_3x1_input[1] - A_3x1_input[1] * B_3x1_input[0];
 }
 
+template <typename scalar_t>
+OPEN3D_HOST_DEVICE OPEN3D_FORCE_INLINE scalar_t
+cross_mag_3x1(const scalar_t* A_3x1_input, const scalar_t* B_3x1_input) {
+    scalar_t temp_0 =
+            A_3x1_input[1] * B_3x1_input[2] - A_3x1_input[2] * B_3x1_input[1];
+    scalar_t temp_1 =
+            A_3x1_input[2] * B_3x1_input[0] - A_3x1_input[0] * B_3x1_input[2];
+    scalar_t temp_2 =
+            A_3x1_input[0] * B_3x1_input[1] - A_3x1_input[1] * B_3x1_input[0];
+    return sqrt(temp_0 * temp_0 + temp_1 * temp_1 + temp_2 * temp_2);
+}
+
 template <typename scalar_t>
 OPEN3D_HOST_DEVICE OPEN3D_FORCE_INLINE scalar_t
 dot_3x1(const scalar_t* A_3x1_input, const scalar_t* B_3x1_input) {

cpp/open3d/t/geometry/TriangleMesh.cpp

@@ -283,6 +283,39 @@ TriangleMesh &TriangleMesh::ComputeVertexNormals(bool normalized) {
     return *this;
 }
 
+double TriangleMesh::GetSurfaceArea() const {
+    double surface_area = 0;
+    if (IsEmpty()) {
+        utility::LogWarning("TriangleMesh is empty.");
+        return surface_area;
+    }
+
+    if (!HasTriangleIndices()) {
+        utility::LogWarning("TriangleMesh has no triangle indices.");
+        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");
+    }
+
+    surface_area = triangle_areas.Sum({0}).To(core::Float64).Item<double>();
+
+    return surface_area;
+}
+
 geometry::TriangleMesh TriangleMesh::FromLegacy(
         const open3d::geometry::TriangleMesh &mesh_legacy,
         core::Dtype float_dtype,

cpp/open3d/t/geometry/TriangleMesh.h

@@ -665,6 +665,10 @@ class TriangleMesh : public Geometry, public DrawableGeometry {
     /// rendering.
     TriangleMesh &ComputeVertexNormals(bool normalized = true);
 
+    /// \brief Function that computes the surface area of the mesh, i.e. the sum
+    /// of the individual triangle surfaces.
+    double GetSurfaceArea() const;
+
     /// \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/open3d/t/geometry/kernel/TriangleMesh.h

@@ -25,6 +25,10 @@ void ComputeVertexNormalsCPU(const core::Tensor& triangles,
                              const core::Tensor& triangle_normals,
                              core::Tensor& vertex_normals);
 
+void ComputeTriangleAreasCPU(const core::Tensor& vertices,
+                             const core::Tensor& triangles,
+                             core::Tensor& triangle_areas);
+
 #ifdef BUILD_CUDA_MODULE
 void NormalizeNormalsCUDA(core::Tensor& normals);
 
@@ -35,6 +39,10 @@ void ComputeTriangleNormalsCUDA(const core::Tensor& vertices,
 void ComputeVertexNormalsCUDA(const core::Tensor& triangles,
                               const core::Tensor& triangle_normals,
                               core::Tensor& vertex_normals);
+
+void ComputeTriangleAreasCUDA(const core::Tensor& vertices,
+                              const core::Tensor& triangles,
+                              core::Tensor& triangle_areas);
 #endif
 
 }  // namespace trianglemesh

cpp/open3d/t/geometry/kernel/TriangleMeshImpl.h

@@ -107,6 +107,52 @@ void ComputeTriangleNormalsCPU
     });
 }
 
+#if defined(__CUDACC__)
+void ComputeTriangleAreasCUDA
+#else
+void ComputeTriangleAreasCPU
+#endif
+        (const core::Tensor& vertices,
+         const core::Tensor& triangles,
+         core::Tensor& triangle_areas) {
+    const int64_t n = triangle_areas.GetLength();
+    const core::Dtype dtype = triangle_areas.GetDtype();
+    const core::Tensor triangles_d = triangles.To(core::Int64);
+
+    DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(dtype, [&]() {
+        scalar_t* area_ptr = triangle_areas.GetDataPtr<scalar_t>();
+        const int64_t* triangle_ptr = triangles_d.GetDataPtr<int64_t>();
+        const scalar_t* vertex_ptr = vertices.GetDataPtr<scalar_t>();
+
+        core::ParallelFor(triangle_areas.GetDevice(), n,
+                          [=] OPEN3D_DEVICE(int64_t workload_idx) {
+                              int64_t idx = 3 * workload_idx;
+
+                              int64_t triangle_id1 = triangle_ptr[idx];
+                              int64_t triangle_id2 = triangle_ptr[idx + 1];
+                              int64_t triangle_id3 = triangle_ptr[idx + 2];
+
+                              scalar_t v01[3], v02[3];
+                              v01[0] = vertex_ptr[3 * triangle_id2] -
+                                       vertex_ptr[3 * triangle_id1];
+                              v01[1] = vertex_ptr[3 * triangle_id2 + 1] -
+                                       vertex_ptr[3 * triangle_id1 + 1];
+                              v01[2] = vertex_ptr[3 * triangle_id2 + 2] -
+                                       vertex_ptr[3 * triangle_id1 + 2];
+                              v02[0] = vertex_ptr[3 * triangle_id3] -
+                                       vertex_ptr[3 * triangle_id1];
+                              v02[1] = vertex_ptr[3 * triangle_id3 + 1] -
+                                       vertex_ptr[3 * triangle_id1 + 1];
+                              v02[2] = vertex_ptr[3 * triangle_id3 + 2] -
+                                       vertex_ptr[3 * triangle_id1 + 2];
+
+                              area_ptr[workload_idx] =
+                                      0.5 * core::linalg::kernel::cross_mag_3x1(
+                                                    v01, v02);
+                          });
+    });
+}
+
 }  // namespace trianglemesh
 }  // namespace kernel
 }  // namespace geometry

cpp/pybind/t/geometry/trianglemesh.cpp

@@ -191,6 +191,20 @@ The attributes of the triangle mesh have different levels::
                       "before rendering.",
                       "normalized"_a = true);
 
+    triangle_mesh.def(
+            "get_surface_area", &TriangleMesh::GetSurfaceArea,
+            R"(Computes the surface area of the mesh, i.e., the sum of the individual triangle surfaces.
+
+Example:
+    This computes the surface area of the Stanford Bunny::
+        bunny = o3d.data.BunnyMesh()
+        mesh = o3d.t.io.read_triangle_mesh(bunny.path)
+        print('The surface area is', mesh.get_surface_area())
+
+Returns:
+    A scalar describing the surface area of the mesh.
+)");
+
     triangle_mesh.def(
             "compute_convex_hull", &TriangleMesh::ComputeConvexHull,
             "joggle_inputs"_a = false,

python/test/t/geometry/test_trianglemesh.py

@@ -419,6 +419,27 @@ def test_pickle(device):
                                 mesh.triangle.indices.cpu().numpy())
 
 
+@pytest.mark.parametrize("device", list_devices())
+def test_get_surface_area(device):
+    # Test with custom parameters.
+    cube = o3d.t.geometry.TriangleMesh.create_box(float_dtype=o3c.float64,
+                                                  int_dtype=o3c.int32,
+                                                  device=device)
+    np.testing.assert_equal(cube.get_surface_area(), 6)
+
+    empty = o3d.t.geometry.TriangleMesh(device=device)
+    empty.get_surface_area()
+    np.testing.assert_equal(empty.get_surface_area(), 0)
+
+    # test noncontiguous
+    sphere = o3d.t.geometry.TriangleMesh.create_sphere(device=device)
+    area1 = sphere.get_surface_area()
+    sphere.vertex.positions = sphere.vertex.positions.T().contiguous().T()
+    sphere.triangle.indices = sphere.triangle.indices.T().contiguous().T()
+    area2 = sphere.get_surface_area()
+    np.testing.assert_almost_equal(area1, area2)
+
+
 @pytest.mark.parametrize("device", list_devices())
 def test_select_faces_by_mask_32(device):
     sphere_custom = o3d.t.geometry.TriangleMesh.create_sphere(