☰
VTK实现Delaunay三角剖分与等距切割:篱笆图生成全攻略
2026/10/17 10:01:05 网站建设 项目流程

篱笆图这个词,乍一听像工程制图里的东西,但搞过地质可视化或者三维表面重建的人,应该都遇到过类似的诉求:有一堆散落的采样点,想看看沿着某个方向、在固定的间距上,模型的轮廓到底是怎么变化的。传统做法是把数据插值成规则网格,然后再一个面一个面地切,麻烦不说,碰到点分布不均匀的时候,插出来的面常常“名不副实”。我自己这些年处理这种需求,最顺手的方案就是VTK里Delaunay2D做三角剖分,再用一组等距平面去切这个三角网,把切出来的交线全部叠到一起。这就是所谓的篱笆图,英文叫fence diagram,本质上就是等距切割。

这篇文章就把这条完整的路线拆开讲清楚:散点怎么变成三角网,等距切割的平面怎么选,VTK的Cutter用起来有哪些隐藏的细节,以及实际操作时最容易踩的坑。无论你是刚接触VTK,还是已经在项目里用了一段时间想优化方案,应该都能从这里拿到一套可以直接复现的思路。当然,老规矩,所有东西都是我实际跑过的用法,不是说明书式的罗列。

1. 为什么先做Delaunay三角剖分,再谈等距切割

篱笆图的核心不是“绘图”,而是“切面”。要切面,前提是你得先有一个连续的、可求交的几何表面。而手头的原始数据往往只是一堆零散的点——来自钻孔、测点或者扫描仪——这些点本身根本没有拓扑关系,切不了。所以第一步,必须让这些点“长”成一个完整的表面。

1.1 不规则的散点,用什么办法生成表面

对规则网格数据,最省事的做法是直接生成结构网格,四边形或者六面体都行,因为它自带行列关系,天然知道哪些点相邻。但散点数据不一样,点与点之间没有隐含的邻居关系,必须由算法去判断“哪些点应该连在一起”。

这里最常见的两个派别,一是反距离加权或者克里金插值到规则网格,二是直接做三角剖分。前者虽然好理解,但有个老毛病:插值出来的网格会平滑掉原始数据的细节,而且在数据稀疏的地方会出现夸张的“牛眼”或者振荡。三角剖分则完全沿袭原始点的位置,不新造数据,只是把点连成三角形,因此能最大程度保留原始采样的真实性。

Delaunay2D这个滤波器,做的就是把这堆散点连成符合Delaunay准则的三角形网格。它在输入上接受一个vtkPolyData,输出一个覆盖这些点的三角化表面。这个表面虽然只存在于二维投影平面,但完全足够作为切割的对象。

1.2 Delaunay三角剖分好在哪

很多教程会把Delaunay三角剖分降级成“一种三角化方法”,但实际上,它的数学性质决定了它非常适合做切割前的前置处理。

关键性质有两个。第一是空外接圆特性:任意一个三角形的外接圆内,不能包含其他散点。这个特性保证了生成的三角形不会有大角度或者极度狭长的形状,网格整体形态比较均衡。第二是最大化最小角:在所有可能的三角化方案里,Delaunay方案能让所有三角形里最小的那个角尽可能大,换句话说,它能尽量避免“刺状”三角形。

为什么这对篱笆图很关键?因为切割平面和三角形面片求交时,如果三角形特别狭长,交线的位置会非常敏感,很小的数据浮动都会导致切割线抖动明显。均衡的三角形,能够让切割结果稳定、可解释。

提示:第2点虽然表面上看只是网格质量好坏的问题,但放到大量等距切割的场景下,它直接决定最终篱笆图是否有说服力。我在处理一组山区地形点时就遇到过极端瘦长三角形导致某几条切面线“穿模”的情况,后来手动控制了边界点密度才解决。

1.3 篱笆图与等距切割的语义对应

