How to triangulate a surface from point cloud

T4mmi · December 9, 2019, 3:58pm

Hi all,

Quite new to vtk and meshing in general (got a raster processing background).

I understand the concept of mesh (and polydata) as a tuple of :

vertices/points (i.e. space coordinates): (x,y,z) in 3D
cells (i.e. list of vertices/points defining the geometric primitives): (pt1,pt2,pt3) for triangles

I did manage to convert Numpy array of coordinates to VTK Points, PointSet or even PolyData (yeah, working with python sorry)

import numpy as np
import vtk
from vtk.util.numpy_support import numpy_to_vtk, vtk_to_numpy
from vtk.numpy_interface import dataset_adapter as dsa

# lets create dummy {x,y,z} coordinates
pts = np.random.rand(512*3).reshape(-1,3)

# build points & polydata from numpy_to_vtk
points = vtk.vtkPoints()
points.SetData(numpy_to_vtk(pts))
poly = vtk.vtkPolyData()
poly.SetPoints(points)

# build points & polydata from dataset_adapter
point_set = dsa.PointSet(pts)
poly_data = dsa.PolyData(pts)

Dunno which way is the correct one … but the dataset_adapter way comes with more problems:

print(f"There are {point_set.GetNumberOfPoints()} points.")
AttributeError: 'numpy.ndarray' object has no attribute 'GetNumberOfPoints'

Anyway I want to just build the cells list from these vertices

what is the correct and most efficient way to do this ?

lassoan · December 9, 2019, 5:39pm

We usually create a VTK array that has the appropriate size, get read-write access to it using vtk_to_numpy, then copy values from numpy array (element-wise, using a[:]=b syntax). This way, arrays created in VTK are owned by VTK and arrays created by numpy are owned by numpy objects, so we don’t need to worry about reallocation or deleting these arrays.

You can get read-write access as numpy array to point coordinates like this:

github.com

Slicer/Slicer/blob/0abd1fa219628f0dfd5166e430afeb1598b6fcf8/Base/Python/slicer/util.py#L976-L988


def arrayFromModelPoints(modelNode):
  """Return point positions of a model node as numpy array.
  Point coordinates can be modified by modifying the numpy array.
  After all modifications has been completed, call :py:meth:`arrayFromModelPointsModified`.


  .. warning:: Important: memory area of the returned array is managed by VTK,
    therefore values in the array may be changed, but the array must not be reallocated.
    See :py:meth:`arrayFromVolume` for details.
  """
  import vtk.util.numpy_support
  pointData = modelNode.GetPolyData().GetPoints().GetData()
  narray = vtk.util.numpy_support.vtk_to_numpy(pointData)
  return narray

You can access cells similarly.

Surface reconstruction from point cloud is a hard problem. VTK has a few filters that can do this, but of course each has its limitations. If your source data contains cells (not just a point cloud) then import those cells into VTK, too.

T4mmi · December 10, 2019, 8:20am

Hi @lassoan,
Thanks for your answer. It brings more questions actually …

I understand the concerns about sharing memory between VTK <> Numpy,
but when aware of the pitfalls it avoids HUGE data copies in both ways (both RAM and time consuming basically for nothing else than safety).
Anyway, I don’t have any other choice since VTK adapters does not seems to work as expected…

Surface reconstruction from point cloud is a hard problem. VTK has a few filters that can do this, but of course each has its limitations.

Yeah ok, but where can we find the doc about these methods ?!
And in my particular case, I dont wanna extract the surface, I just wanna mesh/triangulate the points (i.e. create cells from a uniform grid array of z coordinates)

If your source data contains cells (not just a point cloud) then import those cells into VTK, too.

Actually thats why I’m trying to triangulate these points, because I don’t have the cells

Paulo_Carvalho · December 10, 2019, 2:24pm

Hello, Thomas. Your problem is called “mesh generation” (a.k.a tesselation) which is an entire research field in itself. Please, take a look at this: https://stackoverflow.com/questions/4882993/mesh-generation-from-points-with-x-y-and-z-coordinates for suggestions of libraries specialized on mesh generation.

VTK of course has such an algorithm in a class called vtkDelaunay3D but, as stated in the topic, Delaunay 3D may yield poor results depending, of course, on your input data set (point cloud).

lassoan · December 10, 2019, 9:05pm

There are a couple of more, see for example:

You may use vtkDepthImageToPointCloud filter to create point cloud from your depth image.

They work as expected. You just have to be careful not to delete or modify memory layout of a buffer that is used in both VTK and numpy (I don’t think they can notify each other about the change). In general, it is not worth spending a lot of time with maintaining a shared buffer and easier to just deep-copy the data instead. However, if performance is critical then using shared buffers may be the better choice, despite the additional complexity.

banesullivan · December 12, 2019, 5:39am

@T4mmi, you may benefit from checking out PyVista’s interface to VTK to make working with your VTK dataset a bit easier. Here’s an example with the numpy array you shared:

import pyvista as pv
import numpy as np

pts = np.random.rand(512*3).reshape(-1,3)

# Make vtkPolyData of the points array
point_cloud = pv.PolyData(pts)
point_cloud.plot(render_points_as_spheres=True, point_size=10)

# runs the delaunay 3D algorithm
mesh = point_cloud.delaunay_3d(alpha=0.25)

# Make a wireframe with some data
wires = mesh.compute_cell_sizes(length=False, area=False, volume=True).wireframe()
# Plot it
wires.plot(line_width=2)