ii“jgt” — 2011/12/20 — 15:51 — page 1 — #1

iiVol. [VOL], No. [ISS]: 1–??

Golden Ratio Sequences
For Low-Discrepancy Sampling
Colas Schretter and Leif Kobbelt
RWTH Aachen University, Germany

Paul-Olivier Dehaye ETH Zurich, Switzerland

Abstract. Most classical constructions of low-discrepancy point sets are based on
generalizations of the one-dimensional binary van der Corput sequence whose implementation requires non-trivial bit-operations. As an alternative we introduce thequasi-regular golden ratio sequences which are based on the fractional part of successive integer multiples of the golden ratio. By leveraging results from number theory
we show that point sets which evenly cover the unit square or disc can be computedby a simple incremental permutation of a generator golden ratio sequence. We compare ambient occlusion images generated with a Monte Carlo ray-tracer based on
random, Hammersley, blue noise and golden ratio point sets. The source code of
the ray-tracer used for our experiments is available online.

1.

Low-discrepancy point sets

Hammersley A classical approach for generating quasi-random sequencesof point samples in the unit square [0, 1)2 relies on generalizations of the vander Corput low-discrepancy sequence based on small prime numbers [Niederreiter 92]. For instance, the Hammersley point set of size N is{(H2 (i), i/N )}N
i=1(1)© A K Peters, Ltd.

11086-7651/06 $0.50 per page

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 2 — #2

i

2

i

journal of graphics tools

where Hb (i) enumerates the elements of the van der Corput sequence in baseb. The value Hb (i) is computed by reversing the digits of the integer representation of i in base b and taking the result as the fractional part of a fixed
point number. This operation can be implemented efficiently with simple bitwise operations for the special cases of bases equal to a power of two, whereb = 2p , ∀p > 1.
Halton When the number of samples is not known beforehand, Halton proposed a low-discrepancy sequence of self-avoiding points using different small
prime numbers as the bases of van der Corput sequences for successive dimensions. For instance, the points(H2 (i), H3 (i)),

∀i≥1

(2)

are elements of the Halton sequence in two dimensions. The rationale for thismethod is that indices i are expressed in base of prime numbers; therefore,
the digits of representations in base 2 and 3 are less likely to be correlated.
However, the discrepancy of the points {Hb (i)}N
i=1 is far from optimal when
N is not close to a power of b; hence, N needs to be close to a power of each of
the different primes, which is a prohibitive restriction even in two dimensions.
Blue noise Intuitively, clustering of points should be avoided and therefore
distances to the nearest neighbor point should be constant. This property isensured only by regular grids such as an hexagonal packing of disks. Since
regular sampling causes strong aliasing artifacts, it is often wanted that the
point set exhibits a blue noise frequency spectrum. A blue noise pattern isa random point distribution with attenuated response in the low-frequency
band. Several geometrical constructions such as Penrose tiling followed by
iterative Lloyd’s relaxation [Balzer et al. 09] or dart-throwing [Dunbar andHumphreys 06] have been applied for creating point sets with blue noiseproperty. Unfortunately, those heuristic methods are rather slow and often
complicated to implement.
Rank-1 lattices Alternatively, one can construct a rank-1 Fibonacci latticefor N equal to the Fibonacci number Fk , ∀k ≥ 3 as follows,{(hi · Fk−1 /Fk i, i/Fk )}N
i=1

(3)

where hxi = x − bxc. The orthogonal projection of points on the second axis
yields regular equispaced samples while the projection on the first axis shows
the quasi-regular distribution generated by fractional parts of multiples of
the ratios of successive Fibonacci numbers. This construction generalizes tohigher dimension [Sloan and Reztsov 02]; however, this lattice rule has provenoptimal discrepancy in two dimensions only [Niederreiter 92].

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 3 — #3

i

Schretter et al.: Golden Ratio Sequences For Low-Discrepancy Samplingi3