如果已经得到了三角化表面,接下来就很简单:沿着坐标轴的某个方向,每隔固定间距放一个平面,让这个平面与三角网求交,得到的交线就是这一层剖面的轮廓。所有平面切出来的线放到同一个三维场景里,看起来像一排篱笆,错落有致,每一根“栅栏”都是数据在某一个位置上的剖面轮廓。

这里要注意一个概念区分:等距切割的平面不一定是水平面。虽然最常见的是用Z轴作为切割方向,也就是地形上从低到高一层层来切,但只要数据有主轴方向,也可以沿着X轴或者Y轴切,甚至可以旋转任意角度去切。重点是这组平面互相平行,间距相等,这正是“等距”的含义。

2. 管线架构:从散点到篱笆图,每个环节该怎么选

VTK的管线思想在这一整套流程里体现得淋漓尽致:数据流从一个模块流向下一个模块,每个模块只做一件事。要搭出篱笆图,本质上就是串联几条不同的管线,最后把结果合到同一个渲染场景里。

2.1 数据准备:什么样的输入才能直接跑通

Delaunay2D的输入是vtkPolyData,最常见的构造方式是把vtkPoints塞进去。有一点要注意,Delaunay2D是在二维平面上进行三角剖分的,默认工作平面是XY平面。如果原始点的Z值正是你要保留的高程或者场值,那没问题,Z值会作为点属性被保留下来,后续切割时依然能读到。但如果数据本身是三维空间里任意取向的,比如一个倾斜的岩层面,你就得先决定剖分在哪个投影平面做。

处理三维点的通用做法是:先把点变换到一个合适的局部坐标系,让Z值变为0,做完三角剖分,再变换回原坐标系。这个操作可以用vtkTransformFilter实现,也可以用vtkDelaunay2D的SetTransform直接指定一个变换。后者更简单,因为Delaunay2D内部会先应用变换做剖分,再自动逆变换回原始坐标系,输出的点和输入点完全对得上。

代码层面构造点集的时候,还有个容易被忽视的点:vtkPoints没有所谓“执行”的概念,它只是一个数组容器。要把它变成vtkPolyData,得调用SetPoints。如果某些点没有属性值,切出来的交线就只带着空间坐标;如果每个点带了一个标量值(比如温度、硬度),切割结果自动会做插值,交线上每个点的标量值是根据三角形三个顶点的值插出来的,这个特性很有用,后面会讲到怎么利用它做染色。

2.2 切割方向与等距步长怎么定

选切割方向,第一原则是看数据的形态。地形类数据,Z是高度,自然选法向量(0,0,1)的平面去切,得到的是不同高度上的水平剖面;三维模型如果沿某个主轴拉伸,就沿着该主轴方向切,才能看出内部结构变化。

步长的设定,取决于你想在结果里表达多少层信息。我在项目里的一般经验是:先取数据包围盒最长边的1/30到1/50作为初始间距,然后根据切出来的线条数微调。切得太密,场景会变得杂乱,切得太疏,会错过关键变化。比如一个20×20范围的地形,Z跨度大约10,我通常会切12到15层,每层之间的视觉间距刚刚好。

对于斜向切割,vtkPlane有两个参数要设置:法向量和原点。法向量决定了切割方向,原点决定了平面位置。等距切割循环里,每次只改原点位置,法向量保持不变。

2.3 管线串联:Cutter、TubeFilter和AppendPolyData

等距切割过程中,核心模块是vtkCutter。它的工作就是拿一个隐式函数(通常是vtkPlane)去和输入几何体求交。Cutter的输入是Delaunay2D生成的三角网,CutFunction是那个平面,输出是交线(vtkPolyData)。

切出来的线是纯线数据,直接显示的话只有一像素宽,在三维场景里很不显眼,也不方便观察层次。我一般会在Cutter后面接一个vtkTubeFilter,把每一条线“撑”成一根圆管,这样切面轮廓在空间里就有了体积感。

