Search code examples
pythonmeshcontouropen3d

Python Implementation for creating a triangular mesh from an array of closed loop planar contours

I'm a wee bit stuck.

I have a 3D point cloud (an array of (n,3) vertices), in which I am trying to generate a 3D triangular mesh from. So far I have had no luck.

The format my data comes in:

  • (x,y) values in regularly spaced (z) intervals. Think of the data as closed loop planar contours stored slice by slice in the z direction.
  • The vertices in my data must be absolute positions for the mesh triangles (i.e. I don't want them to be smoothed out such that the volume begins to change shape, but linear interpolation between the layers is fine).

Illustration:

Z=2. : ..x-------x...   <- Contour 2
Z=1.5: ...\......|...   <- Join the two contours into a mesh.
Z=1. : .....x----x...   <- Contour 1
Repeat for n slices, end up with an enclosed 3D triangular mesh.

Things I have tried:

  • Using Open3D:
    • The rolling ball (pivot) method can only get 75% of the mesh completed and leaves large areas incomplete (despite a range of ball sizes). It has particular problems at the top and bottom slices where there tends to be large gaps in the middle (i.e. a flat face).
    • The Poisson reconstruction method smooths out the volume too much and I no longer have an accurate representation of the volume. This occurs at all depths from 3-12.
  • CGAL:
    • I cannot get this to work for the life of me. SWIG is not very good, the implementation of CGAL using SWIG is also not very good.
    • There are two PyBind implementations of CGAL however they have not incorporated the 3D triangulation libraries from CGAL.
  • Explored other modules like PyMesh, TriMesh, TetGen, Scikit-Geometry, Shapely etc. etc. I may have missed the answer somewhere along the line.

Given that my data is a list of closed-loop planar contours, it seems as though there must be some simple solution to just "joining" adjacent slice contours into one big 3d mesh. Kind of like you would in blender. There are non-python solutions (like MeshLab) that may well solve these problems, but I require a python solution. Does anyone have any ideas? I've had a bit of a look into VTK and ITK but haven't found exactly what I'm looking for as of yet.

I'm also starting to consider that maybe I can interpolate intermediate contours between slices, and fill the contours on the top and bottom with vertices to make the data a bit more "pivot ball" method friendly.

Thank you in advance for any help, it is appreciated. If there is a good way of doing this that isn't coded yet, I promise to code it and make it available for people in my situation :)

Solution

  • There is a very neat paper titled "Technical Note: an algorithm and software for conversion of radiotherapy contour‐sequence data to ready‐to‐print 3D structures" in the Journal of Medical Physics that describes this problem quite nicely. No python packages are required, however it is more easily implemented with numpy. No need for any 3D packages.

    A useful excerpt is provided: ...

    1. The number of slices (2D contours) constituting the specified structure is determined.
    2. The number of points in each slice is determined.
    3. Cartesian coordinates of each of the points in each slice are extracted and stored within dedicated data structures...
    4. Numbers of points in each slice (curve) are re‐arranged in such a way, that the starting points (points with indices 0) are the closest points between the subsequent slices. Renumeration starts at point 0, slice 0 (slice with the lowest z coordinate).
    5. Orientation (i.e., the direction determined by the increasing indices of points with relation to the interior/exterior of the curve) of each curve is determined. If differences between slices are found, numbering of points in non‐matching curves (and thus, orientation) is reversed.
    6. The lateral surface of the considered structure is discretized. Points at the neighboring layers are arranged into threes, constituting triangular facets for the STL file. For each triangle the closest points with the subsequent indices from each layer are connected.
    7. Lower and upper base surfaces of the considered structure are discretized. The program iterates over every subsequent three points on the curve and checks if they belong to a convex part of the edge. If yes, they are connected into a facet, and the middle point is removed from further iterations.

    So basically it's a problem of aligning datasets in each slice to the nearest value of each slice. Then aligning the orientation of each contour. Then joining the points between two layers based on distance.

    The paper also provides code to do this (for a DICOM file), however I re-wrote it myself and it works a charm.

    I hope this helps others! Make sure you credit the author's in any work you do that uses this.