Closed-Form Construction of Voronoi Diagrams with Star-Shaped Metrics

概述

这篇工作研究如何在平面与三维曲面上生成丰富、可平滑渐变的胞状图案（cellular patterns）。这类图案在装饰纹样、建筑表面和机械超材料中都很常见，兼具美观与功能。传统的同胚镶嵌（isohedral tilings）虽然规整对称，但缺乏灵活性：不同镶嵌之间一般不存在平滑插值，也无法在双曲率曲面上无畸变地铺展。

Voronoi 图提供了更灵活的图案化手段。给定一组站点（sites），Voronoi 图由到最近站点等距的点集定义。将标准欧氏距离替换为广义距离，尤其是星形（star-shaped）度量，可以极大扩展可能的胞形空间，并通过插值度量参数实现图案的连续渐变。此前 Martínez 等人用基于栅格化（rasterization）的前沿推进算法在二维上探索了这一思路，效果可观，但其离散本质无法提供梯度，限制了对图案质量的优化控制。

本文提出一个全新的、闭式（closed-form）、完全可微的构造方法，针对分段线性星形度量的 Voronoi 图。它支持对站点位置和度量参数进行基于梯度的优化，以满足美学与功能目标，并可自然推广到任意维度，包括三维曲面。

星形度量与广义 Voronoi 图

站点 \(c_i\) 的 Voronoi 胞定义为

\[R_i = \{p : d_i(p) < d_j(p)\ \forall j \neq i\}\]

其中 \(d_i(p) = \lVert p - c_i \rVert\) 是点 \(p\) 到站点 \(c_i\) 的欧氏距离。用其他距离度量替换即得到广义 Voronoi 图。

星形距离由一个包含原点的星形集合 \(M_i \subset \mathbb{R}^n\) 定义。当 \(M_i\) 边界上任意一点 \(p\) 到原点的线段完全落在集合内部时，\(M_i\) 就是星形的。\(M_i\) 类比于欧氏距离下的单位球：位于 \(c_i + \partial M_i\) 上的点到站点的距离为 1，位于 \(c_i + \alpha \partial M_i\) 上的点距离为 \(\alpha\)。当度量集合是多边形或多面体时，得到的就是分段线性星形距离。

闭式构造的核心思想

算法的输入包括：每个 Voronoi 站点的位置 \(c_i\)；对应的度量多面体 \(\partial M_i\)（由若干度量顶点定义，三维时还需其三角化，且必须是星形的）；以及定义域（二维多边形域或三维曲面网格，允许细长三角形、孔洞、非流形与不连通网格）。

关键观察是：站点 \(c_i\) 周围的空间可以被划分为一组广义锥（generalized cones）\(K_i^j\)。锥 \(K_i^j\) 是从 \(c_i\) 出发、穿过度量多面体某个面 \(F_i^j\) 的所有射线之并。锥内的星形距离是线性的：

\[d_i^j(p) = \frac{n_i^j \cdot (p - c_i)}{b_i^j}\]

其中 \(n_i^j \cdot x + b_i^j = 0\) 是包含 \(F_i^j\) 的超平面。由此，两站点 \(c_0\) 与 \(c_1\) 在锥交集 \(K_0^j \cap K_1^k\) 上的二分线（bisector）也是线性的，令 \(d_0^j(x) = d_1^k(x)\) 得

\[\left(b_1^k n_0^j - b_0^j n_1^k\right) \cdot x + \left(b_0^j n_1^k \cdot c_1 - b_1^k n_0^j \cdot c_0\right) = 0\]

因此完整的二分线是分段线性的，其折点只出现在锥的过渡处（距离在此仅为 \(C^0\) 连续）。当存在第三个站点更近时，相应的二分线段被排除在图之外；若无线段留存，则两站点在 Voronoi 图中不相邻。

Voronoi 顶点与可微性

Voronoi 图的几何完全由其顶点及连接边确定。在 \(m\) 维空间中，任一 Voronoi 顶点都是 \(m\) 个超平面的交点，这些平面可以是 Voronoi 二分线、锥分隔面，或边界面。确定方程组后，顶点位置以闭式求解，相当于求 2×2 或 3×3 矩阵的逆；顶点位置及其对站点位置和度量顶点的导数，均借助计算机代数系统得到表达式。

