[3D基础]Delaunay与Triangulate构网

Delaunay三角剖分（Delaunay Triangulation）和Triangulate构网是两种常用于生成三角网格的方法，它们都有其独特的特点和应用场景。

Delaunay三角剖分：

Delaunay三角剖分是一种经典的三角剖分方法，其特点是任意三角形的外接圆不包含任何其他顶点。换句话说，对于给定的点集，Delaunay三角剖分是使得生成的三角形尽可能接近等边三角形，并且最大化了三角形的最小角。这种特性使得Delaunay三角剖分在许多应用中都具有优越性，例如地理信息系统、有限元分析、计算机图形学等领域。

Delaunay三角剖分的算法有多种实现方式，包括增量构建法、分治法、凸壳法等。其中增量构建法是最常见的一种，其基本思想是逐步将点集中的点加入到一个初始空的三角形网格中，使得生成的三角形满足Delaunay条件。

Triangulate构网：

Triangulate构网是一种简单而直观的三角剖分方法，其基本思想是将给定的区域分割成许多小三角形。Triangulate方法通常用于离散化几何图形或表面，生成用于渲染或分析的三角网格。

Triangulate构网的实现方式通常比Delaunay三角剖分简单，因为它不需要满足Delaunay条件。常见的Triangulate构网方法包括：最近邻三角化、四边形划分、有限元分析中的特定方法等。这些方法可以根据应用的需求和性能要求选择合适的算法。

总的来说，Delaunay三角剖分适用于需要满足特定条件（如Delaunay条件）的应用，而Triangulate构网适用于一般的离散化和三角化需求。选择合适的方法取决于具体的应用场景、性能要求和计算资源等因素。

Delaunay.js 是一个 JavaScript 库，用于在浏览器中执行 Delaunay 三角剖分。它提供了一种方便的方法来将点集转换为三角形网格，这在许多领域中都很有用，例如计算机图形学、地理信息系统、物理模拟等。

以下是对 Delaunay.js 库的详细解释：

功能特点：
- Delaunay.js 库提供了一个方便的接口，允许你将给定的点集转换为 Delaunay 三角剖分。
- 它可以处理任意数量的点，即使是大型点集也可以快速处理。
- 该库还提供了用于检查给定点是否在三角形内部的方法，以及查找包含给定点的三角形的方法。
- Delaunay.js 还支持 Voronoi 图的计算，允许你根据 Delaunay 三角形网格构建 Voronoi 图。
使用方法：
- 首先，在你的 HTML 页面中引入 Delaunay.js 库的 JavaScript 文件。
- 使用 Delaunay.from() 方法将给定的点集转换为 Delaunay 三角剖分。
- 一旦得到了三角剖分对象，你可以使用它的方法执行各种操作，比如查找包含给定点的三角形、获取三角形的顶点等。
- 如果需要构建 Voronoi 图，可以使用 Delaunay.voronoi() 方法来获取 Voronoi 图对象。
示例用途：
- Delaunay.js 可以用于创建网格化地理数据的可视化效果，例如地图和地形图。
- 它可以用于模拟物理系统中的三角形网格，例如流体模拟、粒子系统等。
- 该库还可以用于计算三角形网格上的几何性质，例如计算三角形的面积、周长等。

总的来说，Delaunay.js 是一个功能强大且易于使用的 JavaScript 库，适用于执行 Delaunay 三角剖分和 Voronoi 图计算。它在许多领域中都有广泛的应用，并且提供了丰富的功能和灵活的接口，使其成为处理三角网格数据的理想选择。


var Delaunay;