多根管道要放在同一个场景里,可以分别添加到不同的Actor,也可以用一个vtkAppendPolyData把它们合并成一个数据集再显示。合并的好处是后续可以统一用一个Actor渲染,省资源;缺点是如果想给不同高度的线用不同颜色,合并后处理起来就要靠数据里的标量值来区分。

下图是一个简化版的管线示意,说起来就是这个链路:

  • 散点(vtkPolyData) → vtkDelaunay2D → 三角网格
  • 三角网格 + vtkPlane → vtkCutter → 等距交线
  • 等距交线 → vtkTubeFilter → 带体积的管道
  • 多组管道 → vtkAppendPolyData → 单个Actor显示

这套管线看起来平淡无奇,但顺着走下来,你就能很清晰地把“原始点”和“结果线”对应起来,排查问题也方便:哪一步结果不对,就在哪一步后设置一个临时Actor输出看状态。

3. 直接可跑的示例:C++实现等距切割篱笆图

下面给一段我自己常用的示例,基于C++与VTK原生接口。用Python的朋友也别觉得陌生,VTK的Python绑定只是把类名保持原样,流程和参数完全一样,甚至更简洁。

3.1 一个模拟的双峰曲面作为测试数据

先用代码构造一个模拟的“双峰”表面:在40×40的点阵上,生成两个高斯型的隆起。这样做的好处是曲面变化丰富,切出来的篱笆图层次分明,观察效果直观。

#include <vtkSmartPointer.h> #include <vtkPoints.h> #include <vtkPolyData.h> #include <vtkDelaunay2D.h> #include <vtkPlane.h> #include <vtkCutter.h> #include <vtkTubeFilter.h> #include <vtkAppendPolyData.h> #include <vtkPolyDataMapper.h> #include <vtkActor.h> #include <vtkRenderer.h> #include <vtkRenderWindow.h> #include <vtkRenderWindowInteractor.h> #include <cmath> int main(int, char*[]) { // 1. 构造散点数据 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New(); const int n = 40; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { double x = i * 0.5; double y = j * 0.5; double r1 = (x - 10) * (x - 10) + (y - 10) * (y - 10); double r2 = (x - 5) * (x - 5) + (y - 5) * (y - 5); double z = 6.0 * std::exp(-r1 / 20.0) + 4.0 * std::exp(-r2 / 10.0); points->InsertNextPoint(x, y, z); } } vtkSmartPointer<vtkPolyData> cloud = vtkSmartPointer<vtkPolyData>::New(); cloud->SetPoints(points); // 2. 执行Delaunay三角剖分 vtkSmartPointer<vtkDelaunay2D> del = vtkSmartPointer<vtkDelaunay2D>::New(); del->SetInputData(cloud); del->SetTolerance(0.001); del->Update(); // 3. 计算Z的范围,确定等距切面的位置 double bounds[6]; del->GetOutput()->GetBounds(bounds); double zMin = bounds[4]; double zMax = bounds[5]; const int slices = 15; double dz = (zMax - zMin) / slices; // 4. 循环执行等距切割 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New(); for (int s = 0; s <= slices; ++s) { double z = zMin + s * dz; vtkSmartPointer<vtkPlane> plane = vtkSmartPointer<vtkPlane>::New(); plane->SetNormal(0, 0, 1); plane->SetOrigin(0, 0, z); vtkSmartPointer<vtkCutter> cutter = vtkSmartPointer<vtkCutter>::New(); cutter->SetInputConnection(del->GetOutputPort()); cutter->SetCutFunction(plane); cutter->SetGenerateCutScalars(1); cutter->Update(); append->AddInputData(cutter->GetOutput()); } // 5. 把切出来的线变成管道,便于三维观察 vtkSmartPointer<vtkTubeFilter> tube = vtkSmartPointer<vtkTubeFilter>::New(); tube->SetInputConnection(append->GetOutputPort()); tube->SetRadius(0.08); tube->SetNumberOfSides(8); // 6. 渲染 vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New(); mapper->SetInputConnection(tube->GetOutputPort()); vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New(); actor->SetMapper(mapper); vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New(); renderer->AddActor(actor); renderer->SetBackground(0.1, 0.1, 0.2); vtkSmartPointer<vtkRenderWindow> window = vtkSmartPointer<vtkRenderWindow>::New(); window->AddRenderer(renderer); window->SetSize(800, 600); vtkSmartPointer<vtkRenderWindowInteractor> interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New(); interactor->SetRenderWindow(window); interactor->Start(); return 0; }