在二维与曲面受限情形下，作者归纳出五类顶点：二分线交汇点（Bisector Junction）、锥过渡点（Cone Transition）、边界二分点（Boundary Bisector）、边界锥过渡点（Boundary Cone Transition）以及边界角点（Boundary Corner）。通过二维中每个顶点都是两条线（Voronoi 二分线、锥分隔线或边界边）交点这一事实，可以证明这组类型是完备的——唯一的例外是两条锥分隔线相交，因为锥分隔线本身不产生 Voronoi 边。曲面受限情形下这套类型同样适用，因为该构造可分解为逐三角形的二维构造。

由于顶点位置完全定义了图的几何，构造的可微性等价于顶点位置关于输入的可微性。单个顶点的导数总能计算，但在拓扑变化（某些顶点消失并被新顶点替代）处，构造一般不可微。不过某些定义在 Voronoi 图上的量（如胞体积）即便跨越拓扑变化也是可微的。

追踪式构造算法

算法通过遍历已知 Voronoi 边来增量构建整个图。维护一个「追踪」（trace）队列，每个追踪代表一条尚不知另一端点的顶点-边对。对每个追踪，沿边搜索最近的另一端点，若该顶点未曾访问则从中派生新追踪。给定合适的初始化，当整个定义域遍历完毕时算法终止。

初始化上，二维多边形域的完整边界必然存在于图中，故用边界顶点沿边界边的追踪初始化队列；三维曲面则把所有网格边视为边界边，先将全部网格顶点与边加入，再在后处理中移除（非闭合网格的外边保留），从而保证在拓扑不连通时仍覆盖整张网格。

追踪时，把被追踪边写成 \(p + t d,\ t > 0\)，在候选端点中寻找使 \(t\) 最小且为正的端点。边界追踪和二分线追踪各有几类候选端点，分别对应上面的顶点类型。总的追踪数（即 Voronoi 边数）与站点数和度量多边形的度成正比，瓶颈情形需遍历所有其他站点及其锥，故整体运行时间关于站点数和度量顶点数是二次的；实践中用锥的包围体层次结构来缩减候选端点集合。

边的方向由边界边或二分线方向确定，但只到符号。多数情况下可依据边界方向、指向锥内部或指向域内的要求来定号；否则通过比较决定追踪起点的锥来选择方向，例如从边界二分点出发的边界追踪，正确方向满足

\[\frac{d \cdot n_i^p}{\lvert b_i^p \rvert} < \frac{d \cdot n_j^q}{\lvert b_j^q \rvert}\]

去除不连通分量

星形距离 Voronoi 图有个特性：一个胞可能由多个不连通分量组成，这常带来独特而美观的镶嵌，但用于制造的结构不能有断开的曲线。作者提出一个精确算法，从既有 Voronoi 图中逐个移除不连通分量（即不包含自身站点的分量），每次移除后沿此前连向被删分量的边重新局部构造，且追踪时只考虑与被删分量共享过边的站点之间的二分线。

移除一个分量可能生成新的不连通分量，但作者证明该递归过程一定终止：新分量的邻居数严格少于被删分量，且可能顶点数有限（每个顶点由三个锥与边界边的组合唯一确定），故可能的分量数有限。实践中很少产生新分量，去除过程远快于初始构造。需要注意的是，结果镶嵌相对星形距离不再是严格的 Voronoi 图，但所有顶点仍属前述类型，且构造保持可微。

三维曲面上的图案生成

在三维曲面上生成准均匀镶嵌需要精心构造度量多面体。作者对所有三维结果使用挤出度量（extruded metrics），即把二维星形多边形沿曲面法向挤出，这比固定朝向的多面体度量能得到更均匀的镶嵌（例如在圆柱面上，固定度量会因朝向失配而使胞畸变）。

