Find what is material and what is air in 2D cross sections from 3d .stl files using python

Question

I would like to create 2D cross-sections from 3D files, without losing the information of what is material and what is air.

In the end, I would like to obtain a dictionary which contains the outer-most points that make up the material and air inclusions (can be multiple), i.e.

"material" : [[x1,y1],[x2,y2]...]
"air_inclusions": [[[x11,y11],[x12,y12],...],[[x21,y21],[x22,y22],...],[[x31,y31],[x32,y32],...]]

Here is an example of how I am trying to do it:

I have the following .stl file, which you can download here https://filebin.net/c9o0zy4bnv8dvuew

Using the amazing python package trimesh, I can import the .stl file

import trimesh
import numpy as np   

mesh = trimesh.load_mesh(r"PATH_TO_FILE")
# give it a color
mesh.visual.face_colors = [100, 100, 100, 255]
# and show it
mesh.show(viewer='gl')

Creating the 2D slide

# I can create a 2D slice of the geometry at origin [0,0,5] and slice-plane with normal direction [0,0,1] 
slice = mesh.section(plane_origin=[0,0,5],
                     plane_normal=[0,0,1])
slice.show(viewer='gl')

Extracting the vertices

# take 2D slice (before was still 3D) 
slice_2D, to_3D = slice.to_planar()
# get vertices
vertices =  np.asanyarray(slice_2D.vertices)

# plot 
import matplotlib.pyplot as plt
x,y = vertices.T
plt.scatter(x,y,s=0.4)
plt.show()

My approach to retrieve the information about what is material and what is air

My assumption

The outer most points define the border for the material. All points within define the border of air inclusions.

I get the outer most point -> convex hull

from scipy.spatial import ConvexHull

# compute the hull
hull = ConvexHull(vertices)
# plot
plt.plot(vertices[:,0], vertices[:,1], 'o')
for simplex in hull.simplices:
  plt.plot(vertices[simplex, 0], vertices[simplex, 1], 'k-')

To know all the points inside of the hull, I use this answer What's an efficient way to find if a point lies in the convex hull of a point cloud?

# Source: https://stackoverflow.com/questions/16750618/whats-an-efficient-way-to-find-if-a-point-lies-in-the-convex-hull-of-a-point-cl
def in_hull(p, hull):
    """
    Test if points in `p` are in `hull`

    `p` should be a `NxK` coordinates of `N` points in `K` dimensions
    `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the 
    coordinates of `M` points in `K`dimensions for which Delaunay triangulation
    will be computed
    """
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    return hull.find_simplex(p)>=0

I gather the remaining points

# Remaining points

remaining = []
for i,in_hull in enumerate(in_hull(vertices,hull.simplices)):
    if in_hull:
        remaining.append(vertices[i])

The issues

1. The remaining points are only two points, but it should be many more, as can be seen in the plot above. Why is that and how can I fix it?

   [TrackedArray([21.60581633, 8.99397324]), TrackedArray([12.95590211, 23.97608075])]

2. Do you have any idea of how I could find all air inclusions if there are more than one? i.e.

you can find the file here: https://filebin.net/6blzvrrwhanv0jib

Answer 1

Thanks to shapely polygons, on which trimesh is built you don't need to go over the vertices. You can use some of the inbuilt functions. After creating your 2D path you are almost there. The path has two attributes: polygons_full and polygons_closed. The first one is the outmost polygon without the inner ones and the second one are all polygons of your path.
You can simply do:

slice_2D, to_3D = slice.to_planar()
# create a new figure to which you can attach patches
fig = plt.figure(1)
ax = fig.add_subplot(111)
# Here you get your outmost polygons and add them as patches to your plot
for p in slice_2D.polygons_full:
  ax.add_patch(PolygonPatch(p))
# this is needed due to the differences of polygons_full and polygons_closed to check, if the polygon is one of the outer polygons
outer_polys = [x.exterior for x in slice_2D.polygons_full]
# iterate over all polygons and check, whether they are one of the outmost polygons. If not plot it (the outmost ones have already been added as patches).
for p in (slice_2D.polygons_closed):
  if p.exterior not in outer_polys:
    plt.plot(*(p.exterior.xy), 'r')
# show the plot
plt.show()

Plot:

Or you can shorten that down using the interior property of polygons:

slice_2D, to_3D = slice.to_planar()
for p in slice_2D.polygons_full:
  plt.plot(*(p.exterior.xy),'k')
  for r in p.interiors:
    plt.plot(*zip(*r.coords), 'b')

Interiors are the "holes" of a filled polygon, so this should be exactly what you want. They are a LinearRing so you can't directly use Polygon properties.