任意维度的快速泊松碟采样

发表于 | 分类于
| 字数统计 2218

原文地址: https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf

原文提及的附加材料及源代码: https://dl.acm.org/doi/10.1145/1278780.1278807

罗伯特·布里德森 (Robert Bridson) 英属哥伦比亚大学 (University of British Columbia)

译者注:
前置知识: 蓝噪声, 泊松碟

摘要

在许多图形应用程序中，尤其是在渲染中，从蓝噪声 (blue noise) 分布中生成样本非常重要。但是，现有的高效技术并不能轻松生成二维以外的范围。这里我对飞镖投掷方法进行了简单的改造，可以在 O(N) 时长内在任意维度轻松生成泊松碟样本。

CR类别： I.3.0 [计算机图形]：综合

关键词：采样(sampling)，蓝噪声(blue noise)，泊松碟(Poisson disk)

1. 前言

蓝噪声采样模式 —— 以泊松碟分布生成为例，在这种模式中，用户提供密度参数 r ，所有样本至少相距 r —— 被许多应用程序视为理想的渲染（见库克的里程碑式论文【1986】）。很遗憾，基于舍选的朴素的泊松碟样本生成方法飞镖投掷法，效率很低。

不少论文已经在二维上克服了这一低效的缺陷；我着重研究了邓巴和汉弗莱【2006】的方法。他们的算法保持了“扇形”数据结构，这种结构可以高效地编码距离现有样本 r 至 2r 之间的平面区域的几何结构，并从该区域中高效均匀地采样。由于该区域中的每一个点是一个容许样本，而且包含现有样本的每个最大泊松碟采样必须也包含该区域内的一个点，他们的算法可以非常高效地生成所需分布。

但许多采样应用程序采用了三维或更高维度：用动态模糊（motion blur）或景深效果（epth-of-field）、许多专用动画系统等进行渲染。现有的二维快速蓝噪声采样器不能轻松推及到更高维度，因此从业者们想用均匀随机分布（尽管存在不如人意的聚集现象）、抖动/分层采样（可以减少但不能消除聚集）或者更结构化的分布诱导各向异性。

在这份初稿中，我呈现了一种新算法，很容易在任何维度实施，可保证在 O(N) 时长内生成 N 个泊松碟样本。与邓巴和汉弗莱的方法相似，候选样本只从现有样本附近的区域选择，但不用准确计算允许区域，而是采用舍选采样发现样本。

2. 算法

该算法将样本域的范围 Rn、样本之间的最小距离 r 以及常数 k 输入作为算法舍选之前选择样本的极限（一般 k = 30）。

图1：从算法得出的二维样本模式及相关多行平均周期图谱

step 0. 初始化一个 n 维背景网格，用于储存样本，加速空间搜索。我们选择单元格尺寸受限于 r / \(\sqrt{n}\)，这样每个单元格都将最多包含一个样本，网格可以用一个简单的 n 维数组来实现：默认 -1 表示没有样本，一个非负整数表示一个单元格内样本的索引。

step 1. 从样本域里均匀随机选择初始样本 X0。将其插入背景网格，用该索引（0）初始化 “活动列表”（样本索引的数组）。

step 2. 当 活动列表 不为空时，从中选择一个随机索引（称为 i）。从距离 Xi 半径 r 到 2r 之间的球形环中均匀地生成 k 个点。依次检查每个点是否位于现有样本的 r 距离内（背景网格只用于检测附近样本）。如果某个点距离现有样本足够远，将其作为下一个样本，并添加至活动列表中。如果尝试 k 次之后未发现这种点，则从活动列表中移除 i。

3. 分析

将第2步执行 2N - 1 次，产生 N 个样本：每执行一次要么产生一个新样本，将其加入活动列表，要么从活动列表中移除一个现有样本。每次执行第2步需要 O(k) 时长，由于 k 保持不变（一般比较小），所以算法是线性算法。

图1显示了二维结果；随附材料显示了三维结果。随附原型代码介绍了如何把空间维度作为一个参数的简单操作。

致谢

本研究获得加拿大自然科学与工程研究理事会资金的支持。

参考文献

库克R.L.,1986年，“计算机图形随机采样”，《ACM Trans.Graph》，第5,1篇。

邓巴D.和汉弗莱G.，2006年，“快速泊松碟样本生成的空间数据结构”，《ACM Trans.Graph.》，25，3，503-508。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#ifndef BLUENOISE_H 
#define BLUENOISE_H 

