VTK based implementation of extract trianglemesh

cpp/open3d/t/geometry/TriangleMesh.h

@@ -594,6 +594,20 @@ class TriangleMesh : public Geometry, public DrawableGeometry {
             core::Dtype int_dtype = core::Int64,
             const core::Device &device = core::Device("CPU:0"));
 
+    /// Create a mesh from a 3D scalar field (volume) by computing the
+    /// isosurface.  This method uses the Flying Edges dual contouring method
+    /// that computes the isosurface similar to Marching Cubes.  The center of
+    /// the first voxel of the volume is at the origin (0,0,0).  The center of
+    /// the voxel at index [z,y,x] will be at (x,y,z).
+    /// \param volume 3D tensor with the volume.
+    /// \param contour_values A list of contour values at which isosurfaces will
+    /// be generated. The default value is 0.
+    /// \param device The device for the returned mesh.
+    static TriangleMesh CreateIsosurfaces(
+            const core::Tensor &volume,
+            const std::vector<double> contour_values = {0.0},
+            const core::Device &device = core::Device("CPU:0"));
+
 public:
     /// Clear all data in the trianglemesh.
     TriangleMesh &Clear() override {

cpp/open3d/t/geometry/TriangleMeshFactory.cpp

@@ -5,6 +5,7 @@
 // SPDX-License-Identifier: MIT
 // ----------------------------------------------------------------------------
 
+#include <vtkFlyingEdges3D.h>
 #include <vtkLinearExtrusionFilter.h>
 #include <vtkNew.h>
 #include <vtkTextSource.h>
@@ -285,6 +286,28 @@ TriangleMesh TriangleMesh::CreateText(const std::string &text,
     return tmesh;
 }
 
+TriangleMesh TriangleMesh::CreateIsosurfaces(
+        const core::Tensor &volume,
+        const std::vector<double> contour_values,
+        const core::Device &device) {
+    using namespace vtkutils;
+    core::AssertTensorShape(volume, {core::None, core::None, core::None});
+    core::AssertTensorDtypes(volume, {core::Float32, core::Float64});
+
+    auto image_data = vtkutils::CreateVtkImageDataFromTensor(
+            const_cast<core::Tensor &>(volume));
+    vtkNew<vtkFlyingEdges3D> method;
+    method->SetNumberOfContours(contour_values.size());
+    for (int i = 0; i < int(contour_values.size()); ++i) {
+        method->SetValue(i, contour_values[i]);
+    }
+    method->SetInputData(image_data);
+    method->Update();
+    auto polydata = method->GetOutput();
+    auto tmesh = CreateTriangleMeshFromVtkPolyData(polydata);
+    return tmesh.To(device);
+}
+
 }  // namespace geometry
 }  // namespace t
 }  // namespace open3d

cpp/open3d/t/geometry/VtkUtils.cpp

@@ -12,6 +12,7 @@
 #include <vtkCellData.h>
 #include <vtkDoubleArray.h>
 #include <vtkFloatArray.h>
+#include <vtkImageData.h>
 #include <vtkLinearExtrusionFilter.h>
 #include <vtkPointData.h>
 #include <vtkPoints.h>
@@ -182,6 +183,34 @@ static vtkSmartPointer<vtkPoints> CreateVtkPointsFromTensor(
     return pts;
 }
 
