Paper

Deterministic Blue Noise Samplingby Solving Largest Empty Circle Problems
Yoshihiro KANAMORI† , Zoltan SZEGO†† , Tomoyuki NISHITA†† (Member)† University of Tsukuba†† The University of Tokyo〈Summary〉 Sampling patterns with a blue noise distribution are widely used in many

areas of computer graphics, yet their eﬃcient generation remains a diﬃcult problem. We
propose a method to generate point sets with a blue noise distribution using a deterministic
algorithm with no preprocessing. We insert each new sample at the center of the largest
empty circle in the point set, which is obtained by calculating the Delaunay-triangulationof the set and ﬁnding the triangle with the largest circumcircle. Our method supports
adaptive sampling according to a user-speciﬁed density function, as well as specifying theexact number of required samples. It can also be extended to perform sampling on athree-dimensional curved surface.
Key words: Blue noise sampling, Delaunay triangulation, halftoning

than some of the alternatives, the proposed method

1.

Introduction

results in patterns with better spectral properties.
The proposed method supports adaptive sampling

Sampling is a fundamental part of various graph-with a simple modiﬁcation to the algorithm. Fur-ics applications such as image processing or render-

thermore, it can be naturally extended to handle

ing. With the number of samples allowed often being

arbitrary sampling domains, including volumes andlimited in these applications, getting the best-quality

three-dimensional curved surfaces.

results out of as few samples as possible is essential.

2.

For this purpose, one of the most often used sam-Related Workpling patterns are Poisson Disk patterns, deﬁned bya minimum distance between all neighboring sample
points. These have the spectral properties of blue
noise, deﬁned as having a spectrum with low energy
in the low frequencies. This property has been shown
to be very eﬀective at providing good sampling witha low amount of noise or aliasing, as evidenced by the

In this section, we ﬁrst introduce the basic framework for evaluating the quality of a point set basedon its spectral properties, and then describe the twocategories of methods for generating point sets withblue noise properties: methods based on dart throwing and tile-based methods.

fact that it also occurs naturally, for example, in the

2.1

placement of photoreceptor cells in the eye11).

Here we brieﬂy explains how the quality of a blue

Quality evaluation framework

This paper presents a novel method that produces

noise distributed point set is evaluated, as detailed

point sets deterministically, with high-quality bluein an in-depth survey5). We assume the samples arenoise distributions without any precomputation, us-generated on a two-dimensional domain [0, 1)2 , anding a geometrically-based algorithm to determine theconsider N point samples xi as a signal given by theexact position for each new sample. While slower

sum of Dirac-deltas (δ) located at the position of each1The Journal of the IIEEJ vol. 40 no. 210(2011)
point:

costly iterative process that constructs the Voronoi-Nδ(x − xi ).

(1)

i=1

diagram for the entire point set repeatedly.
In 2008, Wei10)proposed a method to generate

Taking the Fourier-transform of this signal gives theblue noise distributed samples at high speed on pro-

spectrum for a single point set, and averaging many

grammable GPUs based on the observation that ar-of these gives the characteristic power spectrum P (f )eas of the sampling domain that are suﬃciently far

for the method used, with f being the frequency. Theapart do not aﬀect each other, and can therefore be

typical blue noise spectrum has a low-energy low-

handled independently in parallel. The method pro-frequency band surrounding the center, followed byduces good quality results at high speed thanks toa sharp transition to higher energies, ﬁnally ﬂatten-

the high performance provided by GPUs; however, it

ing out for high frequencies. The spectrum is radiallyis diﬃcult to control the exact number of produced

symmetric, therefore it can be represented by a one-samples due to the random nature of the algorithm.

dimensional radially averaged power spectrum by av-

Our method can deterministically generate any given

eraging P (f ) in concentric rings, each corresponding

number of samples, which is useful for several applica-

to a frequency:

tions including object placement. Additionally, while
N (fr )

Pr (fr ) =1
P (C(i)),N (fr ) i=1

extending Wei’s method to non-Euclidean domains is(2)

where N (fr ) is the number of samples that are takenfrom the spectrum in each ring with radius fr . C(i)not trivial, our method can naturally handle samplingon spherical or arbitrary polygonal surfaces.
2.3

Tile-based methods

gives the position for the ith sample in the ring. The