#include <cmath> 
#include <cstdlib> 

// Provides type Vec<N,T>, a C++ wrapper around T[N] to 
// encapsulate N-dimensional vectors. 
#include "vec.h"

// Produce a sample uniformly chosen from the annulus between 
// (radius,2*radius] around center x. Uses random() for underlying PRNG.
template<unsigned int N, class T> 
void sample_annulus(T radius, const Vec<N,T> &centre, Vec<N,T> &x)
{
    Vec<N,T> r;
    for(;;){ // simple rejection sampling approach
        for(unsigned int i=0; i<N; ++i){
            r[i]=4*(random()/(T)2147483647-(T)0.5);
        }
        T r2=mag2(r); // magnitude squared of r
        if(r2>1 && r2<=4)
            break;
    }
    x=centre+radius*r;
}

// Translate a multidimensional array index into a one dimensional array index.
// (rounds towards zero if T is float or double)
template<unsigned int N, class T> 
unsigned long int n_dimensional_array_index(const Vec<N,unsigned int> &dimensions,
                                            const Vec<N,T> &x)
{
    unsigned long int k=0;
    for(unsigned int i=N; 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;
}

// Generate blue noise samples at least radius apart in N dimensions 
// within box bounded by xmin and xmax. T should either be float or double. 
template<unsigned int N, class T>
void blue_noise_sample(T radius, Vec<N,T> xmin, Vec<N,T> xmax,
                       std::vector<Vec<N,T> > &sample, // contains samples on return
                       int max_sample_attempts=30)
{
    // initialize acceleration grid (step 0)
    T grid_dx=T(0.99999)*radius/std::sqrt((T)N);
    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 indicates no sample there;
                                                  // otherwise index of sample point
    // first sample (step 1)
    Vec<N,T> x;
    for(unsigned int i=0; i<N; ++i){
        x[i]=(xmax[i]-xmin[i])*(random()/(T)2147483647) + xmin[i];
    }
    sample.clear();
    sample.push_back(x);
    std::vector<unsigned int> active_list;
    active_list.push_back(0);
    unsigned int k=n_dimensional_array_index(dimensions, (x-xmin)/grid_dx);
    accel[k]=0;
    
    // generate the remaining samples (step 2)
    while(!active_list.empty()){
        unsigned int r=(unsigned int)((random()/(T)2147483648u)*active_list.size());
        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){
            sample_annulus(radius, sample[p], x);
             // check this sample is within bounds
             for(unsigned int i=0; i<N; ++i){
                 if(x[i]<xmin[i] || x[i]>xmax[i])
                     goto reject_sample;
             }
             // find range in acceleration grid that may contain interfering samples
             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;
            }
            // loop over the selected grid cells (a little obfuscated since this is,
            // in effect, a nested N-dimensional loop)
            for(j=jmin;;){
                // check if there's a sample at j that's too close to x
                k=n_dimensional_array_index(dimensions, j);
                if(accel[k]>=0 && accel[k]!=p){ // if there is a sample different from p
                    if(dist2(x, sample[accel[k]]) < radius*radius)
                        goto reject_sample; // proposed sample is too close to accel[k]
                }
                // move on to next j (N-dimensional nested version of ++j)
                for(unsigned int i=0; i<N; ++i){
                    ++j[i];
                    if(j[i]<=jmax[i]){
                        break;
                    }else{
                        if(i==N-1) goto done_j_loop;
                        else j[i]=jmin[i]; // and try incrementing the next dimension along
                    }
                }
            }
            done_j_loop:
            // if we made it here, we're good!
            found_sample=true;
            break;
            // if we goto here, x is too close to an existing sample
            reject_sample:
            ; // nothing to do except go to the next iteration in this loop
        }
        if(found_sample){
            unsigned int q=sample.size(); // the index of the new sample
            sample.push_back(x);
            active_list.push_back(q);
            k=n_dimensional_array_index(dimensions, (x-xmin)/grid_dx);
            accel[k]=(int)q;
        }else{
            // since couldn't find a sample on p's disk, remove p from the active list
            active_list[r]=active_list.back(); // overwrite with the last entry here
            active_list.pop_back(); // and delete the last entry
        }
    }
}
#endif

译者注:
javascript 实现: https://github.com/kchapelier/poisson-disk-sampling