Unfortunately, equispacing of the coordinates in the second dimension isinherent to lattice rules based on generator vectors [Dammertz et al. 08] aswell as Hammersley sets. Alignments and systematic regularities will generatealiasing artifacts in image rendering applications. In his Instant Radiosity
paper, Keller proposed to jitter the Hammersley points as a workaround to
avoid aliasing [Keller 97]. The jitter is a small random displacement appliedin both dimensions since the van der Corput sequence in base 2 also generates
grid-aligned coordinates.

2.

Golden ratio sequences

This paper introduces the one-dimensional golden ratio sequences as an alternative to the dyadic scheme of van der Corput. We define the golden ratio
sequences by the fractional parts of a seed constant s ∈ [0, 1) plus an integer
multiple of the golden ratio φ, i.e.,
Gs (i) = hs + i · φi,

∀i≥1

(4)

with

√
1+ 5
≈ 1.618034 . . .
(5)2Note that the conjugate golden section τ = φ − 1 can be used in place of
φ since only fractional parts are retained. Hence, computing Gs (i + 1) given
Gs (i) only involves an addition and a test for integer overflow.
Golden ratio sequences are much faster and simpler to compute than van der
Corput sequences, even in base 2. Furthermore, computations can be carriedout most easily in the integer domain. For instance, if unsigned integers are
stored in 32 bits registers, a sequence of quasi-uniform numbers in [0, 232 )is produced by successive additions of bτ · 232 c = 2654435769. The registercan be initialized with any seed value and implicit overflows of the carry bit
alleviate the need for additional comparisons.
Figure 1 compares the five first values of the van der Corput in base 2 andthe golden ratio sequence with s = 0. In fact, the van der Corput sequencealways splits the largest interval between points in two and the ratio of thelargest interval to the smallest is two. In particular, when the number of
elements is a power of the base prime, all large intervals have been broken into
two small ones and the ratio is one, yielding a perfectly regular distribution.
Samples from the golden ratio sequence cover more evenly the unit interval for
any number of elements. In comparison to the van der Corput sequence, thedistance between the two extremes points (5 and 3) is larger. Simultaneously,
the largest gap between neighbor points (2 and 4) is smaller, resulting in
better avoidance between samples.
φ=

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 4 — #4

i

4

i

journal of graphics tools

4

5

2

1
5
3
van der Corput sequence2

41
Golden ratio sequence3

Figure 1. The five first elements of the van der Corput sequence in base 2 (top)and the golden ratio sequence with s = 0 (bottom).

The golden ratio conjugate τ is the most irrational1 from all irrational
numbers smaller than one [Coxeter 53]. As a consequence, the variability of
pairwise distances between points is minimized when using this constant infractional sequences. To our knowledge, the golden ratio has not yet beenused as a mean for the generation of low-discrepancy sequences in computer
graphics. The next section further exploits known results about this sequenceand describes a very fast algorithm for generating quasi-regular point sets.

3.

Golden point sets

A very fast method using permutations of an initial golden ratio sequence
allows the construction of low-discrepancy quasi-lattice point sets. The firstcoordinates of points in a golden point set are given by a golden ratio sequence,
while the second coordinates are a permutation of that same sequence. The
result is a irregular distribution whose orthogonal projections yields a golden
ratio sequence, as shown in figure 2. This approach is inspired by the Faure
point sets that are constructed from permutations of a set of coordinates thatare initialized with van der Corput in base 2 [Faure 92]. Golden point sets
are thus defined as{(Gs (i), Gs (σs (i)))}N(6)i=1
where the second coordinates are obtained by the permutation σs of [1, . . . , N ]that sorts the elements of {Gs (i)}N
i=1 .
While this description involves sorting a set of values, there is an alternativealgorithm that generates the same permutation without explicit sorting. It
relies on the fact that φ can be arbitrarily well approximated by ratios of consecutive Fibonacci numbers. Therefore, as long as the denominator is greater
1 Number theorists have related the quality of approximations of a real number by rational numbers to the size of the coefficients that appear in its continued fraction expansion.
Smaller coefficients in this expansion lead to reals with worse rational approximations. The
golden section φ, with a continued fraction expansion with constant coefficients 1 (minimal)thus has the worst approximations.

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 5 — #5

