65.9K
CodeProject 正在变化。 阅读更多。
Home

Delaunay 三角剖分实现快速网格生成

HoshiKata
starIconstarIconstarIconstarIconstarIcon

5.00/5 (23投票s)

2012 年 11 月 12 日

CPOL

8分钟阅读

viewsIcon

115247

downloadIcon

6003

实际即时从任意点快速生成网格。

下载 Delaunay.zip

介绍  

十年前，实时计算曲面网格尚不现实，并且也无法获得可定制的源代码模块。有更快的版本，但它们是大型实现，难以阅读和修改。  借助现代计算机和现代语言，不仅可以实时生成简单的网格，而且代码的编写方式也易于理解。作者认为，计算几何和计算机视觉正进入一个新时代，实时处理对于越来越多的问题来说是可行的。  该代码设计用于可重用性，以便一个顶点可以同时用于生成网格，还可以作为其他数据结构的一部分并在给定点跟踪其他方面。 

为了形成一个可以被显卡着色的曲面，目标是创建一个顶点列表，以及一个构成逆时针三角形的顶点三元索引列表。无论是可视化有限元建模 (FEM) 求解的微分方程，还是查看 3D 物体的形状，第一步都是形成三角形网格。  通常选择一组点很容易，对于圆柱体或球体等简单物体，点生成直接导致三角形生成。从任意点集生成网格是 Delaunay 三角剖分发挥作用的地方。  代码的目标不是计算理想的三角剖分，而是为大多数实际的 3D 或 FEM 问题提供一个“足够好”的解决方案，并快速完成。因此，代码在添加每个点后递归地改进三角剖分，直到达到递归计数或它是一个理想的 Delaunay 三角剖分。这限制了网格的计算时间，并避免了在出现软件错误时无限递归。

背景

Delaunay 三角剖分以 Boris Delaunay 命名。其基本思想是形成一个网格，其中每个三角形的三个点位于一个不包含任何其他点的圆的边缘上。这意味着给定任意两个相邻的三角形（四边形），位于分隔线相对角的度数之和小于 180 度。


三角形网格 a,b,d 和 b, c, d 不满足 Delaunay 条件。但 b,c,a 和 c,d,a 满足。这是因为角 bcd + dab > 180，而 abc+cda < 180。这迫使网格中的三角形尽可能均匀分布。

Bowyer Watson 算法 通过一次添加一个点来迭代地将网格从一个 Delaunay 网格转换为另一个。通用算法通过删除所有其外接圆包含新点的三角形并重新形成网格来工作。提供的代码采用一种更简单的方法。 

从上到下，三角形被翻转以满足 3 个点的 Delaunay 条件

这做了同样的事情，但有七个点。请注意，最初三角形大多是细长条，而其他三角形非常大。

 在每次插入后进行一次翻转后，网格更加规则（如下）。

经过第二次翻转后，第一次翻转的效果将被翻转（如下）。

设计 

网格生成

Bower-Watson 算法有点复杂，因为： 

  • 添加一个点可能会导致一组三角形被删除并重新形成。
  • 一次重新形成一组会尝试形成许多可能的三角剖分。
  • 必须搜索受影响的三角形。 
  • 调整网格没有一个有序的模式。  

常用的“凸包”技术也从一个点螺旋向外，从内到外重新形成三角剖分。这种方法很好，因为它对计算和重新计算网格采取了有序的方法。缺点是软件正在实时重新计算每个三角形的圆并对三角形进行排序。  

提供的代码采用一种混合方法，它形成一个三角形的链表。每个三角形由三个顶点组成，并且每条边都有一个指向相邻三角形（或 null）的指针。代码具有以下优点：  

  • 每个三角形都有一个指向其邻居的链接。   
  • 如果一个三角形被修改，递归会向外探索，它的邻居可能会受到影响（没有圆）。
  • 三角形通常不会被删除，而是被添加或修改。 
  1. 计算边界矩形： 

    2. 请注意，a->b->c 是逆时针的，因此如果它有一个面由 direct3D 等渲染，它将从页面朝外。三角形 c->b->a 是顺时针的，因此它的面朝向页面。只要方向保持不变，起始点就无关紧要，因此 bdc 与 dcb 或 cbd 相同。抽象代数中的这个概念通常称为交换律。

    此时内存中： 

    三角形 1 = a,b,c，边 b,c = 三角形 2。

    三角形 2 = bdc，边 c,b = 三角形 1。

  2. 找到包含新点的三角形（例如，abc），然后用三个新三角形替换旧的三角形 1。 

  3. 修复指针： 

    4. abe 边 be 应指向 ebc。

    Abe 边 ea 应指向 aec。

    等等。 

    然后修复那些以前指向 abc 的旧相邻三角形：bdc 边 bc 以前指向 abc，使其指向 ebc。

  4. 检查每个三角形及其外部邻居：例如 bec 与 bdc 共边。

  5. 检查相对相邻边的角度。如果大于 180 度，则翻转相邻线。

  6. 如果三角形未翻转，请跳至步骤 10。  

  7. 将旧三角形重新形成新三角形，修复三角形边。 

  8. 对于每个三角形及其相邻三角形，递归重复步骤 5-8。 

