效果图
应用效果图
算法
关键代码如下：
// 以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); } } }
红点是得到的采样点；
黄圈是采样点的圆域；
