Найдите, что такое материал и что такое воздух в 2D сечениях из файлов 3D .stl, используя python

Question

Я бы хотел создать 2D-сечения из 3D-файлов, не теряя информацию о том, что является материалом, а что - воздухом.

В конце я хотел бы получить словарь, который содержит самые внешние точки, которые составляют материальные и воздушные включения (могут быть множественными), то есть

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

Вот пример того, как я пытаюсь это сделать:

У меня есть следующий файл .stl, который вы можете скачать здесь https://filebin.net/c9o0zy4bnv8dvuew

Использование удивительного python пакета trime sh, я могу импортировать файл .stl

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')

Создание 2D слайда

# 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')

Извлечение вершин

# 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()

Мой подход для получения информации о что такое материал и что такое воздух

Мое предположение

Внешние точки определяют границы материала. Все точки внутри определяют границы воздушных включений.

Я получаю самую крайнюю точку -> выпуклый корпус

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-')

Чтобы узнать все точки внутри корпуса, я использую этот ответ Какой эффективный способ определить, лежит ли точка в выпуклой оболочке облака точек?

# 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

Я набираю оставшиеся баллы

# Remaining points

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

Проблемы

1. Остальные баллы - это только два балла, но их должно быть гораздо больше, как видно на графике выше. Почему это так и как я могу это исправить?

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

2. У вас есть представление о том, как я могу найти все воздушные включения, если их несколько? то есть

Вы можете найти файл здесь: https://filebin.net/6blzvrrwhanv0jib

Answer 1

Благодаря стройным многоугольникам, на которых построен трим sh, вам не нужно go над вершинами. Вы можете использовать некоторые из встроенных функций. После создания 2D-пути вы почти у цели. Путь имеет два атрибута: polygons_full и polygons_closed. Первый - это крайний многоугольник без внутренних, а второй - многоугольники вашего пути.
Вы можете просто сделать:

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()

Сюжет:

Или вы можете сократить это значение, используя свойство interior полигонов:

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')

Интерьеры - это «дыры» заполненного многоугольника, поэтому это должно быть именно то, что вы хотите. Они LinearRing, поэтому вы не можете напрямую использовать свойства Polygon.