Delaunator 库指南：快速平面 Delaunay 三角剖分

作者：Aud42nce
2023-07-05
11 6 0

编者按

本文来自游戏设计博客 Red Blob Games，由 indienova 取得授权并译制发表，原文链接见文末。

  • 原文作者：Amit Patel
  • 译者：Aud42nce

正文

JavaScript 库 Delaunator 能够快速进行平面上的 Delaunay 三角剖分。它的输入是平面内的点集：

Delaunator 的输入

输出是平面内的一种三角剖分：

Delaunator 的输出

输出结果会存储为紧凑数组（compact array），数组类型为整型。相较其他表示方法， 这样似乎不太方便，但这正是 Delaunator 库能够高效运行的原因。

Delaunay 三角形

我们通过语句 delaunay = Delaunator.from(points) 可以构建出一个 delaunay 对象，它包含 triangles （三角形）和 halfedges （半边）两个数组，两个数组都以半边的 id 为索引。什么是半边呢？

半边

三角形的一条边有可能同时属于另一个三角形。我们将同一条边分别表示成 A→B 和 B→A 两个半边，而不考虑 A↔︎B。使用两个半边是实现 Delaunator 库各种功能的关键所在。

delaunator 输出中的两个数组都以半边 e 为索引：

  • delaunay.triangles[e] 返回半边 e 起点的 id
  • delaunay.halfedges[e] 返回相邻三角形中对应的半边，如没有相邻三角形则返回 -1

同时，三角形的 id 和构成它的半边的 id 之间存在如下关系：

  • 对于三角形 t，构成它的半边分别为 3t, 3t +1 以及 3*t+2。
  • 对于半边 e，它所属的三角形为 floor(e/3)。

为此我们可以构建如下辅助函数：

function edgesOfTriangle(t) { return [3 * t, 3 * t + 1, 3 * t + 2]; }
function triangleOfEdge(e)  { return Math.floor(e / 3); }

此外，我们也可以构建辅助函数，帮助我们通过一条半边找到同一三角形内的前一条及后一条半边：

function nextHalfedge(e) { return (e % 3 === 2) ? e - 2 : e + 1; }
function prevHalfedge(e) { return (e % 3 === 0) ? e + 2 : e - 1; }

注：页面提供的例码仅作为阅读参考，并不用于实际运行。


Delaunay 边

我们无需构建三角形对象，就能绘制出三角形的所有边。每条边都可以被表示为两个半边。半边 e 的起点为 points[delaunay.triangles[e]]，与它反向的另一个半边 delaunay.halfedges[e] 则以另外一端为起点，我们以此得知该边的两个端点。然而，沿凸包外侧的半边不会有相反的半边，此时 delaunay.halfedges[e] 为 -1，因此无法找到 points[delaunay.halfedges[e]] 对应的点。为了确保能顺利找到边的另一个端点，我们需要使用 points[nextHalfedge(e)] 而非 points[delaunay.halfedges[e]] 。我们可以遍历所有半边，并选择其中一半来绘制出图像：

function forEachTriangleEdge(points, delaunay, callback) {
    for (let e = 0; e < delaunay.triangles.length; e++) {
        if (e > delaunay.halfedges[e]) {
            const p = points[delaunay.triangles[e]];
            const q = points[delaunay.triangles[nextHalfedge(e)]];
            callback(e, p, q);
        }
    }
}
绘制三角形的边


构建三角形

一个三角形由三个连续半边 3*t, 3*t + 1, 3*t + 2 组成。每条半边 e 的起点为 points[e]，因此，将这三个点相连就能得到一个三角形。

function edgesOfTriangle(t) { return [3 * t, 3 * t + 1, 3 * t + 2]; }

function pointsOfTriangle(delaunay, t) {
    return edgesOfTriangle(t)
        .map(e => delaunay.triangles[e]);
}

function forEachTriangle(points, delaunay, callback) {
    for (let t = 0; t < delaunay.triangles.length / 3; t++) {
        callback(t, pointsOfTriangle(delaunay, t).map(p => points[p]));
    }
}
绘制三角形


邻接三角形

我们还可以使用三角形的半边来找到相邻的三角形。由于每个半边对应的另一半边（若存在）处于相邻的三角形中，因此我们可以借助辅助函数 edgeIdToTriangleId 找到另一半边所处的三角形：

function triangleOfEdge(e)  { return Math.floor(e / 3); }

function trianglesAdjacentToTriangle(delaunay, t) {
    const adjacentTriangles = [];
    for (const e of edgesOfTriangle(t)) {
        const opposite = delaunay.halfedges[e];
        if (opposite >= 0) {
            adjacentTriangles.push(triangleOfEdge(opposite));
        }
    }
    return adjacentTriangles;
}

Voronoi 胞

将 Delaunay 三角形的外心相连，就能得到对应的 Voronoi 图。我们可以利用 Delaunay 图的对偶来完成这一工作。步骤如下：

  1. 计算出每一个三角形的外心
  2. 根据两个外心构建 Voronoi 边
  3. 将上述边连接成 Voronoi 胞


三角形的外心

关于三角形外心的计算方法见维基百科。外心通常在三角形内部，但并不一定。

function circumcenter(a, b, c) {
    const ad = a[0] * a[0] + a[1] * a[1];
    const bd = b[0] * b[0] + b[1] * b[1];
    const cd = c[0] * c[0] + c[1] * c[1];
    const D = 2 * (a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]));
    return [
        1 / D * (ad * (b[1] - c[1]) + bd * (c[1] - a[1]) + cd * (a[1] - b[1])),
        1 / D * (ad * (c[0] - b[0]) + bd * (a[0] - c[0]) + cd * (b[0] - a[0])),
    ];
}
Delaunay 三角形的外心