i

Schretter et al.: Golden Ratio Sequences For Low-Discrepancy Sampling

Lattice (regular)i5Golden set (irregular)

Figure 2. 16 points placed by a lattice rule and the similar but irregular golden
point set with s = 0 and N = 16.

or equal to N , comparison of the rational approximations to the fractionalparts is enough to compare the fractional parts themselves. Thereby, we canrecover σs with a simple incremental rule as follows.
Let the increment a = F2i+1 and the decrement b = F2i be the two smallest
Fibonacci numbers with successive even and odd index such that a + b ≥ N .
An algorithm for mapping the elements of a golden ratio sequence towards
the second dimension starts with the index σs (1) initialized to the position of
the minimum fractional part in the first dimension. This initial value is foundby retaining the minimum when generating the first dimension.
The second coordinate of the first point is assigned with Gs (σs (1)). The
remaining assignations are determined using a very simple iteration defined
as follows:(σs (i) + a if σs (i) ≤ b
σs (i + 1) =(7)
σs (i) − b otherwise
Note that the result can overflow after an increment by a and the assignation
should simply be skipped in this case. A facultative third pass could generate
another golden sequence for the first dimension, using a second seed value.
This variant produces different sets of coordinates for each dimensions.
For illustration, we give an example for N = 6 and s = 0. Then, a = F5 = 5,
b = F4 = 3 and σ0 (1) = 5. The sequence produced by the above recurrence
relation is [5, 2, 7, 4, 1, 6, 3]. The boldfaced number above N should be skipped
and we get the permutation σ0 = [5, 2, 4, 1, 6, 3]. In our implementation, thepermutation vector is not explicitly stored but assignments of coordinates forthe second dimension are directly executed. The execution time is linear with
N , as only one pass per dimension is required with a test and one addition or
subtraction per element.

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 6 — #6

i

6

i

journal of graphics tools

RandomHammersley

Blue noise

Golden set

Figure 3. Comparison of 100 irregular (see Figure 2) golden points with various
schemes and plots in both Cartesian (top) and polar (bottom) coordinates.

Density images

Hammersley

Golden set

Figure 4. The Hammersley point set and a golden point set are warped according toa clamped linear ramp (top) and a mixture of two multivariate Gaussians (bottom).

4.

Experiments

Figure 3 plots point sets in Cartesian and polar coordinate systems. One
dimension is interpreted as a normalized radius and the second is multipliedby 2π to yield an angle. The random sequence fails to fill gaps between
clusters. The other point sets exhibit the low-discrepancy property by more
uniformly covering the unit square. Note that the blue noise pattern contains
spurious alignments that are visible by strings of points in polar coordinates.
In a phyllotaxis interpretation [Coxeter 53], golden point sets in cylindricaland polar coordinates are realistic models for distributing scales on the skinof pineapples and seeds in the head of sunflowers.

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 7 — #7

i

Schretter et al.: Golden Ratio Sequences For Low-Discrepancy Sampling

RandomHammersley

Blue noisei7Golden set

Figure 5. Ambient occlusion integration with 64 (top) and 16 (middle) samples
per pixel. Golden point sets yield banding-free images with low-variance.

Golden point sets inherit from the many theoretical properties of the golden
ratio number. In particular, such sets are quasi-regular and thus provide
natural stratification while avoiding aliasing that would arise from regularalignements. Golden point sets have been compared with the Hammersley
set for adaptive sampling driven by a discrete probability density functions.
Figure 4 shows warped points. A sub-pixel accurate implementation of theinversion method with guide tables for hashing [Devroye 86, Page 96] was used
for warping points in amortized constant time.
Figure 5 compares ambient occlusion images computed by Monte Carlointegration. Aliasing (banding) artifacts impair the image computed with theHammersley point set and Cranley-Patterson rotations [Kollig and Keller 02].
The blue noise point set was generated with a variant of the iterative Lloyd’s
algorithm [Balzer et al. 09]. The experiment demonstrates that while the blue
noise model is more appropriate than Hammersley for antialiasing purpose, it
performs slightly worse than our quasi-lattice rule for numerical integration.
One reason is that even if blue noise point sets covers evenly the unit square,
they can loosen their stratification after warping (see Fig. 3, 2nd row, 3rd col.).

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 8 — #8

i

8

i

journal of graphics tools

References[Balzer et al. 09] M. Balzer, T. Schlömer, and O. Deussen. “CapacityConstrained Point Distributions: A Variant of Lloyd’s Method.” In Proc.
ACM SIGGRAPH 2009, pp. 86–93, 2009.
[Coxeter 53] H. S. M. Coxeter. “The Golden Section, Phyllotaxis, andWythoff’s Game.” Scripta Mathematica 19 (1953), 135–143.
[Dammertz et al. 08] S. Dammertz, H. Dammertz, and A. Keller. “Efficient
Search for Low-Dimensional Rank-1 Lattices with Applications in Graphics.” In Proc. Monte Carlo and Quasi-Monte Carlo Methods, pp. 271–287,
2008.
[Devroye 86] L. Devroye. Non-Uniform Random Variate Generation. SpringerVerlag, New York, 1986.
[Dunbar and Humphreys 06] D. Dunbar and G. Humphreys. “A spatial data
structure for fast Poisson-disk sample generation.” In Proc. ACM SIGGRAPH 2006, pp. 503–508, 2006.
[Faure 92] H. Faure. “Good permutations for extreme discrepancy.” Journalof Number Theory 42 (1992), 47–56.
[Keller 97] A. Keller. “Instant Radiosity.” In Proc. ACM SIGGRAPH ’97,
pp. 49–56, 1997.
[Kollig and Keller 02] T. Kollig and A. Keller. “Efficient Multidimensional
Sampling.” Computer Graphics Forum 21:3 (2002), 557–563.
[Niederreiter 92] H. Niederreiter. Random number generation and quasi-Monte
Carlo methods, CBMS-NSF Regional Conference Series in Applied Mathematics, 63. SIAM, 1992.
[Sloan and Reztsov 02] I. H. Sloan and A. V. Reztsov. “Component-bycomponent construction of good lattice rules.” Mathematics of Computation 71:237 (2002), 263–273.

Web Information:
C++ ray-tracer source code: http://jgt.akpeters.com/papers/Schretter2012/

Received [DATE]; accepted [DATE].

i

i
i

i

i

i
“jgt” — 2011/12/20 — 15:51 — page 9 — #9i

Schretter et al.: Golden Ratio Sequences For Low-Discrepancy Sampling

A.

i

9

Generating golden point sets in C++

1 void golden set(double points[][2], const unsigned int N) {2// set the initial first coordinate
3
double x = drand48();
4
double min = x;5
unsigned int idx = 0;789// set the first coordinatesfor(unsigned int i = 0; i < N; ++i) {points[i][1] = x;

11
12
13

// keep the minimum
if(x < min)min = x, idx = i;

15
16
17
18

// increment the coordinate
x += 0.618033988749894;
if(x >= 1) −−x;
}

20
21
22
23
24
25

// find the first Fibonacci >= N
unsigned int f = 1, fp = 1, parity = 0;
while(f + fp < N) {unsigned int tmp = f; f += fp, fp = tmp;++parity;
}

27
28
29
30

// set the increment and decrement
unsigned int inc = fp, dec = f;if(parity & 1)
inc = f, dec = fp;

32
33
34
35
36
37
38
39
40
41

// permute the first coordinates
points[0][0] = points[idx][1];for(unsigned int i = 1; i < N; ++i) {if(idx < dec) {idx += inc;if(idx >= N) idx −= dec;
} elseidx −= dec;
points[i][0] = points[idx][1];
}

43
44

// set the initial second coordinate
double y = drand48();

46
47
48

// set the second coordinatesfor(unsigned int i = 0; i < N; ++i) {points[i][1] = y;

50
51
52
53
54 }

// increment the coordinate
y += 0.618033988749894;
if(y >= 1) −−y;
}

i

i
i

i