python - VTK 不能用 vtkClipClosedSurface 构造一个合适的闭合曲面

Question

这是我在 vtk 中所做的粗略解释:

创建一个表面(一个最小的表面，它是什么不太相关，但几何形状很重要:gyroid 有两个完全相互隔离的迷宫)。

使用 vtkClipClosedSurface关闭其中一个迷宫，以便我得到一个不再有开放表面的物体。 regular surface looks like this ，它看起来是一个封闭的表面 like this .

这是我的问题:对于我的结构的更复杂版本，我得到了这个:
你能看到左上角它如何正常工作，而在右下角附近它停止创建曲面吗？有时我也会在最后一部分得到非常奇怪的三角形。

据我了解 vtkClipClosedSurface从表面法线知道在哪里关闭表面，在哪里不关闭。问题是:我的结构的法线很好，它们都指向正确的方向。如果您仔细观察结构，您会注意到下部基本上是顶部的倒置，逐渐变化，全部在一个表面上。

我试图在切割之前修改我的结构，比如 vtkSmoothPolyDataFilter , vtkCleanPolyData或 vtkPolyDataNormals .我什至尝试用 vtkFeatureEdges 提取边界面，这导致了更糟糕的结果。偶 vtkFillHolesFilter没有产生任何可接受的结果。我的表面看起来完美无瑕，很容易创建一个边界。

我不知道还有什么可尝试的。其他结构也会发生这种情况。用 CAD 工具修复它是不可能的，因为它应该是开箱即用的。请帮我!

这是未正确闭合曲面的几何体的另一个示例。这次我用了 vtkFillHolesFilter这导致结构内部的表面，而它们应该只占据 te 对象的边界。


如果您需要更详细地了解我的管道，请看这里:

使用 mayavi.mlab.contour3d 创建曲面

获取 PolyData通过提取 actor.mapper.input

从 tvtk 转换格式到正规 vtk

vtkClipClosedSurface带有切掉部分结构的平面集合(平面集合与结构边界相同时会出现错误)

想象一下

编辑:好的，这个没有得到足够的关注，所以我构建了一个 minimal, complete and verifiable working example重现行为:

import numpy as np
import vtk  # VTK version 7.0
from mayavi import mlab  # mayavi version 4.4.4
from mayavi.api import Engine, OffScreenEngine
from tvtk.api import tvtk


def schwarz_D(x, y, z, linear_term=0):
        """This is the function for the Schwarz Diamond level surface."""
        return (np.sin(x) * np.sin(y) * np.sin(z) + np.sin(x) * np.cos(y) * np.cos(z) +
                np.cos(x) * np.sin(y) * np.cos(z) + np.cos(x) * np.cos(y) * np.sin(z)) - linear_term * z


def plane_collection(xn, x, yn, y, zn, z):
        """Defines the 6 planes for cutting rectangular objects to the right size."""
        plane1 = vtk.vtkPlane()
        plane1.SetOrigin(x, 0, 0)
        plane1.SetNormal(-1, 0, 0)
        plane2 = vtk.vtkPlane()
        plane2.SetOrigin(0, y, 0)
        plane2.SetNormal(0, -1, 0)
        plane3 = vtk.vtkPlane()
        plane3.SetOrigin(0, 0, z)
        plane3.SetNormal(0, 0, -1)
        plane4 = vtk.vtkPlane()
        plane4.SetOrigin(xn, 0, 0)
        plane4.SetNormal(1, 0, 0)
        plane5 = vtk.vtkPlane()
        plane5.SetOrigin(0, yn, 0)
        plane5.SetNormal(0, 1, 0)
        plane6 = vtk.vtkPlane()
        plane6.SetOrigin(0, 0, zn)
        plane6.SetNormal(0, 0, 1)

        plane_list = [plane4, plane1, plane5, plane2, plane6, plane3]
        planes = vtk.vtkPlaneCollection()
        for item in plane_list:
            planes.AddItem(item)
        return planes

[nx, ny, nz] = [2, 2, 8]  # amount of unit cells
cell_size = 1
gradient_value = 0.04  # only values below 0.1 produce the desired geometry; this term is essential
x, y, z = np.mgrid[-cell_size*(nx + 1)/2:cell_size*(nx + 1)/2:100j,
                   -cell_size*(ny + 1)/2:cell_size*(ny + 1)/2:100j,
                   -cell_size*(nz + 1)/2:cell_size*(nz + 1)/2:100*2j] * np.pi / (cell_size/2)

# engine = Engine()
engine = OffScreenEngine()  # do not start mayavi GUI
engine.start()
fig = mlab.figure(figure=None, engine=engine)
contour3d = mlab.contour3d(x, y, z, schwarz_D(x, y, z, gradient_value), figure=fig)

scene = engine.scenes[0]
actor = contour3d.actor.actors[0]
iso_surface = scene.children[0].children[0].children[0]
iso_surface.contour.minimum_contour = 0
iso_surface.contour.number_of_contours = 1
iso_surface.compute_normals = False

iso_surface.contour.auto_update_range = False

mlab.draw(fig)
# mlab.show()  # enable if you want to see the mayavi GUI

polydata = tvtk.to_vtk(actor.mapper.input)  # convert tvtkPolyData to vtkPolyData

# Move object to the coordinate center to make clipping easier later on.
center_coords = np.array(polydata.GetCenter())
center = vtk.vtkTransform()
center.Translate(-center_coords[0], -center_coords[1], -center_coords[2])
centerFilter = vtk.vtkTransformPolyDataFilter()
centerFilter.SetTransform(center)
centerFilter.SetInputData(polydata)
centerFilter.Update()

# Reverse normals in order to receive a closed surface after clipping
reverse = vtk.vtkReverseSense()
reverse.SetInputConnection(centerFilter.GetOutputPort())
reverse.ReverseNormalsOn()
reverse.ReverseCellsOn()
reverse.Update()

bounds = np.asarray(reverse.GetOutput().GetBounds())
clip = vtk.vtkClipClosedSurface()
clip.SetInputConnection(reverse.GetOutputPort())
clip.SetTolerance(10e-3)
# clip.TriangulationErrorDisplayOn()  # enable to see errors for not watertight surfaces
clip.SetClippingPlanes(plane_collection(bounds[0] + cell_size/2, bounds[1] - cell_size/2,
                                        bounds[2] + cell_size/2, bounds[3] - cell_size/2,
                                        bounds[4] + cell_size/2, bounds[5] - cell_size/2))
clip.Update()

# Render the result
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(clip.GetOutputPort())
actor = vtk.vtkActor()
actor.SetMapper(mapper)
renderer = vtk.vtkRenderer()
renderWindow = vtk.vtkRenderWindow()
renderWindow.AddRenderer(renderer)
renderWindowInteractor = vtk.vtkRenderWindowInteractor()
renderWindowInteractor.SetRenderWindow(renderWindow)
renderer.AddActor(actor)
renderWindow.Render()
renderWindowInteractor.Start()

这真的很短，我尽可能多地剥离。问题仍然存在，我想不出解决方案。

Answer 1

尝试使用 pymeshfix .我在生成的一些低分辨率 mandelbulb 中遇到了非常相似的问题。
您可能还想查看 pyvista，它是 vtk 的一个很好的 Python 包装器。