Methods based on various tiling schemes use a dif-variance of the power for each frequency, s2 (fr ), is

ferent approach than dart throwing. Instead of rely-used to derive the anisotropy of the spectrum:ing on randomness, they generally precalculate dis-s2 (fr )Ar (fr ) = 2,
Pr (fr )tributions for a given number of predeﬁned tiles, and
(3)

where

then arrange those on the sampling domain, allowing
for better performance at the expense of quality and

N (fr )1
s2 (fr ) =(P (C(i)) − Pr (fr ))2 .
N (fr ) i=1(4)

precomputation time.
For example, Wang tiles are used by various meth-

High-quality blue noise spectra are characterized

ods. Hiller et al.3)introduced the use of Wang tiles

by a wide central low-energy ring and low anisotropy

for generating a point set. However, the power spec-(around -10dB).

trum of the method is of low quality. Kopf et al.4)

2.2

Methods based on dart throwing

used Wang tiles recursively to produce hierarchicalDart throwing, introduced by Cook2), generates

sampling patterns at high speed from a large precom-

Poisson disk distributions by ﬁrst generating uni-

puted dataset, but the spectrum still has unwanted

formly distributed points, then rejecting points thatspikes and noise.

are too close to each other given a minimum sepa-

Ostromoukhov9)proposed a diﬀerent tiling scheme

ration distance. While easy to implement, this algo-

based on hierarchically subdivided polyominoes. The

rithm is slow and diﬃcult to control because a min-method produces reasonably good results at highimum distance needs to be given instead of the re-speed, however, it still requires a rather complicatedquired number of samples. Dart throwing has beenand expensive precomputing step. The variety andoptimized and extended in a number of ways, one

randomness of the results are limited by the numberof which is called relaxation dart throwing7), whereof pregenerated tiling variations.

the minimum separation distance starts out large

Our method produces point sets with high-quality

and is gradually reduced, producing hierarchical sam-

spectra, and is not limited in the variety of outcomes

ples somewhat faster than traditional dart throwing.

by any kind of precomputed data. Note that our

A method called Lloyd’s relaxation6)is often applied

method does not ﬁt into the two categories presented

to these results as a post-process, which is a rather

in Sections 2.3 and 2.2, it should instead be classi-

2

Paper : Deterministic Blue Noise Samplingﬁed as a geometrically-based method that selects theexact position for each sample point.

3.

Proposed Algorithm

In this section, we present the algorithm used togenerate blue noise distributed point samples in detail, going into speciﬁcs for both the two-dimensional
and three-dimensional cases.
3.1

Fig. 2

Overview

Our method generates samples deterministically

Illustration of a completely triangulated
toroidal domain. A few triangles and theirneighbors are highlighted to emphasize thewrap-around eﬀect.

and sequentially. Each new sample is placed at thecenter of the largest empty circle in the already ex-gles with a circumcircle whose center lies outside of

isting set of samples. The intuitive reasoning for thethe domain [0, 1)2 are excluded from the search.

algorithm used is that given a set of samples, the next

3.2

one should go in the most sparsely sampled area soUsing the above algorithm as-is results in a point

far, and should be equidistant from the other samples

set whose Delaunay triangulation is constrained by

8)

Using a toroidal domain

surrounding it. Finding the largest empty circle isthe boundaries of the 2D domain [0, 1)2 . When aequivalent to ﬁnding the triangle with the largest cir-triangle’s circumcenter lies outside this domain, it iscumcircle in the Delaunay triangulation of the point

ignored and the next largest triangle is used instead.

set, a triangulation which guarantees by deﬁnition

As a consequence, elongated triangles tend to clus-that no other point lies within any triangle’s circum-

ter up near the edges, resulting in undesirable gaps

circle. The process is repeated, gradually ﬁlling up

around boundaries.

the sampling domain with points.

To alleviate this problem, the sampling domain can

Initially, the input point set consists of a few seed

be changed from a simple two-dimensional domain topoints, which serve as a means of providing determin-a toroidal domain. That is, the opposite boundaries

istic control over the pseudorandom outcome of theof the domain “wrap around”, the top boundary be-process. The procedure to generate new points is as

ing connected to the bottom boundary, and similarly

follows:

