泊松分布采样 （Poisson-Disk-Sample）代码及详细注释【OpenCV】

// 以center为圆心radius为半径的圆环范围内随机产生新的采样点
template<unsigned int N, class T>
void sample_annulus(T radius, const Vec<N, T> ¢re, unsigned int &seed, Vec<N, T> &x)
{
    Vec<N, T> r;
    for(;;)
    {
        for(unsigned int i = 0; i < N; ++i)
        {
      // r为 [-2, 2] 之间的随机数
            r[i] = 4 * (randhash(seed++) / (T)UINT_MAX - (T)0.5);
        }
    // 计算 r 的平方和
        T r2 = mag2(r);
    // 圆环范围在[radius+1, radius+2]
        if(r2 > 1 && r2 <= 4)
            break;
    }
  // 新采样点坐标
    x = centre + radius * r;
}
// 函数：将原先二维/三维的网格索引转换为一维
// 参数：dimension表示行列的网格数
// x 表示当前采样点所在的网格索引
template<unsigned int N, class T>
unsigned long int n_dimensional_array_index(const Vec<N, unsigned int> &dimensions, const Vec<N, T> &x)
{
  // k保存了全局网格的一维索引（列主序）
    unsigned long int k = 0;
    if(x[N - 1] >= 0)
    {
    //  保存网格行索引（二维的情况下）
        k = (unsigned long int)x[N - 1];
    // 如果超出网格索引范围则截断
        if(k >= dimensions[N - 1])
      k = dimensions[N - 1] - 1;
    }
    for(unsigned int i = N - 1; i > 0; --i)
    {
    // 得到网格的列数，并乘上之前保存的网格行索引（二维的情况下）
        k *= dimensions[i - 1];
        if(x[i - 1] >= 0)
        {
      // 保存网格的列索引（二维的情况下）
            unsigned int j = (int)x[i - 1];
      // 如果超出网格索引范围则截断
            if(j >= dimensions[i - 1]) 
        j = dimensions[i - 1] - 1;
      // 加上当前所在列的索引
            k += j;
        }
    }
    return k;
}
template<unsigned int N, class T>
void bluenoise_sample(T radius, Vec<N, T> xmin, Vec<N, T> xmax, std::vector<Vec<N, T> > &sample,
                      unsigned int seed = 0, int max_sample_attempts = 30)
{
  // 保存生成的有效采样点的坐标
    sample.clear();
  // 存储当前“激活”的采样点索引
    std::vector<unsigned int> active_list;
    // 用于加速比较的网格
  // 圆的半径略大于网格的对角线；一个网格存储一个采样点
    T grid_dx = T(0.999) * radius / std::sqrt((T)N); 
  // [0]存储了网格的列数，[1]存储了网格的行数
    Vec<N, unsigned int> dimensions;
  // 总的网格数
    unsigned long int total_array_size = 1;
    for(unsigned int i = 0; i < N; ++i)
    {
    // 计算可分配多少个网格
        dimensions[i] = (unsigned int)std::ceil((xmax[i] - xmin[i]) / grid_dx);
        total_array_size *= dimensions[i];
    }
    std::vector<int> accel(total_array_size, -1); // 开辟大小为总网格数的空间，用于存储其中采样点的索引；每个元素的初始值为-1
    // ！第一个采样点！
    Vec<N, T> x;
  // 在 xmin 与 xmax 范围内随机一个坐标
    for(unsigned int i = 0; i < N; ++i)
    {
        x[i] = (xmax[i] - xmin[i]) * (randhash(seed++) / (T)UINT_MAX) + xmin[i];
    }
    sample.push_back(x);
    active_list.push_back(0);
  // 将第一个采样点的网格索引（一维）
    unsigned int k = n_dimensional_array_index(dimensions, (x - xmin) / grid_dx);
  // 起始采样点的索引为0
    accel[k] = 0;
  // 当“激活”采样点列表不为空时
    while(!active_list.empty())
    {
    // 在 0 到 active_list.size() - 0.0001f 之间随机一个非负整数（从当前“激活”点列表中随机取一个采样点）
        unsigned int r = int(randhashf(seed++, 0, active_list.size() - 0.0001f));
    // 取得当前随机“激活”采样点的索引 p
        int p = active_list[r];
        bool found_sample = false;
        Vec<N, unsigned int> j, jmin, jmax;
    // 在当前“激活”点附近迭代搜索
        for(int attempt = 0; attempt < max_sample_attempts; ++attempt)
        {
      // 在“激活”点附近产生新的采样点坐标 x
            sample_annulus(radius, sample[p], seed, x);
            // 检查新采样点是否在边界范围内，否则拒绝该采样
            for(unsigned int i = 0; i < N; ++i)
            {
                if(x[i] < xmin[i] || x[i] > xmax[i])
                    goto reject_sample;
            }
      //  计算新采样点圆域内的网格左右、上下边界
            for(unsigned int i = 0; i < N; ++i)
            {
        // 检测新采样点的左侧（上侧）是否超过网格边界，如果超出，则截断
                int thismin = (int)((x[i] - radius - xmin[i]) / grid_dx);
                if(thismin < 0) 
          thismin = 0;
                else if(thismin >= (int)dimensions[i]) 
          thismin = dimensions[i] - 1;
                jmin[i] = (unsigned int)thismin;
        // 检测新采样点的右侧（下侧）是否超过网格边界，如果超出，则截断
                int thismax = (int)((x[i] + radius - xmin[i]) / grid_dx);
                if(thismax < 0) 
          thismax = 0;
                else if(thismax >= (int)dimensions[i]) 
          thismax = dimensions[i] - 1;
                jmax[i] = (unsigned int)thismax;
            }
      // 对新采样点附近的网格进行遍历（从最左上侧开始）
            for(j = jmin;;)
            {
        // 二维网格索引 j 转为一维网格索引 k
                k = n_dimensional_array_index(dimensions, j);
        // 如果该范围内存在一个不同于p（当前“激活”采样点的索引）的有效采样点
                if(accel[k] >= 0 && accel[k] != p) // if there is a sample point different from p
                {
          // 通过 一维网格索引 k 得到 采样点的索引 accel[k]
          // 根据采样点的索引 得到 采样点的坐标 ；
          // 如果该有效采样点距离 x 太近，则拒绝
                    if(dist2(x, sample[accel[k]]) < sqr(radius))
                        goto reject_sample;
                }
                for(unsigned int i = 0; i < N; ++i)
                {
          // 右（下）一个网格
                    ++j[i];
          // 如果其未超出边界范围
                    if(j[i] <= jmax[i])
                    {
            // 未超出则跳出本层循环， 继续上层循环的距离判断
            // 注意：下次再进入本for循环，i还是0，即++j[i]表示向右遍历
                        break;
                    }
          // 如果其超出边界范围
                    else
                    {
            // 范围内的最后一个网格，则结束范围内的网格迭代
                        if(i == N - 1) 
              goto done_j_loop;
                        else 
              // 恢复j[i]为起始位置（比如，j[0]的话就是最左侧）
              j[i] = jmin[i]; 
              // 注意：接着for循环++j[i]表示往下一个网格
                    }
                }
            }
done_j_loop:
            // 找到新的采样点
            found_sample = true;
            break;
    reject_sample:
      //如果x 离已存在的采样点太近
            ; 
        }
        if(found_sample)
        {
            size_t q = sample.size(); // 新采样点的索引
            sample.push_back(x);
            active_list.push_back(q);// 更新为新“激活”采样点
            k = n_dimensional_array_index(dimensions, (x - xmin) / grid_dx);
            accel[k] = (int)q;//  保存新“激活”采样点的索引
        }
        else
        {
      // 如果在p的圆环范围内没有找到新的采样点，则将当前激活列表中的第 r 项移除（即移除采样点索引p）
            erase_unordered(active_list, r);
        }
    }
}