+OPEN3D_LOCAL vtkSmartPointer<vtkImageData> CreateVtkImageDataFromTensor(
+        core::Tensor& tensor, bool copy) {
+    core::AssertTensorDtypes(tensor,
+                             {core::UInt8, core::Float32, core::Float64});
+    if (tensor.NumDims() != 2 && tensor.NumDims() != 3) {
+        utility::LogError(
+                "Cannot convert Tensor to vtkImageData. The number of "
+                "dimensions must be 2 or 3 but is {}",
+                tensor.NumDims());
+    }
+
+    // Create a flat tensor that can be converted to a vtkDataArray
+    auto tensor_flat = tensor.Reshape({tensor.NumElements(), 1});
+    if (tensor.GetDataPtr() != tensor_flat.GetDataPtr()) {
+        copy = true;
+    }
+    auto data_array = CreateVtkDataArrayFromTensor(tensor_flat, copy);
+
+    vtkSmartPointer<vtkImageData> im = vtkSmartPointer<vtkImageData>::New();
+    im->GetPointData()->SetScalars(data_array);
+    std::array<int, 3> size{1, 1, 1};
+    for (int i = 0; i < tensor.NumDims(); ++i) {
+        size[i] = tensor.GetShape(tensor.NumDims() - i - 1);
+    }
+    im->SetDimensions(size.data());
+    return im;
+}
+
 namespace {
 // Helper for creating the offset array from Common/DataModel/vtkCellArray.cxx
 struct GenerateOffsetsImpl {
@@ -213,7 +242,9 @@ static vtkSmartPointer<vtkCellArray> CreateVtkCellArrayFromTensor(
     const int cell_size = tensor.GetShape()[1];
 
     auto tensor_flat = tensor.Reshape({tensor.NumElements(), 1}).Contiguous();
-    copy = copy && tensor.GetDataPtr() == tensor_flat.GetDataPtr();
+    if (tensor.GetDataPtr() != tensor_flat.GetDataPtr()) {
+        copy = true;
+    }
     auto connectivity = CreateVtkDataArrayFromTensor(tensor_flat, copy);
 
     // vtk nightly build (9.1.20220520) has a function cells->SetData(cell_size,

cpp/open3d/t/geometry/VtkUtils.h

@@ -5,6 +5,7 @@
 // SPDX-License-Identifier: MIT
 // ----------------------------------------------------------------------------
 
+#include <vtkImageData.h>
 #include <vtkPolyData.h>
 #include <vtkSmartPointer.h>
 
@@ -22,6 +23,14 @@ namespace vtkutils {
 /// Logs an error if no conversion exists.
 int DtypeToVtkType(const core::Dtype& dtype);
 
+/// Creates a vtkImageData object from a Tensor.
+/// The returned object may directly use the memory of the tensor and the tensor
+/// must be kept alive until the returned vtkImageData is deleted.
+/// \param tensor The source tensor.
+/// \param copy If true always create a copy of the data.
+vtkSmartPointer<vtkImageData> CreateVtkImageDataFromTensor(core::Tensor& tensor,
+                                                           bool copy = false);
+
 /// Creates a vtkPolyData object from a point cloud or triangle mesh.
 /// The returned vtkPolyData object may directly use the memory of the tensors
 /// stored inside the Geometry object. Therefore, the Geometry object must be

cpp/pybind/t/geometry/trianglemesh.cpp

@@ -519,6 +519,48 @@ This example shows how to create a hemisphere from a sphere::
 
         mesh = o3d.t.geometry.TriangleMesh.create_text('Open3D', depth=1)
         o3d.visualization.draw([{'name': 'text', 'geometry': mesh}])
+)");
+
+    triangle_mesh.def_static(
+            "create_isosurfaces",
+            // Accept anything for contour_values that pybind can convert to
+            // std::list. This also avoids o3d.utility.DoubleVector.
+            [](const core::Tensor& volume, std::list<double> contour_values,
+               const core::Device& device) {
+                std::vector<double> cv(contour_values.begin(),
+                                       contour_values.end());
+                return TriangleMesh::CreateIsosurfaces(volume, cv, device);
+            },
+            "volume"_a, "contour_values"_a = std::list<double>{0.0},
+            "device"_a = core::Device("CPU:0"),
+            R"(Create a mesh from a 3D scalar field (volume) by computing the
+isosurface.
+
+This method uses the Flying Edges dual contouring method that computes the
+isosurface similar to Marching Cubes. The center of the first voxel of the
+volume is at the origin (0,0,0). The center of the voxel at index [z,y,x]
+will be at (x,y,z).
+
+Args:
+    volume (open3d.core.Tensor): 3D tensor with the volume.
+    contour_values (list): A list of contour values at which isosurfaces will
+        be generated. The default value is 0.
+    device (o3d.core.Device): The device for the returned mesh.
+
+Returns:
+    A TriangleMesh with the extracted isosurfaces.
+
+
+This example shows how to create a sphere from a volume::
+
+    import open3d as o3d
+    import numpy as np
+
+    coords = np.stack(np.meshgrid(*3*[np.linspace(-1,1,num=64)], indexing='ij'), axis=-1)
+    vol = np.linalg.norm(coords, axis=-1) - 0.5
+    mesh = o3d.t.geometry.TriangleMesh.create_isosurfaces(vol)
+    o3d.visualization.draw(mesh)
+
 )");
 
     triangle_mesh.def(

python/test/t/geometry/test_trianglemesh.py

@@ -287,6 +287,19 @@ def test_create_text():
     assert mesh.triangle.indices.shape == (936, 3)
 
 
+def test_create_isosurfaces():
+    """Create signed distance field for sphere of radius 0.5 and extract sphere
+    from it.
+    """
+    coords = np.stack(np.meshgrid(*3 * [np.linspace(-1, 1, num=64)],
+                                  indexing='ij'),
+                      axis=-1)
+    vol = np.linalg.norm(coords, axis=-1) - 0.5
+    mesh = o3d.t.geometry.TriangleMesh.create_isosurfaces(vol)
+    assert mesh.vertex.positions.shape[0] == 4728
+    assert mesh.triangle.indices.shape[0] == 9452
+
+
 def test_simplify_quadric_decimation():
     cube = o3d.t.geometry.TriangleMesh.from_legacy(
         o3d.geometry.TriangleMesh.create_box().subdivide_midpoint(3))