for the left and right boundaries, making the domain
topologically equivalent to the surface of a torus (Fig-

GeneratePoints(N , Points):// N : the total number of samples required// Points : the vector of sample points// (initialized with a small arbitrary set of seed points)
Delaunay := InitializeDelaunayTriangulation(Points);
while (Points.length < N )C := Delaunay.findLargestCircumcircle();Points.append(C.center);Delaunay.insertPointAndUpdateTriangulation(C.center)end

ure 2).
Once the domain wraps around, there are no more“boundaries” to speak of, and therefore the elongated triangles also disappear, resulting in much morebalanced shapes for all generated triangles.

Suchwell-shaped triangles allow us to omit very timeconsuming searches; in the case of an unbalanced triangle, the circumcenter tends to lie outside the trian-See Figure 1 for a graphical overview of one step in
the inner loop of the algorithm described above.

gle in question, and thus, when inserting a new point,
we must specify which triangles should be updated by

In order to ﬁnd the largest circumcircle quickly, insearching for the triangle that the new point belongs

addition to the Delaunay triangulation, we maintainto. In contrast, well-shaped triangles contains theira heap data structure containing all of the triangles,

own circumcenters in a large majority of cases (99%sorted by the radius of their circumcircle. The largestof all triangles in our experiments), and thus we cantriangle can therefore be obtained in constant time,

ﬁnd the triangle containing the point to be inserted

and updates to the triangulation (such as insertingin O(1). In case that a circumcenter lies slightly out-

and modifying triangles) happen in O(log N ). Trian-side the triangle in question, we clamp the point to3

The Journal of the IIEEJ vol. 40 no. 210(2011)(a)Fig. 1(c)(b)

Overview of a single step of the proposed algorithm. In order: (a) an
initial set of points, (b) their Delaunay triangulation with the largest
empty circle shown, (c) the new sample point immediately after being
added, and (d) the updated triangulation. Note that in this example,
only circles with centers within the area bounds were considered. Fora more correct approach, see Section 3.2.

Fig. 4
Fig. 3

Example of the appearance of the undesirable regular patterns mentioned in Section 3.3. The pattern in the generated points(left) causes hexagonally arranged spikes inthe spectrum (right).

be inserted within the triangle.
3.3

(d)

Avoiding regular patterns

While triangles with elongated edges cause loss ofsampling quality and speed, triangles that are too

Adjusting the position of a sample point inside an equilateral triangle.

triangle containing the point is close to being equilateral up to a certain threshold. If it is, the new point
is ﬁrst oﬀset by a small amount (δx, δy), such thatit stays within its triangle, in a pseudorandom direction depending only on the point’s original position,in this case, to keep our algorithm deterministic (seeFigure 4);
θ = (x + y) Cπ,δx = d cos(θ),δy = d sin(θ),

regular and close to equilateral are also undesirable.

where x and y are the original coordinates, C is an

During the execution of the proposed algorithm, it

arbitrary constant and d is the distance to move theproduces results that are often uneven, with notice-point. For the distance, we used 0.2 times the tri-able large areas of regular patterns that show up as

angle’s edge length (determined empirically). This

spikes on the power spectrum. See Figure 3 for an

modiﬁcation eliminates the formation of regular pat-example.

terns and ﬁxes the occasional spikes in the power

The reason these patterns emerge is because of thespectrum.

way the algorithm behaves for equilateral triangles. If

3.4

the Delaunay triangulation happens to include a clus-

With a simple modiﬁcation, the proposed method

ter of equilateral triangles sometime during the pro-

can also perform adaptive sampling based on a user-

cess, and one of them is selected as the one with thespeciﬁed importance function f : [0, 1)2 → [0, 1], thatlargest empty circle, the new point will form smaller

speciﬁes how dense the samples should be in an area.

equilateral triangles, with further subdivision of the

For all of the triangles in the heap data structurearea only preserving the pattern. The result is a regu-

mentioned in Section 3.1, the radius of the circumcir-

lar triangular grid, obviously unsuitable for blue noisecle that acts as the sort key is weighted by a sample

sampling. To prevent this, an additional step is in-from the importance function, taken at the center ofcluded in the algorithm: when inserting a new point

the circumcircle. This way, triangles that lie in areas

into the Delaunay triangulation, we ﬁrst check if therequiring more thorough sampling are given more pri-

