基于劳埃德算法、柏林噪声等的地图生成算法研究
灵感来自于一个突然的想法，我想实现一个可以生成虚构地图的程序。
（ps1，维基百科很好用啊。之前没怎么用过，现在用了发现它对算法的讲解不弱于一些专门讲算法的博客。）
（ps2，没有手动实现还是缺少点感觉，但是不用 VB 又太慢了。）
1. Voronoi 图（沃罗诺伊图）
在数学中，沃罗诺伊图是将平面划分为若干区域的一种方式，每个区域都接近一组给定对象中的某一个。它也可以被归类为一种密铺。在最简单的情况下，这些对象只是平面上的有限个点（称为种子、站点或生成元）。每个种子对应一个区域，称为沃罗诺伊细胞，包含平面上所有离该种子比离其他任何种子更近的点。一组点的沃罗诺伊图与该点集的德劳内三角剖分是对偶的。
简单来说，就是通过在一个度量空间里面生成一组随机点，再通过空间内点与中心点的距离来将空间进行切割的一个算法。这里使用它，主要是为了生成不规则多边形排列。
1.2. Fortune 算法（扫描线算法）
按照定义，我们可以很快给出一套求解方案，对于任意两个点作中垂线。因为中垂线上的点距离两个点相同，所以两个点之间的边界必定是中垂线上的线段、直线或射线。找到所有中垂线的集合，然后按照交点和方向进行划分，就可以简单地求解沃罗诺伊图。但是这样求解的时间复杂度至少为 \(O(n^2\log n)\)。所以在这一章，我将讲一个快速（并不）的求解算法，扫描线算法。
该算法基于抛物线设计。我们将引入一条通过平面移动的直线，而不是考虑不同地点之间的距离，并利用这条直线进行更有效的距离比较。我们称这条线为扫描线。
让我们首先记住，如果我们已知一个点 \(p\) 和一条直线 \(l\)（不包含 \(p\)），那么到 \(p\) 的距离与到 \(l\) 的距离相等的点，会构成一条抛物线。
我们将使用 \(P_{p,l}\) 来表示这条抛物线。和之前一样，是抛物线将平面划分成两个区域，一部分点更接近 \(p\)，另一部分点更接近 \(l\)。
考虑一个点 \(q\)，它的坐标是 \(q=(q_x,q_y)\)，它到 \(p\) 的距离表示为 \(d(q,p)\)。接下来，扫描线是一条水平线，我们称它的纵坐标为 \(l_y\)，所以 \(q\) 和 \(l\) 之间的距离是 \(q_y-l_y\)。因此，位于抛物线 \(P_{p,l}\) 上的点 \(q\) 具有如下关系。
\[ d(q,p)=q_y-l_y \]
更普遍来说，
\[ \begin{aligned} d(q,p) < q_y-l_y && \text{若 } q \text{ 在 } P_{p,l} \text{ 上方} \end{aligned} \]
\[ \begin{aligned} d(q,p) &= q_y-l_y && \text{若 } q \text{ 在 } P_{p,l} \text{ 上} \end{aligned} \]
\[ \begin{aligned} d(q,p) > q_y-l_y && \text{若 } q \text{ 在 } P_{p,l} \text{ 下方} \end{aligned} \]
现在，水平扫描线 \(l\) 将通过平面向下移动。在任何时候，我们都只考虑位于扫描线和由这些点定义的抛物线 \(P_{p,l}\) 之上的站点。如下图所示，因为抛物线上的点到线和点的距离相等，所以两条抛物线的交点到两个焦点（中心点）的距离相等。也就是说，该点位于边界（Voronoi 图的边缘）上。通过这个方法，就可以快速计算出边界。
现在我们可以介绍 \(Fortune's\ algorithm\) 中真正关键的概念，海滩线。它的定义是最低的抛物线弧所形成的曲线，也就是上图中黄色的曲线。随着扫描线不断下降，海滩线上每两条抛物线的交点，正是 Voronoi 图的边缘。这意味着，当扫描线沿平面向下移动时，交点将扫描 Voronoi 图的边缘。因此，要构建 Voronoi 图，我们只需要跟踪交点。
1.2.1. Finding the edges（找边）
当一个新点被添加进来时，它会在原本的海滩线上添加 1 个（焦点在扫描线上时，抛物线就是一条垂直于扫描线的直线）或 2 个交点；当出现两个交点时，它们会构造为一条线段，如下图所示。
在扫描过程中，如下图所示，它最终会收敛到边界或无限远处。
1.2.2. Finding the vertices（找点）
当扫描线向下移动时，交点沿着 Voronoi 图的边缘连续移动，直到到达图的顶点。当沿海滩线的抛物线消失时，也就是两条抛物线相交于一条抛物线的顶点时，就会发生这种情况。这个情况称为圆事件。
在海滩线上，会出现抛物线弧的消失。当一条抛物线收缩到一点 \(q\) 时，这个点就位于三条抛物线上，一条包含消失的弧线，另外两条包含两边的弧线。这意味着 \(q\) 到三个点的距离是相等的，而这三个点对应于这些抛物线，因此这三个点位于以 \(q\) 为圆心的圆上。因此，当扫描线经过那个圆的底部时，我们就找到了顶点。
请注意，圆事件还会创建一条新边。这与一个新的断点形成有关，这个断点来自两条抛物线弧的交点，而这两条抛物线弧正是与消失的弧相邻的弧。这条新边，是除消失的抛物线外，另外两个抛物线中心点的边界。
1.2.3. The algorithm（算法）
Fortune 算法通过扫描线从上到下移动，追踪海滩线上抛物线的顺序，来构建 Voronoi 图。
海滩线的顺序只会在两种事件发生时改变，站点事件（遇到新点）和圆事件（抛物线弧消失）。
遇到站点事件时，将新站点插入海滩线，并记录一条新边。
遇到圆事件时，记录一个顶点（两条边的交点），并记录由新断点产生的新边。
每次事件后，都要检查新增的三条相邻抛物线弧是否会产生未来的圆事件；同时，移除已经失效的圆事件。
最终，算法会逐步找到所有边界和顶点。算法时间复杂度为 \(O(n\log n)\)。
感兴趣可以去看一个 Bilibili 的示例视频，Fortune 算法示例视频。
2. Perlin noise（柏林噪声）
柏林噪声是一种梯度噪声。它和直接随机数最大的区别是直接随机数每个点之间几乎没有关系，而柏林噪声会让相邻位置的数值连续变化，所以它看起来不像白噪声那样碎，而更像自然界里的云、山脉、海岸线或者纹理变化。
柏林噪声的计算规则是先在网格点上生成一组随机梯度向量。然后对于空间里的任意一个点，找到它所在的网格单元，分别计算它到周围格点的偏移向量，再把偏移向量和格点上的梯度向量做点积。最后通过平滑插值，把这些点积结果混合起来。
假设当前点为 \(p=(x,y)\)，它所在格子的左下角为 \(c_{00}\)，该角点对应的梯度向量为 \(g_{00}\)。那么这个角点对当前点的影响可以写成。
\[ n_{00}=g_{00}\cdot(p-c_{00}) \]
这里的 \(p-c_{00}\) 是从角点指向当前点的偏移向量，点积则表示当前点在这个梯度方向上的“投影程度”。如果偏移方向和梯度方向比较接近，结果就偏大；如果方向相反，结果就偏小。
在二维情况下，一个点会受到所在网格四个角的影响，所以我们会得到四个点积结果。
\[ \begin{aligned} n_{00} &= g_{00}\cdot(p-c_{00}) \\ n_{10} &= g_{10}\cdot(p-c_{10}) \\ n_{01} &= g_{01}\cdot(p-c_{01}) \\ n_{11} &= g_{11}\cdot(p-c_{11}) \end{aligned} \]
接下来要做的就是插值。设当前点在格子内部的局部坐标为 \(u=x-\lfloor x \rfloor\)，\(v=y-\lfloor y \rfloor\)。先在水平方向插值。
\[ \begin{aligned} a &= lerp(n_{00}, n_{10}, f(u)) \\ b &= lerp(n_{01}, n_{11}, f(u)) \end{aligned} \]
然后再在垂直方向插值。
\[ value = lerp(a,b,f(v)) \]
其中 \(lerp\) 是线性插值。
\[ lerp(a,b,t)=a+t(b-a) \]
不过这里的 \(t\) 一般不会直接使用原始的 \(u\) 或 \(v\)，而是会先经过一个平滑函数。经典实现中常用的是。
\[ f(t)=6t^5-15t^4+10t^3 \]
这个函数的作用是让插值在网格边界附近更加平滑。如果直接线性插值，网格之间容易出现比较明显的人工痕迹；而使用这个平滑函数后，噪声变化会更自然。
在地图生成里，我主要把柏林噪声当作海拔函数使用。也就是说，对于每个 Voronoi 单元，可以取它的中心点坐标 \((x,y)\)，然后计算。
\[ height = noise(x,y) \]
得到的 \(height\) 就可以作为该地块的海拔。低于某个阈值的区域可以设为海洋，高一点的是平原，再高一点就是山地。这样 Voronoi 图负责提供不规则地块，柏林噪声负责提供连续的地形起伏。
但是单层柏林噪声的细节比较单一，所以实际使用时一般会叠加多层噪声。低频噪声控制大陆、大山脉这种大结构，高频噪声负责局部细节。这个叠加过程通常可以写成。
\[ FBM(x,y)=\sum_{i=0}^{k-1} a_i \cdot noise(2^i x,2^i y) \]
其中 \(2^i\) 控制频率，\(a_i\) 控制每一层噪声的权重。一般来说，频率越高，权重越小。这样生成出来的高度图就会同时拥有大尺度结构和小尺度细节。
2.1. 岛屿柏林噪声
单纯使用柏林噪声生成海拔时，会有一个很明显的问题，它只是生成连续的随机起伏，并不会天然知道地图边缘应该是海洋。所以如果直接把噪声值当作高度，就可能出现大陆一直延伸到地图边界的情况。对于世界地图来说这没什么问题，但如果我想生成一座岛，那么就需要额外给噪声加一个越靠近边缘越低的约束。
这个约束一般可以理解为一个遮罩函数。假设地图坐标已经被缩放到 \([-1,1]\) 范围内，也就是中心点附近为 \((0,0)\)，地图边缘接近 \(-1\) 或 \(1\)。我们可以构造一个函数，让它在中心区域接近 1，在边缘区域逐渐接近 0，然后用它去塑造柏林噪声的海拔形状（就是加个权）。
最简单的一种方式是使用距离平方。
\[ m(x,y)=1-(x^2+y^2) \]
这个函数的直觉很简单，离地图中心越远，\(x^2+y^2\) 越大，最终的 \(m(x,y)\) 就越小。所以它比较适合生成接近圆形的岛屿。中心区域保留较高海拔，边缘区域被压低，最后就会自然形成海洋。
如果想让岛屿更贴近方形地图边界，也可以使用 Square Bump。
\[ m(x,y)=(1-x^2)(1-y^2) \]
这个函数在 \(x\) 或 \(y\) 接近边界时都会趋近于 0，所以它会让地图四周都逐渐变低。相比距离平方，它生成的岛屿不会那么圆，而是更容易适配方形画布。
有了遮罩之后，最直接的做法就是把柏林噪声和遮罩相乘。
\[ height(x,y)=noise(x,y)\cdot m(x,y) \]
这样中心区域的噪声基本保留，越靠近边缘，噪声的高度就越被压低。最后再设置一个海平面阈值，例如。
\[ \begin{aligned} height(x,y) < seaLevel && \text{海洋} \\ height(x,y) \ge seaLevel && \text{陆地} \end{aligned} \]
除了这两种基础函数，还可以使用一些形状更柔和的遮罩。例如双曲面形式。
\[ m(x,y)=1-\sqrt{x^2+y^2+c^2} \]
这里的 \(c\) 可以理解为一个平滑参数。普通的距离函数在中心附近可能会比较尖，而加入 \(c^2\) 后，中心区域会更圆滑一点。如果希望中心为 1，边缘中点为 0，还可以对它做一次归一化。
\[ m(x,y)=1-\frac{\sqrt{x^2+y^2+c^2}-c}{\sqrt{1+c^2}-c} \]
还有一种比较有意思的是三角函数乘积。
\[ m(x,y)=\cos\left(\frac{\pi x}{2}\right)\cos\left(\frac{\pi y}{2}\right) \]
它和 Square Bump 有点类似，也会让边界逐渐趋近于 0，但是下降方式更平滑一些。对于想要比较柔和的海岸过渡时，这种函数会更好。
如果想在圆形和方形之间取一个折中，也可以使用 Squircle 形式。
\[ m(x,y)=1-\sqrt{x^4+y^4} \]
它生成的形状不像距离平方那样偏圆，也不像 Square Bump 那样明显贴近方形，而是介于两者之间。对于岛屿生成来说，这种遮罩有时候会比纯圆形更自然，因为真实岛屿本来也不会是一个完美圆。
在实际实现里，我更倾向于把遮罩当作生成倾向，而不是一个绝对规则。也就是说，柏林噪声仍然负责提供地形起伏，遮罩只负责告诉它，越靠近边缘，越应该变成海洋。这样最后得到的岛屿既能保证边界大概率是海，又不会完全失去噪声带来的随机海岸线。
实际实现时存在一定小问题，比如某些遮罩函数会出现负值，那么直接相乘可能会把高度反向拉出奇怪的结果。所以在乘之前，通常需要先做一次截断或者归一化。
最终的岛屿高度可以写成。
\[ height(x,y)=FBM(x,y)\cdot \max(0,m(x,y)) \]
目前的工作做到让柏林噪声负责高低，岛屿遮罩负责大地形的渲染。两者叠在一起，才能生成一张像岛屿的地图，而不是一整块无限延伸的随机大陆。
3. 河流生成
因为比较懒，我的生成逻辑是随机选择海拔高点，然后开始广度搜索，并且通过一个目标函数对结果进行剪枝，最后以海洋或者湖结束。
真实一点的河流生成需要去模拟水流侵蚀、降雨量、汇水面积这些很复杂的东西（地形生成除了柏林噪声其实还有侵蚀算法，感兴趣可以去看 Josh's Channel 的视频）。比较复杂的算法比如先生成降雨，再算每个区域的流向和累积流量，最后只有流量足够大的地方才画成河流。但是我这里的目标只是让地图上出现一些看起来还算合理的河道，所以实现上更像是一个路径搜索问题，从高处出发，沿着比较合适的方向走，直到碰到海洋或者已有河流。
3.1. 为什么河流要沿着 Voronoi 边界走
在这套地图生成里，地形的基础结构是 Voronoi 单元。如果直接在像素级高度图上生成河流，河道虽然会比较自由，但是最后叠到多边形地图上时，很容易出现河流穿过地块中间、和边界系统不一致的问题。所以我的实现选择了一个取巧但很稳定的办法，只允许河流在 Voronoi 单元的边界上移动。
代码里会先根据 \(labels\) 找出所有边界像素。判断方式很简单，如果一个像素和左边或上边的像素属于不同的 Voronoi 单元，那么它就是边界的一部分。
\[ boundary(r,c)= \begin{cases} 1, & labels(r,c)\ne labels(r,c-1)\ \text{或}\ labels(r,c)\ne labels(r-1,c)\\ 0, & \text{否则} \end{cases} \]
然后把这些边界像素当作图上的节点。每个节点记录自己的坐标、海拔 \(elevation\)，以及是否为水域 \(water\)。节点之间使用 8 邻接连接，也就是上下左右和四个斜方向都可以走。这样一来，河流生成就从“在平面上画线”变成了“在边界图上找一条路径”。
3.2. 随机选择河流源头
河流源头不能随便选，否则很容易从海边或者平原开始流，不符合常理，一般而言河流是由于高山上的积雪融化导致的，所以起始点应该较高。代码里的做法是先筛选出所有满足条件的边界点。
\[ elevation(p)\ge h_{min}\quad \text{且}\quad water(p)=false \]
在我的实现中，推荐版本里 \(h_{min}\) 大概取 \(0.58\)，也就是只从比较高的陆地区域开始选。然后再把这些候选点随机打乱，并且要求源头之间保持一定距离，例如 \(28\) 个像素左右。这个限制主要是为了避免所有河流都挤在同一片山地里，看起来像几条线从一个点炸开。
所以源头选择并不是完全随机，而是“先高海拔筛选，再随机抽取，再用距离约束去重”。这样既保留了一点随机性，又不会让结果完全失控。
3.3. 从贪心下降到 Beam Search
最直观的办法是贪心，每一步都从邻居里选一个更低的点，直到走到海洋。这个方法很好写。但是问题也很明显，如果前面一步选错了，后面就容易陷入局部最低点。遇到平台、局部洼地、或者边界走势比较绕的时候（想象一条盘山公路），河流可能很快卡死。
所以后面的版本我改成了 Beam Search。它就是一个带有剪枝的广度搜索，每一轮不是只保留一条路径，而是保留一批当前看起来比较好的路径。下一轮再从这些路径继续扩展，然后继续剪枝，只留下得分最高的前 \(beamWidth\) 条。
简单写成伪代码就是。
\[ beam_0=\{[source]\} \]
\[ beam_{t+1}=TopK\left(Expand(beam_t),\ beamWidth\right) \]
这里的 \(Expand\) 会把当前路径末端的所有可行邻居都加入进来。所谓“可行”，主要是指不能回到已经访问过的点，并且整体上要往低处走。不过为了处理平台和局部小坑，代码里允许一点点上浮，条件类似。
\[ elevation(next)\le elevation(current)+\epsilon \]
其中 \(\epsilon\) 就是 \(downhill\_tolerance\)。推荐版本里它大概是 \(0.02\)。这个值太小河流容易卡死，太大河流就能够爬山。所以它本质上就是一个松驰用的变量，防止边界约束过于严格导致难以求解。
3.4. 路径评分函数
Beam Search 的关键不在搜索本身，而在于怎么判断一条候选路径好不好，及量化指标。比如下降和长度的同时优化。这本质是多目标规划问题，所以代码里给每条候选路径设计了一个目标函数。
搜索过程中使用的是一个中间评分，大致考虑四件事，从源头到当前点下降了多少、路径走了多远、路径本身有多长，以及有没有接触到已有河流。
\[ score_{partial} = w_d\cdot drop +w_s\cdot distance +w_l\cdot length +w_m\cdot merge \]
其中 \(drop\) 表示源头海拔减去当前海拔，\(distance\) 同时参考路径长度和离源头的直线距离，\(merge\) 表示是否已经接入已有河道。这样评分的结果是，路径既倾向于往低处走，也不会太短，还允许后生成的河流汇入前面的河流。
当一条路径成功到达海洋，或者接触到已有河流时，它就会进入候选结果。最后再用完整路径评分函数重新筛选一次。
\[ score = 4.0\cdot totalDrop +0.035\cdot pathLength +1.2\cdot straightness +3.0\cdot \max(avgSlope,0) +mergeBonus \]
这里我额外加了 \(straightness\)，不是为了让河流变成直线，而是为了避免它在很小的区域里绕圈。它的定义大致是。
\[ straightness=\frac{directDistance}{pathLength} \]
如果一条河流走了很长，但是离源头并不远，说明它可能在原地绕了很多圈，这种路径就不应该得太高分。相反，一条河流可以弯曲，但是总体方向最好还是能从山地走向低处。
3.5. 河流终止和合流
河流的结束条件有两个，一个是到达海洋，另一个是接入已有河流。到达海洋很好理解，接入已有河流则是为了形成河网。如果每条河都必须单独流到海里，地图上会出现很多彼此平行的小线，但是现实里的河流经常会合流，所以允许后生成的河流接到前面的河道上，效果会自然很多。
代码里用一个 \(riverMask\) 记录已经生成过的河流像素。每生成一条河，就把它走过的点写入这个 mask。后面的搜索如果碰到了 \(riverMask=true\) 的节点，就可以把这条路径视为成功，并且在评分里给一个合流奖励。
最后，只有评分超过阈值的路径才会真正加入结果。这样可以过滤掉那种虽然找到了终点，但是太短、太绕、下降不明显的河流。整个流程大概可以总结为。
\[ \text{高海拔源头} \rightarrow \text{边界图搜索} \rightarrow \text{目标函数剪枝} \rightarrow \text{海洋或合流} \]
所以这套河流生成并不是严格的水文学模拟，而是一个比较工程化的近似，Voronoi 边界决定河道能走哪里，柏林噪声高度决定河道倾向往哪里走，Beam Search 和目标函数决定最后保留哪几条河。优点是可控、成功率高，而且和前面的多边形地图结构比较搭；缺点也很明显，它不会真的计算流域和水量，所以河流的粗细、支流数量这些还需要后续额外处理。
4. 噪声边缘
到这里为止，地图已经有了 Voronoi 单元、岛屿高度、海岸线和河流。但是还有一个很明显的问题，Voronoi 图的边界太直了。直线边界在调试的时候很舒服，因为每一条边都很清楚；但是放到地图里就会有一种“这个世界是被尺子切出来的”感觉，尤其是海岸线、山脉分界、区域边界这种本来应该比较自然的地方，看起来会有点不自然。
Red Blob Games 这篇文章里的思路很适合解决这个问题，不要直接移动整个 Voronoi 单元，也不要重新生成一套复杂曲线，而是只把每一条 Voronoi 边替换成一条带噪声的折线。这样单元之间的拓扑关系基本不变，但是视觉上直线边界会变成弯曲边界。
简单来说，就是原本一条边只有两个端点。
\[ edge=(corner_A,corner_B) \]
现在把它变成一个点序列。
\[ noisyEdge=\{p_0,p_1,p_2,\dots,p_n\} \]
其中 \(p_0=corner_A\)，\(p_n=corner_B\)，中间的点通过随机扰动生成。最后渲染时不再画原来的直线，而是依次连接这些点。
4.1. 四边形约束
如果只是随便把边的中点往旁边挪，确实也能得到弯曲边界，但是很容易出问题，线可能穿进别的单元，或者两条边互相交叉。所以 Red Blob 的做法不是在无限平面里扰动，而是先给每条边找一个局部四边形。
对一条 Voronoi 边来说，它连接两个角点 \(corner_A\) 和 \(corner_B\)，同时它两侧有两个生成点，也就是 \(center_A\) 和 \(center_B\)。这四个点可以围出一个局部区域。
\[ Q=(center_A,corner_A,center_B,corner_B) \]
这相当于告诉算法，这条边可以在这个四边形附近抖动，但是不要跑得太离谱。我的实现里也是这个思路，`generate_edge` 接收的参数就是 \(center_A,center_B,corner_A,corner_B\)。
下面这张是我在代码里调试四边形约束时输出的图。黑色线条表示每条 Voronoi 边对应的局部四边形结构，红点和蓝点分别对应参与构造的 site（中心点）和 vertex（边界点）。通过这四个点构造的平行四边形就可以将每条边界的波动范围单独约束出来。
4.2. 递归中点扰动
生成噪声边缘的核心是递归细分。最开始只有一条直线，从 \(corner_A\) 到 \(corner_B\)。然后取这条线段的中点。
\[ mid=\frac{corner_A+corner_B}{2} \]
接下来沿着两个中心点的方向做一个随机偏移。代码里的方向是。
\[ dir=\frac{center_B-center_A}{\left\|center_B-center_A\right\|} \]
然后给中点加上一个随机量。
\[ mid'=mid+dir\cdot random(-offset,offset) \]
其中 \(offset\) 和两个中心点之间的距离有关，代码里大致是。
\[ offsetLimit=amplitude\cdot \left\|center_B-center_A\right\| \]
同时还会用 \(maxOffset\) 限制最大偏移，避免某些很长的边被抖得太夸张。我的主程序里相关参数大概是 \(noiseIterations=10\)，\(amplitude=0.25\)，\(maxOffset=0.8\)。也就是说，它会递归很多层，但是每次偏移都被限制住。
递归的过程可以理解为，先把一条线切成两段，然后对左右两段继续做同样的操作。为了让细节逐渐变小，代码里每深入一层，\(amplitude\) 会减半。
更具体一点看，每次生成新的扰动中点后，原来的四边形会被拆成左右两个更小的局部四边形。下一轮递归并不是还在原来的大四边形里乱偏移，而是在各自的小四边形里继续做同样的事。这样 noisy edge 的形状虽然是随机的，但是仍然被局部结构约束住。
\[ amplitude_{k+1}=\frac{amplitude_k}{2} \]
所以大的弯曲在前几层出现，小的毛边在后几层出现。这样生成出来的线不会像纯随机游走那样抖成一团，而是有一种从大形状到小细节逐渐叠加的感觉。
4.3. 从 Voronoi 边到 noisy edge
具体实现时，`NoisyEdgeGenerator.generate_from_edges` 会先根据生成点重新计算一次 Voronoi 图，然后遍历所有 ridge。这里需要跳过无限边，因为无限边没有两个稳定的 Voronoi 顶点。
\[ -1\in ridgeVertices \Rightarrow skip \]
对剩下的有限边，代码会拿到两个 site 和两个 vertex，也就是前面说的两个中心点和两个角点。然后再判断这条边到底要不要生成 noisy edge。我的实现里有几种入口，可以传入显式的 \(boundaryPairs\)，也可以传入 \(boundarySiteLabels\)。如果一条边两侧的 label 不同，就说明它是区域边界，需要被画出来。
\[ label(site_A)\ne label(site_B) \]
这一步很重要。因为我不一定想把所有背景小单元的边都画出来，有时候只想画区域和区域之间的边界，或者只处理海岸线附近的边界。用 label 来筛选，就可以把“是否需要噪声边缘”和“怎么生成噪声边缘”拆开。
下面这张图左边是普通 Voronoi 直边，右边是同一组点生成的 noisy edge。可以看到拓扑结构还在，但视觉上已经不再那么像规则切割出来的网格了。
如果把区域边界也叠上去，就会更清楚一点。下面这张图里，普通的小 Voronoi 单元仍然存在，但是更粗的边界已经开始表现出“哪些边真的会作为区域边界被处理”。这也是我后面做地形区域时比较依赖的一步，底层可以有很多小格子，但是最终展示出来的边界不一定要把所有小格子都暴露出来。
4.4. 重新填充区域
只把线画弯还不够。因为原来的颜色区域是按照直线 Voronoi 边界填充的，如果直接把弯线画在上面，就会出现“线是弯的，颜色块还是直的”的问题。所以我的代码里在生成 noisy edge mask 之后，又做了一次区域重填充。
大概流程是，先把所有 noisy edge 画到一张黑白 mask 上，然后对 \(\neg noisyEdgeMask\) 做连通域标记。每一个连通块就是被弯曲边界切出来的新区域。
\[ components=LabelConnectedComponents(\neg noisyEdgeMask) \]
接下来，代码会看这个连通块内部原来主要属于哪个 site label，然后用这个 label 对应的颜色重新填充整块区域。这样做之后，颜色边界和 noisy edge 就能对齐了。海岸线也类似，会额外生成 \(noisyCoastMask\)，再根据连通区域是否接触地图边缘、内部水域比例等信息，决定哪些地方应该继续作为海洋，哪些地方可以保留为沙滩或陆地。
所以噪声边缘这一步并不是简单地“在图上画一条抖动线”，而是包含了三个阶段，先生成弯曲边界，再根据弯曲边界切分区域，最后重新给区域上色。前两步决定边界形状，最后一步决定它能不能真正融入地图。
下面几张就是把这套逻辑放进完整地图生成流程后的效果。第一张保留了比较多的区域结构，能看出 noisy edge 对海岸线和内部分区的影响；第二张更接近最终的分区上色结果；第三张则把海拔和水域效果叠进去，用来观察噪声边缘在真实地形图里的表现。
\[ VoronoiEdge \rightarrow NoisyPolyline \rightarrow EdgeMask \rightarrow RefillRegions \]
总的来说，这个算法非常适合我的地图生成流程。Voronoi 负责给出稳定的区域拓扑，noisy edge 负责把这些区域的边界变得自然一点。它不改变“谁和谁相邻”这个大结构，只改变“它们之间的边界长什么样”。这种分层处理对程序来说比较好控制，对视觉效果来说也足够明显。
5. 代码实现中的技巧
前面讲的大多是地图长什么样，这一节稍微补一下代码实现里比较实用的两个小技巧。用于加速程序。
5.1. 泊松圆盘算法
最开始生成 Voronoi 图时，最直接的办法当然是随机撒点。比如直接让 \(x\) 和 \(y\) 在地图范围内均匀随机。
\[ p_i=(random(x_{min},x_{max}),random(y_{min},y_{max})) \]
这样写很简单，但是视觉效果经常不太稳定。有些点会挤在一起，有些地方又会空一大片。放到 Voronoi 图里，就会出现一些特别小的碎片单元和特别大的空洞单元。通过这个方法实现的 Voronoi 图分布不均，会让后面的 Lloyd 松弛，依靠更大的迭代次数才收敛。
所以我在新版实现里使用了泊松圆盘采样。它的核心约束很简单，任意两个采样点之间至少保持一个最小距离 \(r\)。
\[ \forall i\ne j,\quad \left\|p_i-p_j\right\|\ge r \]
这个约束带来的效果就是，点仍然是随机的，但是不会随机到挤成一起。它有点像一种比较自然的均匀随机。
具体到我的代码里，`PoissonDiscSampler` 使用的是 `scipy.stats.qmc.PoissonDisk`。代码会先把地图范围缩放到单位正方形，按半径生成一批点，再从里面取需要的数量，最后缩放回真实地图坐标。这里需要注意一点，半径不是越大越好。半径太大时，空间里能塞下的点数量就会变少。
\[ unitRadius=\frac{radius}{min(width,height)} \]
5.2. k-d Tree
k-d Tree 主要解决的是“最近邻查询”问题。Voronoi 图本质上就离不开最近邻，一个像素到底属于哪个 site，就是看它离哪个点最近。
\[ label(q)=argmin_i\left\|q-p_i\right\| \]
如果直接暴力计算，那么每一个像素都要和所有 site 算一遍距离。假设有 \(m\) 个像素、\(n\) 个点，复杂度大概就是。
\[ O(mn) \]
对小图来说这还可以忍，但是我的地图会做像素级归属、Lloyd 近似质心、背景单元到大区域的映射，这些地方都会反复问“离我最近的点是谁”。如果一直暴力算，程序会很快变得不好跑。
k-d Tree 的思路是先把点按坐标轴递归切分成一棵树。查询时，不需要把所有点都看一遍，而是优先进入最可能包含最近点的空间分支，再根据当前最优距离去剪掉一些不可能更近的分支。下面这张图里，灰色和紫色的线就是空间被递归切开的结果，红点是查询点，蓝点是找到的最近点。
例如 Lloyd 松弛时，我会先在地图上铺一层采样网格，然后用 k-d Tree 给每个采样点找最近的 site。找完以后，属于同一个 site 的采样点求平均值，就可以近似得到这个 Voronoi cell 的质心。
\[ centroid_i\approx \frac{1}{|S_i|} \sum_{q\in S_i}q \]
其中 \(S_i\) 表示所有被分配给第 \(i\) 个 site 的采样点集合。也就是说，k-d Tree 在这里不是直接改变地图形状，而是让“把大量像素/采样点分配给最近 site”这件事变得足够快。
小结
本篇文章就先到这里了，因为图片的装备、公式的书写、资料的收集、代码的复现其实挺麻烦的，我花了两周才书写好这篇文章。其实这篇文章有几个部分并没有讲完，一个是区域吞并，一个是大区域划分，以及河流的流量。有关于联通图和聚类分析，展开来就很难讲完了，所以停在这里。如果后面有空，我可能还会尝试大陆的碰撞、人口分布、森林环境、湿度等等的模拟。