多数规则镶嵌（含三角与方形格）在多数三维曲面上无法无奇点存在。由 Euler 特征 \(\chi = V - E + F\) 可知，三角格有 \(E = 3V,\ F = 2V\)，方形格有 \(E = 2V,\ F = V\)，二者都要求 \(\chi = 0\)，仅对亏格为 1 的曲面（如环面）成立；球面等亏格 0 曲面 \(\chi = 2\)，故这些规则镶嵌必然带奇点。为得到近乎规则的站点分布，作者用条纹图案方法：将由最平滑方向场生成的条纹图案与其旋转 \(\frac{\pi}{2}\)（得矩形格）或 \(\frac{\pi}{3}\)（得三角格）的版本相交，交点即站点；三角格分布还可进一步用受限中心 Voronoi 镶嵌（CVT）优化改善。

度量对齐至关重要：在保持站点位置不变的情况下旋转度量多边形往往产生完全不同的镶嵌。作者提出逐扇区（per-sector）参数化，把每个站点度量多边形的扇区映射到实际邻居方向，从而在非规则邻域中生成均匀胞形。对于扇区各异的度量（如三折对称产生的「球-窝」图案），还需保证相邻性约束（若 \(c_i\) 的 A 扇区指向 \(c_j\)，则 \(c_j\) 的 B 扇区应指向 \(c_i\)），这些约束一般无法完全满足，作者将其建模为一个混合整数规划来尽量满足：

\[y^* = \arg\min_y \left( c_1 \sum \mathbb{1}_{y_{ij} = y_{ji}} + c_2 \sum \mathbb{1}_{y_i^k = y_i^{k+1}} \right),\quad y \in [A, B]\]

机械超材料优化

可微构造使得能对镶嵌做基于梯度的高层目标优化。作者在 Voronoi 图的边上生成杆网络（rod network）构造结构化薄片材料，施加单轴拉伸力求解变形平衡态，从沿拉伸方向的变形计算刚度、从垂直方向计算泊松比。沿用 Numerow 等人的平衡约束优化流程，可以优化目标泊松比或一系列拉伸方向上的目标刚度分布。

一个实验中，固定站点、优化三折对称度量多边形的参数，使泊松比从初始的 \(\nu = 0.332\) 降到 \(\nu = -0.392\)，优化后的图案在拉伸平衡态呈现拉胀（auxetic）材料特有的外鼓行为，并做了 3D 打印样件。另一实验中同时优化度量多边形与方形单元内的站点位置，以逼近极坐标图中给定的目标刚度分布。求平衡态用 Newton 法，度量多边形的平衡约束优化用全局收敛的移动渐近线法（MMA）。

实现与性能

方法用 C++ 实现，基于 CGAL 库，配合 Eigen、CHOLMOD 求解器、Gurobi（混合整数规划）、NLopt（MMA）、Geometry Central 与 Geogram（条纹图案与 CVT）以及 Polyscope（可视化）。代码开源在项目仓库。性能测试在一台 Mac mini（64 GB 内存、ARM CPU）上进行：例如球面上数百个站点、上万面片规模的图，构造时间在数百毫秒到数秒量级；更复杂的 Fertility 模型（超过十万面片）构造时间达数十秒；去除不连通分量通常只需数毫秒。

局限与展望

当前实现中，站点非一般位置（如严格规则格）时会出现数值问题，作者用相对格距 \(10^{-5}\) 的微小随机扰动即可完全避免；此外，多于三条 Voronoi 边交汇的角点等特殊情形仍需进一步处理。曲面上的图案化质量可借助最小化奇点的重网格化算法改进；对 CVT 站点放置的依赖使得非三角格（矩形、蜂窝）无法用于一般曲面，限制了可能镶嵌的空间；基于测地距离而非欧氏距离的内蕴表述有望减少高曲率区域的视觉瑕疵。机械优化目前仅限二维薄片材料，可推广到曲面上的结构网络；以美学为导向的目标（如目标胞形）也值得探索，但星形度量下更频繁的拓扑变化使该问题尤为困难。一个完整的三维版本及其在体积超材料上的应用在配套工作中另行描述。