3.2 这段代码里的关键细节

这个示例跑起来后,你会看到一串环形的管道横亘在三维空间里,每个高度一层,形状从低处的完整双环逐渐变化到高处只剩两个尖顶,信息非常直观。整个过程里,有几个地方容易理解错,我单独拿出来讲。

第一,为什么切割后要用Tick(管道化)而不是直接显示线。线在VTK渲染里默认宽度是1个像素,在密密麻麻的多层切割结果中,视觉上会糊成一团。管道化之后,你可以通过调节半径来控制视觉权重,让不同层面的切面轮廓之间有明确的物理距离感。

第二个关键是SetGenerateCutScalars。这一项打开后,Cutter会在输出的每个点上生成一个标量值,这个值代表“该点到原始表面某条边的垂直距离”,简单说就是切割距离。它的实际用途是给交线染色。比如,你可以让越靠近上层的线颜色越偏暖,越靠近下层的线颜色越偏冷。默认情况下这个选项是关闭的,我曾经因为没打开它,导致所有线条都是同一种颜色,结果看不出层次,排查了半天才发现是这个开关的问题。

第三,关于循环里反复创建vtkCutter的做法。每创建一个新cutter,其实是在管线里添加一个新的分支,它们都共享同一个上游del的三角网格。VTK内部有缓存机制,第一次计算Delaunay剖分之后,后续的切割分支不会重复执行三角化,所以性能上完全不用担心。如果你想把代码写得更紧凑,也可以只用一个cutter对象,循环里改掉cutFunction的Origin之后重新Update,效果等价。

第四,处理Z值方向的注意事项。我这里的模拟数据,Z值天然就是高程,所以Delaunay2D默认在XY平面剖分毫无问题。但如果你拿到的是三维点云,比如一个倾斜的矿体表面,直接按上面的代码跑,三角化会忽略Z方向的拓扑,产出错误结果。这时你需要在使用Delaunay2D之前,把点云做一个坐标变换,让它们落在某个投影平面上,等切完再变换回去。这也是前面提到的SetTransform方法的典型应用场景。

3.3 如果你想用Python,代码怎么改

Python里这套流程看起来更精炼,核心逻辑完全照搬:

import vtk import numpy as np # 生成散点 n = 40 xs, ys = np.meshgrid(np.linspace(0, 20, n), np.linspace(0, 20, n)) r1 = (xs - 10)**2 + (ys - 10)**2 r2 = (xs - 5)**2 + (ys - 5)**2 zs = 6 * np.exp(-r1 / 20) + 4 * np.exp(-r2 / 10) points = vtk.vtkPoints() for i in range(n): for j in range(n): points.InsertNextPoint(xs[i, j], ys[i, j], zs[i, j]) cloud = vtk.vtkPolyData() cloud.SetPoints(points) # 三角化 del_ = vtk.vtkDelaunay2D() del_.SetInputData(cloud) del_.SetTolerance(0.001) del_.Update() # 等距切割 z_min, z_max = del_.GetOutput().GetBounds()[4], del_.GetOutput().GetBounds()[5] slices = 15 dz = (z_max - z_min) / slices append = vtk.vtkAppendPolyData() for s in range(slices + 1): plane = vtk.vtkPlane() plane.SetNormal(0, 0, 1) plane.SetOrigin(0, 0, z_min + s * dz) cutter = vtk.vtkCutter() cutter.SetInputData(del_.GetOutput()) cutter.SetCutFunction(plane) cutter.SetGenerateCutScalars(1) cutter.Update() append.AddInputData(cutter.GetOutput()) # 管道化与渲染 tube = vtk.vtkTubeFilter() tube.SetInputConnection(append.GetOutputPort()) tube.SetRadius(0.08) tube.SetNumberOfSides(8)

