泊松分布采样算法（Fast Poisson Disc Sampling）

前(fei)言(hua)

最近在看一些随机地图生成算法，涉及到生成Voronoi图，这需要提前在一个平面内随机生成一堆的点，这些点还要满足随机而且尽量平均分布在平面上。一般文章都提到采用Lloyd Relaxation算法，不过这个算法比较复杂，消耗也比较大，后来看到这个泊松分布采样算法，也是用于生成一堆平均分布的点的，而且算法复杂度在$O(N)$O(N)。
算法的说明论文可以参考这里。个人觉得描述得很不(fan)清(ren)晰(lei)，后来我看了别人的视频说明才真正了解这个算法，其实算法本身很简单，也很直观，下面会讲解下。完整代码我放在个人Github上。

算法说明

下面以2维平面为例：
假设我们需要在一个宽高为$(width,height)$(width,height)的平面内平均生成一堆的点，且这些点之间的距离不能小于$r$r。
我们可以先从平面内随机选一个点，然后在这个点附近随机找一些点，并判断这些点是否合法，合法的话则在这些点附近继续随机寻找，直到找不到合法点为止。

算法简单说起来就是这样，下面详细说下细节。

1. 为了保证能尽量填满整个平面， 随机找点时，采用与中心点距离为$[r,2r)$[r,2r)的圆环内找，这个距离能保证找的点距离自身大于$r$r，且不会离得太远，能填满整个平面。
2. 怎么判断一个点的附近已经找不到合法点了？算法定义了一个常量$k$k，对于每个点，我们尝试在它附近随机找$k$k次，如果都找不到，那么就认为这个点附近已经没有合法点。
3. 怎么快速判定随机找出的点是否合法？这个是算法的关键，可以采用一些空间划分方法来做（游戏场景也会经常用到），首先将平面划分成$m$m行$n$n列的格子，每个格子都保存了格子内部的点。这样当我需要判断一个点是否合法时，我只要和附近的格子内的点做判断即可。
4. 那怎么确定每个格子的大小？我们尽量让每个格子内部最多只能有1个点，这样数据结构就会简单很多。怎么做到呢？我们假设每个格子都是正方形，那正方形内部距离最远的点就是对角线的2个点，所以我们只要保证正方形的对角线长度大于等于$r$r，则正方形内部任意2个点之间的距离肯定小于$r$r，从而保证每个正方形内部肯定最多只能有1个点。假设正方形边长为$a$a，对角线长度为$r$r，那么有：$a^2+a^2=r^2$a2+a2=r2，那么$a=\frac{r}{\sqrt{2}}$a=2​r​。

代码

public static class Algorithm
{
    public static List<Vector2> Sample2D(float width, float height, float r, int k = 30)
    {
        return Sample2D((int)DateTime.Now.Ticks, width, height, r, k);
    }