(function() {
    "use strict";

    var EPSILON = 1.0 / 1048576.0;

    function supertriangle(vertices) {
        var xmin = Number.POSITIVE_INFINITY,
            ymin = Number.POSITIVE_INFINITY,
            xmax = Number.NEGATIVE_INFINITY,
            ymax = Number.NEGATIVE_INFINITY,
            i, dx, dy, dmax, xmid, ymid;

        for(i = vertices.length; i--; ) {
            if(vertices[i][0] < xmin) xmin = vertices[i][0];
            if(vertices[i][0] > xmax) xmax = vertices[i][0];
            if(vertices[i][1] < ymin) ymin = vertices[i][1];
            if(vertices[i][1] > ymax) ymax = vertices[i][1];
        }
        dx = xmax - xmin;
        dy = ymax - ymin;
        dmax = Math.max(dx, dy);
        xmid = xmin + dx * 0.5;
        ymid = ymin + dy * 0.5;
        return [
            [xmid - 20 * dmax, ymid -      dmax],
            [xmid            , ymid + 20 * dmax],
            [xmid + 20 * dmax, ymid -      dmax]
        ];
    }

    function circumcircle(vertices, i, j, k) {
        var x1 = vertices[i][0],
            y1 = vertices[i][1],
            x2 = vertices[j][0],
            y2 = vertices[j][1],
            x3 = vertices[k][0],
            y3 = vertices[k][1],
            fabsy1y2 = Math.abs(y1 - y2),
            fabsy2y3 = Math.abs(y2 - y3),
            xc, yc, m1, m2, mx1, mx2, my1, my2, dx, dy;

        /* Check for coincident points */
        if(fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
            throw new Error("Eek! Coincident points!");

        if(fabsy1y2 < EPSILON) {
            m2  = -((x3 - x2) / (y3 - y2));
            mx2 = (x2 + x3) / 2.0;
            my2 = (y2 + y3) / 2.0;
            xc  = (x2 + x1) / 2.0;
            yc  = m2 * (xc - mx2) + my2;
        }

        else if(fabsy2y3 < EPSILON) {
            m1  = -((x2 - x1) / (y2 - y1));
            mx1 = (x1 + x2) / 2.0;
            my1 = (y1 + y2) / 2.0;
            xc  = (x3 + x2) / 2.0;
            yc  = m1 * (xc - mx1) + my1;
        }

        else {
            m1  = -((x2 - x1) / (y2 - y1));
            m2  = -((x3 - x2) / (y3 - y2));
            mx1 = (x1 + x2) / 2.0;
            mx2 = (x2 + x3) / 2.0;
            my1 = (y1 + y2) / 2.0;
            my2 = (y2 + y3) / 2.0;
            xc  = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
            yc  = (fabsy1y2 > fabsy2y3) ?
                m1 * (xc - mx1) + my1 :
                m2 * (xc - mx2) + my2;
        }

        dx = x2 - xc;
        dy = y2 - yc;
        return {i: i, j: j, k: k, x: xc, y: yc, r: dx * dx + dy * dy};
    }

    function dedup(edges) {
        var i, j, a, b, m, n;

        for(j = edges.length; j; ) {
            b = edges[--j];
            a = edges[--j];

            for(i = j; i; ) {
                n = edges[--i];
                m = edges[--i];

                if((a === m && b === n) || (a === n && b === m)) {
                    edges.splice(j, 2);
                    edges.splice(i, 2);
                    break;
                }
            }
        }
    }

    Delaunay = {
        triangulate: function(vertices, key) {
            var n = vertices.length,
                i, j, indices, st, open, closed, edges, dx, dy, a, b, c;

            /* Bail if there aren't enough vertices to form any triangles. */
            if(n < 3)
                return [];

            /* Slice out the actual vertices from the passed objects. (Duplicate the
             * array even if we don't, though, since we need to make a supertriangle
             * later on!) */
            vertices = vertices.slice(0);

            if(key)
                for(i = n; i--; )
                    vertices[i] = vertices[i][key];

            /* Make an array of indices into the vertex array, sorted by the
             * vertices' x-position. Force stable sorting by comparing indices if
             * the x-positions are equal. */
            indices = new Array(n);

            for(i = n; i--; )
                indices[i] = i;

            indices.sort(function(i, j) {
                var diff = vertices[j][0] - vertices[i][0];
                return diff !== 0 ? diff : i - j;
            });

            /* Next, find the vertices of the supertriangle (which contains all other
             * triangles), and append them onto the end of a (copy of) the vertex
             * array. */
            st = supertriangle(vertices);
            vertices.push(st[0], st[1], st[2]);

            /* Initialize the open list (containing the supertriangle and nothing
             * else) and the closed list (which is empty since we havn't processed
             * any triangles yet). */
            open   = [circumcircle(vertices, n + 0, n + 1, n + 2)];
            closed = [];
            edges  = [];

            /* Incrementally add each vertex to the mesh. */
            for(i = indices.length; i--; edges.length = 0) {
                c = indices[i];

                /* For each open triangle, check to see if the current point is
                 * inside it's circumcircle. If it is, remove the triangle and add
                 * it's edges to an edge list. */
                for(j = open.length; j--; ) {
                    /* If this point is to the right of this triangle's circumcircle,
                     * then this triangle should never get checked again. Remove it
                     * from the open list, add it to the closed list, and skip. */
                    dx = vertices[c][0] - open[j].x;
                    if(dx > 0.0 && dx * dx > open[j].r) {
                        closed.push(open[j]);
                        open.splice(j, 1);
                        continue;
                    }

                    /* If we're outside the circumcircle, skip this triangle. */
                    dy = vertices[c][1] - open[j].y;
                    if(dx * dx + dy * dy - open[j].r > EPSILON)
                        continue;

                    /* Remove the triangle and add it's edges to the edge list. */
                    edges.push(
                        open[j].i, open[j].j,
                        open[j].j, open[j].k,
                        open[j].k, open[j].i
                    );
                    open.splice(j, 1);
                }

                /* Remove any doubled edges. */
                dedup(edges);

                /* Add a new triangle for each edge. */
                for(j = edges.length; j; ) {
                    b = edges[--j];
                    a = edges[--j];
                    open.push(circumcircle(vertices, a, b, c));
                }
            }

            /* Copy any remaining open triangles to the closed list, and then
             * remove any triangles that share a vertex with the supertriangle,
             * building a list of triplets that represent triangles. */
            for(i = open.length; i--; )
                closed.push(open[i]);
            open.length = 0;

            for(i = closed.length; i--; )
                if(closed[i].i < n && closed[i].j < n && closed[i].k < n)
                    open.push(closed[i].i, closed[i].j, closed[i].k);

            /* Yay, we're done! */
            return open;
        },
        contains: function(tri, p) {
            /* Bounding box test first, for quick rejections. */
            if((p[0] < tri[0][0] && p[0] < tri[1][0] && p[0] < tri[2][0]) ||
                (p[0] > tri[0][0] && p[0] > tri[1][0] && p[0] > tri[2][0]) ||
                (p[1] < tri[0][1] && p[1] < tri[1][1] && p[1] < tri[2][1]) ||
                (p[1] > tri[0][1] && p[1] > tri[1][1] && p[1] > tri[2][1]))
                return null;

            var a = tri[1][0] - tri[0][0],
                b = tri[2][0] - tri[0][0],
                c = tri[1][1] - tri[0][1],
                d = tri[2][1] - tri[0][1],
                i = a * d - b * c;

            /* Degenerate tri. */
            if(i === 0.0)
                return null;

            var u = (d * (p[0] - tri[0][0]) - b * (p[1] - tri[0][1])) / i,
                v = (a * (p[1] - tri[0][1]) - c * (p[0] - tri[0][0])) / i;

            /* If we're outside the tri, fail. */
            if(u < 0.0 || v < 0.0 || (u + v) > 1.0)
                return null;

            return [u, v];
        }
    };

    if(typeof module !== "undefined")
        module.exports = Delaunay;
})();

在 Python 中，Triangle 库是一个用于进行二维三角剖分的第三方库。它是 Triangle 库的 Python 封装版本，Triangle 库本身是一个用 C 语言编写的高效三角剖分工具。Python 的 Triangle 库提供了 Triangle 库的功能，并将其包装成一个 Python 模块，方便在 Python 中调用。

以下是对 Triangle 库的简要解释以及一个使用案例：

功能特点：

Triangle 库用于进行二维三角剖分，它能够将给定的点集转换为三角形网格。
该库支持约束三角剖分、最大最小角限制、区域约束等功能。
Triangle 库提供了高效的算法和数据结构，能够处理大型点集和复杂的几何形状。

使用方法：

首先，需要安装 Triangle 库。你可以在 Python 中使用 pip 安装它：pip install triangle。
导入 triangle 模块：import triangle。
使用 triangle 模块提供的函数来执行三角剖分。最常用的函数是 triangle.triangulate()，它用于执行简单的点集三角剖分。你还可以使用其他函数来执行不同类型的三角剖分，例如约束三角剖分、最大最小角限制等。
一旦完成三角剖分，你可以使用返回的对象来访问三角剖分的结果，例如三角形的顶点、边、面等。

示例用途：

Triangle 库可以用于生成三角形网格，用于计算机图形学、地理信息系统、有限元分析等领域的模拟和可视化。
它还可以用于执行几何计算，例如计算三角形的面积、周长、重心等。
Triangle 库的约束三角剖分功能可以用于生成符合特定要求的三角形网格，例如在给定约束条件下进行网格划分。
import triangle
import matplotlib.pyplot as plt

# 创建一个随机点集
points = triangle.get_data('spiral')

# 执行三角剖分
tri = triangle.triangulate({'vertices': points})

# 绘制三角剖分结果
plt.triplot(points[:, 0], points[:, 1], tri['triangles'])
plt.plot(points[:, 0], points[:, 1], 'o')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Triangle Triangulation')
plt.show()

构网要素:顶点坐标、顶点颜色、索引（三角形由那几个点构建）

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

# 输入经纬度坐标和颜色集合
lon_lat = np.array([[lon1, lat1], [lon2, lat2], [lon3, lat3], ...])  # 输入经纬度坐标
colors = np.array([[r1, g1, b1, a1], [r2, g2, b2, a2], [r3, g3, b3, a3], ...])  # 输入颜色集合

# 使用Delaunay三角剖分算法获取三角形网格
tri = Delaunay(lon_lat)

# 创建绘图对象
fig, ax = plt.subplots()

# 绘制三角形并填充颜色
for vertices, color in zip(tri.simplices, colors):
    ax.fill(lon_lat[vertices, 0], lon_lat[vertices, 1], color=color)

# 隐藏坐标轴
ax.axis('off')

# 保存图像为PNG格式
plt.savefig('triangles.png', bbox_inches='tight', pad_inches=0, transparent=True)

# 显示图像（可选）
plt.show()