Python的接口和C++几乎一一对应,一个模块对应一个类,没有魔法。唯一要留意的是,设置为InputData时,如果上游数据还没Update,切出来的可能是空数据,所以上面的代码里在Delaunay2D后显式调用了Update()。

4. 等距切割的隐藏坑和排错方法

这个流程看着短,但实际跑起来,我见过的坑比想象中的多。有些问题第一次遇到时很让人抓狂,但搞清楚原因之后会发现其实都很简单。我把常见的几种现象和排查思路整理成一张表,顺便挑几个好好说说。

4.1 常见异常现象对照表

现象可能原因解决思路
切出来的线全部在Z=0附近坐标点带有Z,但Delaunay2D把Z当成了纯属性,三角网格的几何Z值被丢失用Transform让点先落到XY平面,或检查SetTransform用法
等距线在边缘处飞出预期范围三角剖分的边界出现极长三角形,边界梯度大控制边界点,或用裁剪滤波器预先剔除边缘三角
切割结果为空Cutter输入未连接或未Update检查过滤器管线的执行时机,确保上游先Update
线条颜色没有层次忘记打开SetGenerateCutScalars打开CutScalars选项
多条切割线在空间里粗细不一致单条线管道半径设置过小,某些区域视觉上明显统一半径,或使用LineWidth渲染模式
切割线在接缝处断开Denaunay网格上有非流形边缘,Cutter求交出现多个线段检查原始点是否有重复点,先通过CleanPolyData合并

此外还有一个很普遍的问题:数据里如果包含大量重复或者几乎重叠的点,Delaunay2D的三角剖分会变得不稳定。重复点会导致三角形退化,切割线也会出现毛刺。遇到这种情况,先给点云过一遍vtkCleanPolyData,把重复点合并掉,基本能消除大部分莫名其妙的问题。

4.2 排查方法:从哪一步下手

VTK的好处是每步都可以可视化。排查篱笆图异常,我个人的习惯是分三段检查。

第一段看三角网格是否正常。直接把Delaunay2D的输出去到一个Mapper渲染,确认三角网格是一个覆盖范围合理、没有大空洞的表面。如果这一步就看到明显的杂乱三角形,说明原始点数据可能有问题,比如点排序、重复点、点坐标范围异常。

第二段看单张切割线是否正常。把循环先改成只切割一个固定平面,比如Z等于中间值的那一层,然后单独查看cutter输出。这时候就能验证切割方向、法向量设置是否正确。我遇到过一种情况:法向量设成了(0,1,0),但数据其实沿着Y方向也有高度起伏,导致切割出来的线与印象中的地形轮廓完全对不上。把这一步单独跑一遍,问题马上会暴露。

第三段检查多个切割线合并后的空间关系。如果单张线都正常,合到一起却显得杂乱,就要检查管道的半径和透明度设置,很可能只是视觉参数不合适,而不是数据出错。

4.3 等距切割时Z向与X/Y向的直觉误区

做地形时,大家都习惯用Z方向切,切出来的是“等高线”的伴侣——等高水平切面。但篱笆图的概念其实比等高线更通用。比如你要研究一条断层在空间上的延展规律,可以沿着断层走向方向做等距切割,这时候平面法向量不再垂直于Z轴,而是垂直于切割方向。

我自己的一个项目里,处理的是河谷沉积层数据。数据的空间主轴大致沿着河流方向弯曲。如果只按Z轴切,得到的剖面和河谷走向斜交,看不出沉积层在横截面上的真实形态。后来改成沿河谷方向等距切,绕数据主轴放了一组法线固定的平面,结果出来的篱笆图清晰反映了沉积层厚度的空间变化。这说明,切割方向的选定本质上是一个“围绕数据形态做投影”的问题。