    public static List<Vector2> Sample2D(int seed, float width, float height, float r, int k = 30)
    {
        // STEP 0

        // 维度，平面就是2维
        var n = 2;

        // 计算出合理的cell大小
        // cell是一个正方形，为了保证每个cell内部不可能出现多个点，那么cell内的任意点最远距离不能大于r
        // 因为cell内最长的距离是对角线，假设对角线长度是r，那边长就是下面的cell_size
        var cell_size = r / Math.Sqrt(n);

        // 计算出有多少行列的cell
        var cols = (int)Math.Ceiling(width / cell_size);
        var rows = (int)Math.Ceiling(height / cell_size);

        // cells记录了所有合法的点
        var cells = new List<Vector2>();

        // grids记录了每个cell内的点在cells里的索引，-1表示没有点
        var grids = new int[rows, cols];
        for (var i = 0; i < rows; ++i) {
            for (var j = 0; j < cols; ++j) {
                grids[i, j] = -1;
            }
        }

        // STEP 1
        var random = new Random(seed);

        // 随机选一个起始点
        var x0 = new Vector2(random.Range(width), random.Range(height));
        var col = (int)Math.Floor(x0.x / cell_size);
        var row = (int)Math.Floor(x0.y / cell_size);

        var x0_idx = cells.Count;
        cells.Add(x0);
        grids[row, col] = x0_idx;

        var active_list = new List<int>();
        active_list.Add(x0_idx);

        // STEP 2
        while (active_list.Count > 0) {
            // 随机选一个待处理的点xi
            var xi_idx = active_list[random.Range(active_list.Count)]; // 区间是[0,1)，不用担心溢出。
            var xi = cells[xi_idx];
            var found = false;

            // 以xi为中点，随机找与xi距离在[r,2r)的点xk，并判断该点的合法性
            // 重复k次，如果都找不到，则把xi从active_list中去掉，认为xi附近已经没有合法点了
            for (var i = 0; i < k; ++i) {
                var dir = random.insideUnitCircle();
                var xk = xi + (dir.normalized * r + dir * r); // [r,2r)
                if (xk.x < 0 || xk.x >= width || xk.y < 0 || xk.y >= height) {
                    continue;
                }

                col = (int)Math.Floor(xk.x / cell_size);
                row = (int)Math.Floor(xk.y / cell_size);

                if (grids[row, col] != -1) {
                    continue;
                }

                // 要判断xk的合法性，就是要判断有附近没有点与xk的距离小于r
                // 由于cell的边长小于r，所以只测试xk所在的cell的九宫格是不够的（考虑xk正好处于cell的边缘的情况）
                // 正确做法是以xk为中心，做一个边长为2r的正方形，测试这个正方形覆盖到所有cell
                var ok = true;
                var min_r = (int)Math.Floor((xk.y - r) / cell_size);
                var max_r = (int)Math.Floor((xk.y + r) / cell_size);
                var min_c = (int)Math.Floor((xk.x - r) / cell_size);
                var max_c = (int)Math.Floor((xk.x + r) / cell_size);
                for (var or = min_r; or <= max_r; ++or) {
                    if (or < 0 || or >= rows) {
                        continue;
                    }

                    for (var oc = min_c; oc <= max_c; ++oc) {
                        if (oc < 0 || oc >= cols) {
                            continue;
                        }

                        var xj_idx = grids[or, oc];
                        if (xj_idx != -1) {
                            var xj = cells[xj_idx];
                            var dist = (xj - xk).magnitude;
                            if (dist < r) {
                                ok = false;
                                goto end_of_distance_check;
                            }
                        }
                    }
                }

                end_of_distance_check:
                if (ok) {
                    var xk_idx = cells.Count;
                    cells.Add(xk);

                    grids[row, col] = xk_idx;
                    active_list.Add(xk_idx);

                    found = true;
                    break;
                }
            }

            if (!found) {
                active_list.Remove(xi_idx);
            }
        }

        return cells;
    }
}

// 测试代码
class Program
{
    static void Main(string[] args)
    {
        var width = 1024;
        var height = 1024;
        var r = 50f;
        var points = Algorithm.Sample2D(width, height, r);

        var image = new Bitmap(width, height);
        using (var graphics = Graphics.FromImage(image)) {
            graphics.FillRectangle(Brushes.Black, 0f, 0f, width, height);

            var dot_r = 3f;
            var pen = new Pen(Color.DarkRed, 2f);
            foreach (var p in points) {
                graphics.FillEllipse(Brushes.Yellow, p.x - dot_r, p.y - dot_r, 2f * dot_r, 2f * dot_r);
                graphics.DrawEllipse(pen, p.x - r / 2f, p.y - r / 2f, r, r);
            }
        }

        image.Save("out.png");
    }
}

输出的图片：可以看到黄点分布随机且比较平均，且任意2个黄点之间的距离都小于$r$r（红色圆的半径是$r/2$r/2）