三角形包含一个点

算法的关键是找到要修改的正确三角形。

经典上，一个三角形包含一个点如果它在内部侧。如果一条线通过两个点，该线通常表示为 y = mx + b，其中 m 是斜率，b 是 x=0 截距。线上方区域可以表示为 y > mx+b。因此，如果一个三角形有边 s, r，而另一个点是 t。那么如果测试点 p 位于不等式的同一侧，则测试点 p 位于与另一个顶点同一侧的线上。

这意味着要计算一个点 p 是否在内部，代码会将其与边对应顶点进行比较，看它是否在同一侧。这需要大量计算。缓存信息可以加快比较速度，但不足以完全解决问题。  大多数三角形不包含该点，并且该点也不在三角形附近。为了加速计算，只有靠近三角形的点才需要使用经典方法进行测试。软件计算每个顶点与三角形中心之间的距离。如果点离三角形中心不够近，则无需测试。

加速

三角形对象会缓存按需计算的信息以形成网格。当一个三角形与其邻居“翻转”时，它的两条边会发生变化，但两个点保持不变。当一个顶点移动时，它会使这些边的现有信息失效。

例如，如果任何一个顶点发生变化，三角形的中心就会发生变化，设置点 A、B、C 会将计算值设置为 false，请求中心时会计算它并设置标志。

Vertex m_Center;

public Vertex Center
{
    get
    {
        if (m_centerComputed) return m_center;
        m_center = new Vertex(
            (A.X + B.X + C.X) / 3f,
            (A.Y + B.Y + C.Y) / 3f,
            (A.Z + B.Z + C.Z) / 3f);

        float delta = m_center.DeltaSquared(A);
        float tmp = m_center.DeltaSquared(B);
        delta = delta > tmp ? delta : tmp;
        tmp = m_center.DeltaSquared(C);
        delta = delta > tmp ? delta : tmp;
        FarthestFromCenter = delta;
        m_centerComputed = true;

        return m_center;
    }
}

一旦计算出中心，软件就可以重用这些信息，并且仅在需要时重新计算。

使用代码

该表单提供了一个生成半聚集种子点的演示。Mesh 类具有以下 API：

  • Points – 网格中使用的 Vertex 对象列表。
  • Facets – 网格的三角形面的列表。
  • Bounds – 用于形成初始方形区域的矩形边界。
  • Recursion – 从初始更改的三角形开始迭代的次数。对于大多数小于 100 个点的网格，需要递归重新三角剖分超过 5 或 6 次的几率非常低。这允许用户平衡计算时间和网格质量。
  • Compute() - 从顶点列表中计算网格。
  • Append() - 向三角剖分添加一个新顶点。 
  • Setup() - 设置边界。 
  • Draw() - 将网格绘制到 System.Drawing.Graphics 对象以进行调试。 
  • GetVertexIndicies() - 获取构成逆时针面的索引集（每个三角形 3 个）。 

示例： 

List<Vertex>
set = new List<Vertex>();
Random.Random r = new Random.Random(seed);
for(int i = 0; i < count; i++)
{
   set.Add(new Vertex(r.Float(0,width), r.Float(0, height), 0));
} 
Mesh m = new Mesh();
m.Recursion = 6;
m.Compute(set, new RectangleF(0, 0, 640, 480));            
 
m.Draw(g, 0, 0, 640, 680);
List<Vertex> p = m.Points;               // Points for directX vertex buffer.
int
[] indicies = m.GetVertexIndicies(); // Indexes for directX surface.

性能  

在我的笔记本电脑（一台不错的 I7，带 Windows 7）上，该算法在小型网格上的运行速度足够快，可以达到视频速度。每个点都需要搜索需要拆分成三个的点，因此时间随搜索而增长。我的目标是将其用于另一个项目，在该项目中我需要能够重用代码、维护它并使其运行得相当快。虽然这个速度相对较快，但也有更快的，如果您需要超快的速度，OpenCV 提供了该算法的一个版本。

50 Vertexes 5 msec
100 Vertexes 11 msec
500 Vertexes 129 msec
1k Vertexes 380 msec
5k Vertexes 7.8 seconds
10k Vertexes 31 seconds
60 Hz 16 msec 150 vertexes
30 Hz 33 msec 230 vertexes 

20 Hz 50 msec 280 vertexes
© . All rights reserved.