4.4 另一个容易翻车的细节:等距步长的枚举例

循环写s从0到slices时,很多人会纠结到底切slices层还是slices+1层。这取决于你的语义。如果数据Z范围是[zMin, zMax],等距切N层,dz = (zMax-zMin)/N,那你会得到N+1条线,因为把两端都包含进去了。如果你希望刚好N条线,应该让起点从zMin+dz/2开始,这样两端各留半个步长。

要不要包含两端,看你的应用场景。做地质剖面展示时,我通常包含两端,从最顶层一直切到最底层,视觉上有“盖帽”的感觉,边界清晰。做统计特征分析时,用中间位置错开端点更合理,避免端点处的异常值干扰结果。

5. 进一步提高篱笆图的表现力

完成等距切割只是第一步,篱笆图真正能发挥生产力,往往在于最后的呈现。从纯技术角度讲,输出一组线已经达到了“能用”的标准;但从工程交付的角度看,一份好的篱笆图应当让人一眼看懂变化趋势,而不是靠在屏幕上旋转半天才能看出门道。

5.1 让每层切面拥有自己的色彩逻辑

最常见的增强方式是利用vtkCutter输出点上的CutScalars给线条做伪彩色映射。你可以接一个vtkLookupTable,自定义颜色范围:低层的线用冷色表示,高层的线用暖色表示,这样不同高度切面的分布一目了然。

如果有更复杂的业务语义,比如每条切割线代表一种岩性层,那就最好给每一层切面单独生成一个vtkActor,并手动设置颜色。注意,多条切割线合到一个AppendPolyData里之后,再用vtkCellData去区分层次比较麻烦;更稳妥的方法是循环里就为每个切面生成独立的Actor,分别上色,最后全部AddActor进渲染器,这样每条线之间可以设置透明度,叠加起来还能看到前后遮挡关系。

5.2 结合原始表面的半透明渲染

篱笆图本身是“骨架”性质的结果。让它更丰满的一个经典玩法,是将Delaunay2D得到的原始三角面片以半透明的形式一起渲染。三角网格用半透明材质,等距切割线用实心管道,两者叠加,既能看整体曲面形态,又能看到每个剖面的局部轮廓。

代码上,只需要增加一个actor给原始的del输出,设置透明度在0.2到0.4之间,然后和切出来的管道一起放进同一个Renderer。这个做法对调试和汇报都特别好用,听众能直观地理解“这些线是从哪来的”。

5.3 从篱笆图继续拓展的方向

如果这套流程已经稳定跑通了,接下来有几个很自然的扩展方向。第一,把等距线进一步填充成真正的纵剖面网格面,也就是把两条相邻的切面轮廓用四边形拼接起来,形成完整的剖面图。这在很多专业软件里叫纵向剖面模型。第二,将点云上的标量属性(比如硬度、密度)传递到切面线上,再以切面线作为边界做栅格填充,就能生成连续剖面云图。第三,对多组切面做动画生成,按顺序显示切面位置,便于演示数据在空间中的渐进变化。

这些扩展都不会偏离等距切割这条主线。这也说明篱笆图并不是一个孤立的工具,它更像是把复杂三维数据“降维”成若干可读剖面的一套通用思路。

在实际操作中,我最后常做的一个小优化是把切割线的管径调整到和数据尺度匹配。在我常用的那套地形数据里,包围盒对角线大约30个单位,半径设成0.08刚好合适;但如果你换成毫米级的零件扫描数据,0.08的管道就会把整个画面糊死。一个稳妥经验是先跑一次看整体,再调整TubeFilter的半径,直到线条既不突兀又能一眼看清为止。这个微调虽然不起眼,却是整个流程里最体现手感的一步。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询