4

Adaptive samplingPaper : Deterministic Blue Noise Samplingority, resulting in sample points that are distributed

Table 1

according to the importance function.
3.5

Sampling on a 3D surface

Time to generate a given number of samples (msec)
Number of samples

Method20,00050,000100,000

Recursive Wang tiles

0.59

1.35

2.66

Bridson et al.’s dart throwing

165

420

829

cess is mostly the same: starting from a sparse set of

Proposed method (1 thread)

132

378

853

points and their points are added in triangles with the

Proposed method (4 threads)

92

241

502

Our method can be extended to generate samples on three-dimensional surfaces. The basic pro-largest circumcircle. The new points are projected
onto the surface and the triangulation is updated.
In case the surface is speciﬁed by a complex polyg-

Table 2

Scaling by parallelization for 100,000 samples

onal mesh, the initial set of points and their con-

Subdivision

Time

Scaling

nectivity can be obtained with a mesh simpliﬁcationnone

853 ms

1.00×

2×2

559 ms1.52×

4×4

502 ms

1.69×

8×8

502 ms

1.69×

method such as edge decimation. Adaptive samplingis also possible if an importance function is deﬁned
on the surface.
3.6

Improving performance by parallelization

4.

Results

In order to further improve the performance of theWe implemented a prototype system that uses our

proposed method, we can parallelize the process us-algorithm in C++. The environment used for our

ing modern multicore CPUs. One observation whenexperiments was a standard desktop PC with an Intel

testing our algorithm was that the updating of theCore 2 Quad Q6700 CPU at 2.66GHz and 2GB of

Delaunay property generally does not propagate very

RAM. We used the Intel Compiler version 10.1 and

far in the triangulation. In fact, the number of trian-

the OpenMP extension for parallelization.

gle ﬂips after a point insertion was always less than

We compared the performance and quality of the

10. In other words, when the number of points is suf-

proposed method to two other related techniques,ﬁciently large, individual steps of the algorithm op-

Recursive Wang tiles4)and an optimized variation oferate fairly locally. A similar observation is also ex-

the dart throwing algorithm utilized by Bridson et10)

in his highly parallelized sample gen-al.1)We performed experiments for a varying number

eration technique. We can also observe that the per-of samples, namely, uniform sampling with 20,000,

formance of our algorithm decreases as the number of

50,000 and 100,000 samples generated. Table 1 com-

generated points get bigger, due to the O(N log N )pares the speed measurements for each technique.

runtime of repeatedly updating a heap for N itera-

Figures 5 and 6 show the spectra for each method, astions.

well as the radially averaged energy and anisotropy

ploited by Wei

We therefore employ a divide-and-conquer strategy,

graphs, for 20,000 and 50,000 samples, respectively.

by cutting up the original [0, 1)2 domain into subre-

Performance-wise, Recursive Wang tiles performed

gions. Theoretically, if we divide the region into K

several orders of magnitude faster than the alterna-

subregions, the cost becomes O(K ×N/K log N/K) =tives, however, it did so at the expense of quality. The

O(N log N/K), which is smaller than O(N log N ).

larger the number of samples required, the noisier theNote that no guarantees are made about the Poisson

output becomes, as is evident from the spike patterns

Disk property of samples at the edges of neighboringand noise in the spectrum and anisotropy graphs.

subregions this way, since the subtasks do not take

Our method, while having a speed faster than theconstraints due to other regions into account. How-

optimized dart throwing variation, produced slightly

ever, we found that in the actual results the borders

better quality results, based on the wider inner ringare barely noticeable, and the resulting spectra are in-of the power spectrum.

distinguishable from the single-threaded version. For
details on the scaling achievable, see Section 4.

Table 2 shows the scaling achieved by parallelization. The scaling does not increase as the theory due
5The Journal of the IIEEJ vol. 40 no. 210(2011)㪉㪇

frequency
anisotropy [dB]anisotropy [dB]

㪇

㪇

㪈㪇

㪈

㪇

(b) Dart throwing(c) Proposed method㪉㪇anisotropy [dB]frequency

㪇

㪈㪇

㪈

㪇

frequency

㪇

㪈㪇

㪈㪇

frequency(a) Recursive Wang tiles㪈㪇

frequency
anisotropy [dB]

㪈