利用下列函数，能够很方便通过三角形的 id 找到它的外心。

function triangleCenter(points, delaunay, t) {
    const vertices = pointsOfTriangle(delaunay, t).map(p => points[p]);
    return circumcenter(vertices[0], vertices[1], vertices[2]);
}


Voronoi 边

有了外心，我们可以在不构建多边形的情况下绘制 Voronoi 边。每个 Delaunay 三角形的半边对应一个 Voronoi 多边形的半边。Delaunay 半边连接 delaunay.triangles[e] 和 delaunay.triangles[nextHalfedge(e)] 这两个点。 Voronoi 半边连接两个三角形的外心 triangleOfEdge(e) 和 triangleOfEdge(delaunay.halfedges[e])。我们可以遍历半边并构建这些线段：

function forEachVoronoiEdge(points, delaunay, callback) {
    for (let e = 0; e < delaunay.triangles.length; e++) {
        if (e < delaunay.halfedges[e]) {
            const p = triangleCenter(points, delaunay, triangleOfEdge(e));
            const q = triangleCenter(points, delaunay, triangleOfEdge(delaunay.halfedges[e]));
            callback(e, p, q);
        }
    }
}
绘制 Voronoi 边


构建 Voronoi 胞

想要构建多边形，我们需要找出共用同一顶点的所有三角形。利用半边这一数据结构能够完成我们的目标。假设我们已经有了一条指向共用顶点的半边。我们可以轮流执行下列两步来找出所需的其它部分：

  1. 调用 nextHalfedge(e) 找到当前三角形中下一个从该顶点出发的半边
  2. 调用 halfedges[e] 找到相邻三角形中指向该顶点的半边
围绕顶点执行算法的过程
function edgesAroundPoint(delaunay, start) {
    const result = [];
    let incoming = start;
    do {
        result.push(incoming);
        const outgoing = nextHalfedge(incoming);
        incoming = delaunay.halfedges[outgoing];
    } while (incoming !== -1 && incoming !== start);
    return result;
}

注意，上面的算法至少需要我们先得到任意一条指向该顶点的半边。为了能够快速找到指向特定节点的半边，不妨为所有半边建立索引。比如本文末尾的改版 forEachVoronoiCell 就是一个简单的例子。


绘制 Voronoi 胞

要想绘制出 Voronoi 胞，我们可以通过指向顶点的半边来找到对应三角形，再找到这些三角形的外心。我们可以遍历所有半边，不过，由于多条半边可能会指向同一顶点，我们需要记录哪些顶点已经被访问过。

function forEachVoronoiCell(points, delaunay, callback) {
    const seen = new Set();  // of point ids
    for (let e = 0; e < delaunay.triangles.length; e++) {
        const p = delaunay.triangles[nextHalfedge(e)];
        if (!seen.has(p)) {
            seen.add(p);
            const edges = edgesAroundPoint(delaunay, e);
            const triangles = edges.map(triangleOfEdge);
            const vertices = triangles.map(t => triangleCenter(points, delaunay, t));
            callback(p, vertices);
        }
    }
}
绘制 Voronoi 胞


凸包

上面的 edgesAroundPoint 循环存在一处问题。凸包上的顶点并未完全被三角形包围，所以在循环途中有可能返回-1，导致循环中途停止。以下是三种处理方法：

  1. 忽略这个问题。确保不会围绕凸包上的顶点来执行该算法。
  2. 调整代码。
    • 在所有和 halfedges 有关的代码中找到和 -1 相关的部分
    • 让上述 edgesAroundPoint 循环从整个三角剖分结果的”最左侧“开始运行。此时当算法返回 -1 时，所有的三角形都已经被访问过了。
  3. 调整原始数据。构造网格的"背面"来消除凸包的情形。这时就不会再有返回 -1 的 halfedges 出现了。
    • 在每个返回-1 的半边处增加一个“幽灵（ghost）”半边。
    • 在无限远处增加一个“幽灵”顶点，用以代表整个三角剖分的“背面”。
    • 通过“幽灵“顶点和这些“幽灵”半边来构建“幽灵”三角形。

下面是寻找“最左侧”半边的算法的一则示例：

function forEachVoronoiCell(points, delaunay, callback) {
    const index = new Map(); // point id to half-edge id
    for (let e = 0; e < delaunay.triangles.length; e++) {
        const endpoint = delaunay.triangles[nextHalfedge(e)];
        if (!index.has(endpoint) || delaunay.halfedges[e] === -1) {
            index.set(endpoint, e);
        }
    }
    for (let p = 0; p < points.length; p++) {
        const incoming = index.get(p);
        const edges = edgesAroundPoint(delaunay, incoming);
        const triangles = edges.map(triangleOfEdge);
        const vertices = triangles.map(t => triangleCenter(points, delaunay, t));
        callback(p, vertices);
    }
}

然而，即便有了这些改动，沿凸包构建 Voronoi 胞仍然需要先将边映射到外部再进行剪裁。Delaunator 库并没有提供这些功能。有需要的话可以选择使用 d3-delaunay 库。

小结

Delaunator 库使用半边表示点和三角形之间的关系。本页的示例函数显示了如何在对象类型之间移动：

原文链接：https://mapbox.github.io/delaunator/
*本文内容系作者独立观点，不代表 indienova 立场。未经授权允许，请勿转载。