VTK无法使用vtkClipClosedSurface构造适当的闭合曲面

Question

以下是我所做的粗略解释vtk:

1. 创建一个表面(一个最小的表面,与它不太相关,但几何结构很重要:陀螺仪有两个完全相互关闭的迷宫).
2. 用vtkClipClosedSurface切断迷宫之一,使我得到了没有打开的表面的对象.一个规则的表面看起来是这样的,有一个封闭的表面看起来是这样.

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

根据我的理解,vtkClipClosedSurface从表面法线知道在哪里关闭表面和不在哪里.问题是:我的结构的法线很好,它们都指向正确的方向.如果仔细观察结构,你会注意到下部基本上是顶部的反转,它逐渐变化,都在一个表面上.

我尝试修改我的结构,然后切割许多东西vtkSmoothPolyDataFilter,vtkCleanPolyData或者vtkPolyDataNormals.我甚至试图用边界表面提取vtkFeatureEdges,这导致了更糟糕的结果.甚至vtkFillHolesFilter没有产生任何可接受的结果.我的表面看起来完美无瑕,足以创造一个边界.

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

这是几何不能正确关闭曲面的几何示例.这次我使用的vtkFillHolesFilter结果是结构内部的表面,而它们只应占据te对象的边界.

如果你需要更详细的我的管道纲要,这里有:

1. 使用创建表面 mayavi.mlab.contour3d
2. 获得PolyData通过提取actor.mapper.input
3. 将格式转换tvtk为常规格式vtk
4. vtkClipClosedSurface 使用平面集合切除部分结构(当平面集合与结构边界相同时会发生错误)
5. 想象它

编辑:好的,这没有得到足够的重视,所以我构建了一个最小的,完整的,可验证的工作示例来重现行为:

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。我在生成一些低分辨率曼德尔灯泡时遇到了非常类似的问题。

你可能还想看看 pyvista，它是 vtk 的一个很好的 python 包装器。