㪉

power

power

power
anisotropy [dB]frequency

Spectral quality comparison chart for 20,000 generated samples.

㪉Fig. 6㪇frequency(a) Recursive Wang tiles

㪇frequency㪈㪇㪈㪇frequency

Fig. 5㪈㪇frequency
anisotropy [dB]

㪈

㪉

power

power

power

㪉

frequency

frequency(b) Dart throwing(c) Proposed method

Spectral quality comparison chart for 50,000 generated samples.

to the overheads of the process. Because the subdivision slightly causes quality loss at boundaries, weconclude 4 × 4 subdivision is the optimal.

halftoning of a grayscale input image.
Figure 8 shows the results of producing samples
using our method on spherical HDR light maps. The

Figure 7 shows the examples of adaptive sam-importance function was given by converting the color

pling on a 2D domain, demonstrated by point-basedof each pixel into a grayscale value and biasing the6

Paper : Deterministic Blue Noise Sampling

Fig. 7

(a)

(b)

(c)

(d)

Results of halftoning the image of a mask.
(a) The grayscale original. (b)–(d) Halftoned
results created with the proposed method
consisting of 20,000, 50,000 and 100,000
points, respectively.

Fig. 8generated samples toward bright regions only.
Finally, Figure 9 demonstrates that our method canbe extended to three dimension via Delaunay tetrahe-

Results of adaptively sampling a variety ofspherical HDR radiance environment maps
using our method. In each case, 3,000 samples are placed on the unit sphere, displayed
as dots.

dralization. The three-dimensional spectrum shows
the samples exhibit blue noise property.

Vision Modeling and Visualization Conference 2001, pages
265–272, 2001.

5.

Conclusion and Future WorkWe described a method to generate sampling patterns with a high quality blue noise distribution. Our
method is deterministic, produces the exact numberof required samples, requires no precomputation, andhandles adaptive sampling based on a user-speciﬁed4) Johannes Kopf, Daniel Cohen-Or, Oliver Deussen, andDani Lischinski. Recursive wang tiles for real-time blue
noise. In SIGGRAPH ’06: ACM SIGGRAPH 2006 Papers,
pages 509–518, New York, NY, USA, 2006.
5) Ares Lagae and Philip Dutré. A comparison of methods
for generating Poisson disk distributions.

Report CW

459, Department of Computer Science, K.U.Leuven, Leuven, Belgium, August 2006.
6) S. Lloyd. Least squares quantization in PCM. Information Theory, IEEE Transactions on, 28(2):129–137, Mar

importance function. Furthermore, our method canbe easily extended higher dimensions than two thanks
to the theoretical simplicity.

For future work, wewould like to gain further speed up of our method

1982.
7) Michael McCool and Eugene Fiume. Hierarchical poissondisk sampling distributions. In Proceedings of the conference on Graphics interface ’92, pages 94–105, San Francisco, CA, USA, 1992.
8) Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, andusing CUDA.

Sung Nok Chiu. Spatial Tessellations: Concepts and Ap-

References
1) Robert Bridson, Jim Houriham, and Marcus Nordenstam. Curl-noise for procedural fluid flow. ACM Trans.
Graph., 26(3):46, 2007.
2) Robert L. Cook. Stochastic sampling in computer graphics. ACM Transactions on Graphics, 5(1):51–72, 1986.
3) Stefan Hiller, Oliver Deussen, and Alexander Keller.
Tiled blue noise samples. In VMV ’01: Proceedings of theplications of Voronoi Diagrams. Wiley, 1992.
9) Victor Ostromoukhov. Sampling with polyominoes. ACM
Transactions on Graphics, 26(3):78, 2007.
10) Li-Yi Wei. Parallel poisson disk sampling. ACM Transactions on Graphics, 27(3):1–9, 2008.
11) J. I. Yellott Jr. Spectral consequences of photoreceptor
sampling in the rhesus retina. Science, 221(4608):382–
385, 1983.

(Submitted July 16, 2010)
7

The Journal of the IIEEJ vol. 40 no. 210(2011)

Fig. 9The spectrum of 10,000 3D samples in the unit cube. The threedimensional spectrum is shown as a series of 2D slices. The low-energy
spherical region in the middle indicates a blue noise distribution.

(Received November 5, 2010)8