Wavelet Noise
Robert L. Cook Tony DeRose
Pixar Animation Studios

Abstract
Noise functions are an essential building block for writing procedural shaders in 3D computer graphics. The original noise function
introduced by Ken Perlin is still the most popular because it is simple and fast, and many spectacular images have been made with it.
Nevertheless, it is prone to problems with aliasing and detail loss.
In this paper we analyze these problems and show that they are particularly severe when 3D noise is used to texture a 2D surface. We
use the theory of wavelets to create a new class of simple and fast
noise functions that avoid these problems.
CR Categories: I.3.3 [Picture/Image generation]: Antialiasing—[I.3.7]: Three-Dimensional Graphics and Realism—Color, shading, shadowing, and texture
Keywords: Multiresolution analysis, noise, procedural textures,
rendering, shading, texture synthesis, texturing, wavelets.

1

Introduction

Ever since the introduction of Perlin noise in [Perlin 1985], noisefunctions have been important for creating textures in 3D computer
graphics. Perlin constructs patterns from bands of noise, each ofwhich is limited to a range of frequencies. The bands are used asbuilding blocks to construct complex noise patterns, where the frequency spectrum is controlled by weighting the contributions of different bands. For example, fractal noise can be produced by making
the amplitude of each band inversely proportional to its frequency.
Because this approach uses a single weight per band, it does notsupport detailed spectral shaping. By contrast, [Lewis 1989] introduced a Wiener interpolation method that is capable of constructing
a texture with an arbitrary spectrum. This provides more control
over the characteristics of the noise, but is considerably more expensive.
Despite a considerable amount of research into the definition andconstruction of noise (see [Peachey 2003] for an excellent survey),
Perlin’s original version of noise has continued to be the workhorseof the industry. There are good reasons for this: it is fast, it issimple, and the bands provide sufficient spectral control for mostapplications.
Each Perlin band is intended to be band-limited so that it contains
only frequencies in a power-of-2 range. But each band actuallycontains a much wider range of frequencies; this is discussed indetail in [Lewis 1989] and is evident in the Fourier transforms in
Figure 8. These weak band limits lead to some serious problems.
A Perlin band near the Nyquist limit contains both frequencies thatare low enough to be representable (i.e., they contain detail thatshould be in the image) and frequencies that are high enough to be
unrepresentable (i.e, they can cause aliasing). Excluding the band(a)(b)
Figure 1: A comparison between images created using (a) Perlin
noise and (b) wavelet noise. Image (a) represents best practices
use of Perlin noise at Pixar to achieve the optimal tradeoff between
detail and aliasing; notice how much detail is missing at high spatial
frequencies in the far distance.
causes loss-of-detail artifacts, but including it causes aliasing artifacts. Balancing this tradeoff between loss of detail and aliasing
has been a constant source of frustration for shader writers at Pixarand elsewhere. Because aliasing is usually more unacceptable thanloss of detail in feature film production, bands are attenuated aggressively. An unfortunate consequence of this is that as you zoom
into a scene the texture detail becomes visible later than the geometry detail, so the texture doesn’t appear to be tied to its geometry.
Instead, the texture appears to fade in unnaturally as if there werea haze that obscured only some aspects of the surface appearance
and only when it was farther away. This is illustrated in Figure 1. Asparse convolution method was introduced in [Lewis 1989] to improve the band-limited character of noise, but it too does not completely solve the loss-of-detail vs. aliasing problem.
There is, however, an even more fundamental problem — one that
plagues every existing noise creation method. When rendering, it
is common to texture 2D surfaces by sampling a 3D noise function,
but the resulting 2D texture will in general not be band-limited,even if the 3D function is perfectly band-limited. This means thatthe loss-of-detail vs. aliasing tradeoff cannot be solved simply by
constructing a band-limited 3D function. To our knowledge thisproblem has never even been identified, much less addressed.
In this paper, we introduce a new noise technique called wavelet
noise that addresses all of these issues. Our approach is based on
the observation that Perlin’s construction using sums of scaled andattenuated versions of a band-limited function was an early example of what has come to be known as a multiresolution function.
Subsequent to Perlin’s work, wavelet analysis (cf. [Chui 1992;
Stollnitz et al. 1996]), also known as multiresolution analysis, has
emerged as a powerful way to analyze and construct such functions.
It is therefore natural to use this approach to analyze and constructnoise.
Wavelet noise is almost perfectly band-limited, providing good detail with minimal aliasing, as demonstrated in Figure 1. Moreover,
3D wavelet noise can be used to texture a 2D surface in a way thatmaintains its band-limited character. The technique is also easy to
implement and fast (an implementation is provided in the Appendices).

We construct our multiresolution noise M(x) in a fashion similar
to Perlin, that is, by summing scaled and attenuated versions of anoise band N(x). Specifically,1
bmax

Downsample∑ wb N(2b x)

M(x) =(1)b=bminwhere b indexes the band, and where the weights wb are free variables used to control the spectral character of M(x). We differ from
Perlin in the definition and construction of N(x).

(b)

(a)

Upsample-

We begin by modeling how noise is used in the rendering process.
Renderers typically use a filter kernel to turn a scene function into
pixels. In particular, the value of the i-th pixel isZ

Pixel(i) =S(x)K(x − i)dx(2)where S(x) represents the scene being rendered and K(x − i) is therenderer’s filter kernel K(x) translated so that it is centered at pixel i.
K(x) is usually a non-negative, locally supported function.

(d)

(c)

Figure 2: (a). Image R of random noise, (b) Half-size image R↓ , (c)Half-resolution image R↓↑ , (d) Noise band image N = R − R↓↑ .

2

Overview

To provide a context for the remainder of the paper and to emphasize the simplicity of the basic algorithm, we begin by presenting
a high-level overview, using the two-dimensional case and glossing over the details. The essence of our algorithm consists of thefollowing four steps, which are illustrated in Figure 2:
1. Create an image R filled with random noise. (Figure 2a)
2. Downsample R to create the half-size image R↓ . (Figure 2b)
3. Upsample R↓ to a full size image R↓↑ . (Figure 2c)
4. Subtract R↓↑ from the original R to create N. (Figure 2d)
N is thus created by taking R and removing the part that is representable at half-size. What’s left is the part that is not representableat half-size, i.e., the band-limited part. N is then used just like Perlin noise to construct noise patterns. This construction is similarin spirit to the procedural band-pass pyramids of [Perlin and Velho
1995].
The mathematics for deriving and justifying this approach, particularly the filters used in the downsampling and upsampling steps, ispresented in Section 3. Practical implementation considerations arethen discussed in Section 4, and results are presented in Section 5.

3

TheoryIn this section, we provide the mathematical underpinnings for ouralgorithm. For now, we work in one dimension; the extension tomore dimensions is straightforward and is covered in Section 3.6.
Perhaps the most obvious way to construct band-limited noise
would be to use spectral methods, which are based on sinc filters.
Instead, the approach we take uses wavelets, which are more general than spectral methods in that they work with a wide variety ofbasis functions. This flexibility allows us to use basis functions thathave a finite extent and that more closely match the filters used inrenderers.

The scene function S(x) typically contains a mixture of declarative
data, such as the positions and orientations of geometric primitives,
together with procedural components, such as procedural shaders.
Since the scene can be an arbitrarily complicated function of noise,no noise function can possibly prevent aliasing under all conditions.
(As an extreme example, consider a procedural shader that returns
white if the low order bit of M(x) is one, and black otherwise.)
Instead, we devise a method that guarantees no aliasing or loss ofdetail under ideal conditions, and then show that it behaves well inall but the most pathological of non-ideal conditions.

3.1

Orthogonality

Under our ideal conditions, the scene being rendered, when restricted to the region contributing to pixel i, is locally well approximated by a scaled and translated version of multiresolution noise.
If s is the scale of the scene at pixel i, and k is the offset, then ourassumption is that S(x) ≈ M(2s (x − k)), in which case the value ofi-th pixel is given by:Z

Pixel(i)=Z=

M(2s (x − k))K(x − i)dx

(3)

bmax

∑ wb N(2s+b (x − k))K(x − i)dx

(4)

b=bminbmax +s

=

∑

Z

w j−s

N(2 j x − `)K(x − i)dx

(5)

j=bmin +s

where for notational convenience, in the last equation we substituted s + b → j and 2s+b k → `.
Our objective is to allow shader writers to safely truncate this summation without loss of detail and without aliasing. We assume without loss of generality that the resolution of the image corresponds
to j = −1, so that the noise band at resolution j = −1 is fully representable but noise bands at resolutions j ≥ 0 are not representable.
We want these unrepresentable bands to have no contribution to theimage, which means that:Z

N(2 j x − `)K(x − i)dx = 0(6)

for all j ≥ 0 and for all i and `. This condition means that the finescale versions of N(x) are orthogonal to the renderer’s filter kernel
K(x). In other words, our method is based on the principle thatthe contribution of the noise function to the image should vanish
exactly when aliasing would occur.
1 A band-dependent translation is usually added to 2b x to de-correlate thenoise bands, but we can safely ignore this in our analysis.

We know of no technique that can find such an N(x) for an arbitrary
filter kernel K(x) and arbitrary parameters j and `. However, given
relatively mild restrictions on K(x), we can use wavelet analysis toachieve exact orthogonality when j and ` are integers. With thisobjective in mind, Sections 3.2 and 3.3 review some needed resultsfrom wavelet analysis2 in preparation for constructing our orthogonal noise band N(x) in Section 3.4. We then show in Section 3.5
that the deviation of our noise function from orthogonality is small
in real-world situations.

3.2

Refinability and Upsampling

Much of computer graphics involves representing functions as aweighted sum of basis functions. For example, B-spline curves
are represented as weighted sums of B-spline basis functions, anda monitor approximates an image as a weighted sum of the spot
shapes of the pixels.
Starting with a basis function φ (x) centered at x = 0, other functions
F(x) can be built by taking linear combinations of φ (x) translated
by integer amounts i:
F(x) = ∑ fi φ (x − i)(7)iwhere φ (x − i) is a version of φ (x) centered around x = i, and where
the fi are the coefficients of the representation.
The set of all functions that can be represented by varying the coefficients fi forms a vector space that we call the resolution 0 spaceand denote by S0 :S0 := {F(x)|F(x) = ∑ fi φ (x − i)}(8)i

A larger, resolution 1 space S1 of functions can be represented byscaling down the width of φ (x) by a factor of two (21 ), then translating along the half integers. Functions G(x) in this space take theform
G(x) = ∑ gi φ (2x − i)(9)Figure 3: Refinement of a uniform quadratic B-spline. The black
curve is the basis function φ (x). The red curves are φ (2x − k),integer-offset half-width copies of φ (x) that have been attenuatedby the refinement coefficient pk . The blue dots show the pk values (1/4, 3/4, 3/4, 1/4) for this basis. This basis function is refinable because the sum of the red curves is exactly equal to the black
curve.
The sequence (..., fi↑ , ...) is twice as long as the sequence(..., fk , ...), so Equation 12 represents an upsampling filter; this is
in fact the upsampling filter used in Step 3 of our algorithm. Thisfilter, defined solely by the refinement coefficients and hence by
the basis function, describes how to convert a function represented
at resolution 0 to its resolution 1 representation. For quadratic Bsplines, all of the refinement coefficients pk are 0 except for:p−1 = 1/4, p0 = 3/4 p1 = 3/4, p2 = 1/4
(13)

3.3

Wavelets and Downsampling

Above we established that refinability guarantees that every member of S0 can be represented exactly as a member of the higherresolution space S1 . The converse, however, is not true: not every function G(x) in S1 can be represented exactly in the lowerresolution space S0 ; in general there is some lost detail. However,this loss of detail can be minimized in a least squares sense by dividing G(x) into two parts:
G(x) = G↓ (x) + D(x)
(14)
where
• G↓ (x) is the least squares best approximation to G(x) in S0 .

i

A key idea in wavelet analysis is to ensure that all the functions
in S0 are also in S1 ; that is, that S0 ⊂ S1 . This guarantees that S1enriches the space S0 , rather than replacing it with a completely
different set of functions. This subset relationship is guaranteed if
there exist coefficients pk that allow φ (x) to be written in terms ofφ (2x − k):
φ (x) = ∑ pk φ (2x − k)
(10)k

When such coefficients do exist, φ (x) is said to be refinable. Examples of refinable functions include uniform B-spline basis functions(of any degree), sinc functions, and Daubeschies basis functions.
Figure 3 illustrates the refinability of uniform quadratic B-splines.
Gaussians are good examples of functions that are not refinable.
We henceforth assume that φ (x) is refinable. Given the coefficientsfi that represent a function F(x) in S0 as in Equation 7, there exist
coefficients fi↑ that represent that function exactly in S1 :F(x) = ∑ fi↑ φ (2x − i)

(11)

i

Using the refinement equation, Equation 10, it is straightforward toshow thatfi↑ = ∑ pi−2k fk(12)k
2 A number of texts, such as [Chui 1992] and [Stollnitz et al. 1996], de-scribe this material in more detail.

• D(x) is the least squares residual and contains the informationin S1 that cannot be represented in S0 .
It follows that D(x) and all of its integer translates are orthogonal
to all functions in S0Z. That is,D(x − `)φ (x − i)dx = 0(15)

for all integers ` and i. The set of all such detail functions D(x)forms a vector space called the wavelet space W0 .
Given the coefficients gi for G(x), the coefficients g↓i for G↓ (x)can be determined using a standard result from wavelet analysis.
Specifically, there exist so-called analysis coefficients ak , again depending only on the basis function φ (x), such thatg↓i = ∑ ak−2i gk(16)kThis is a downsampling filter that projects the coefficient sequence(..., gk , ...) into a sequence (..., g↓i , ...) of half the length in a leastsquares optimal way. The analysis coefficients ak comprise thedownsampling filter of Step 2 in our algorithm. For quadratic Bsplines, the sequence (..., ak , ...) is infinitely long, but it decays
quickly and is well approximated by a finite sequence. The variable
a in Appendix 1 contains this sequence for quadratic B-splines.
Because φ (x) is refinable, it can be represented exactly in terms ofφ (2 j x) for all integers j ≥ 0. From this and Equation 15, it canreadily be shown thatZ
D(2 j x − `)φ (x − i)dx = 0(17)

for all integers j ≥ 0 and all integers ` and i. This is exactly thecondition required for our noise band N(x) in Equation 6. Thus, if
we build our noise bands in the wavelet space W0 , then they can bescaled to any resolution j and be guaranteed to have no effect on
images at any resolution less than j.

3.4

(a)

(b)

Constructing noise bands

Now we are ready to construct our noise band N(x) so that it lies in
the wavelet space W0 generated by the filter kernel K(x). At thispoint, we must choose a particular filter kernel K(x), which we taketo be the uniform quadratic B-spline basis function B(x). We do sofor several reasons:• It is a good approximation to the filter kernels used by mostrenderers (including Gaussians);• it has a small support and low degree, which makes the evaluation of N(x) simple and computationally efficient; and• it generates differentiable noise, meaning that noise gradients,
which are also useful in procedural shading, are continuous.

(c)(d)
Figure 4: An example of wavelet noise in one dimension: (a) afunction R(x) built using random coefficients with B-spline basis
functions; (b) the downsampled version R↓ (x); (c) the re-upsampled
version R↓↑ (x); (d) the noise function N(x) = R(x) − R↓↑ (x). Functions appear in black, their coefficients are plotted as blue dots, andthe basis functions multiplied by their corresponding coefficients
appear in red. Note that the basis functions in (b) are twice thewidth of those in the other figures, and that the function R↓ (x) in(b) and the function R↓↑ (x) in (c) are identical.
possible because by upsampling we can represent the S0 partR↓ (x) exactly in S1 as:

One could, however, use this same technique to create noise for any
other refinable basis function, such as uniform cubic B-splines orsincs.

R↓↑ (x) = ∑ ri↓↑ B(2x − i)

(24)

i

What follows here is a more detailed and precise version of thealgorithm presented in Section 2. Code to implement this algorithm(in three dimensions) is provided in Appendix 1.

where the coefficients r↓↑ may be computed by upsampling
the coefficients r↓ using Equation 12:
ri↓↑ = ∑ pi−2k rk↓

1. Create R(x) using a random number generator.
Use a random number generator to create a coefficient sequence R = (..., ri , ...). This sequence defines a function R(x)in S1 . Using the basis function B(x), Equation 9 becomes:
R(x) = ∑ ri B(2x − i)

(25)

k

An example of R↓↑ (x) is shown in Figure 4(c).
4. Compute N(x) by subtracting R↓↑ (x) from R(x).
Using Equation 24, we can now rewrite Equation 23 as:

(18)

i

N(x)

An example of such a function is shown in Figure 4(a).

=∑ ri B(2x − i) − ∑ ri↓↑ B(2x − i)(26)

∑ ni B(2x − i)

(27)

i

2. Compute R↓ (x) by downsampling R(x).

=

i

i

As shown in Equation 14, R(x) can be decomposed into twoparts: a part in S0 , denoted by R↓ (x), and a detail part in W0 ,which is the function N(x) whose coefficients we seek:R(x) = R↓ (x) + N(x)↓N(x) = R(x) − R (x)

(19)
(20)

Using Equation 7, we get:R↓ (x) = ∑ ri↓ B(x − i)

(21)

i

The coefficients r↓ may be computed from the coefficients rusing Equation 16:ri↓ = ∑ ak−2i rk

(22)

k

An example of R↓ (x) is shown in Figure 4(b).
3. Compute R↓↑ (x) by upsampling R↓ (x).
Using Equations 18 and 21, we can express Equation 20 as:
N(x) = ∑ ri B(2x − i) − ∑ ri↓ B(x − i)
i

(23)

i

The first sum in this equation is an S1 representation that uses
the basis function B(2x − i); the second sum is an S0 representation that uses the basis function B(x − i). To combine
these sums, we need them to have a common basis, which iswhere the coefficients of N(x) are given by ni = ri − ri↓↑ . An
example of the resulting noise band N(x) is shown in Figure 4(d).
Once the coefficients ni have been determined, a value of N(x) for agiven x can be computed using any evaluation method for quadratic
B-splines (cf. [Farin 2002]). The code in Appendix 2 shows thiscomputation in 3 dimensions.
As established in Section 3.3, this construction guarantees that all
bands of noise N(2 j x) for j ≥ 0 are orthogonal to S0 , and hence do
not contribute to the image. Additionally, noise bands are orthogonal to each other, which makes spectral shaping more controllable.

3.5

Fractional scales and translates

The noise band N(x) has been constructed so that N(2 j x − `) contributes nothing to pixel values when ` is an arbitrary integer andwhen j is a non-negative integer. However, in practice j and ` arenot always integers, and in these cases noise bands for which j ≥ 0will have some contribution to the pixel; omitting them from thesummation will cause some detail to be lost. Fortunately the contribution falls off very rapidly as j increases.
As an indication of this falloff, we have numerically studied thecontribution functionZ
C( j, `) :=N(2 j x − `)B(x)dx

(28)

-2

-1

0

1

2

3

Figure 5: The contribution function C( j, 1/2) for −2 ≤ j ≤ 3.
as a continuous function of j and `, using various instances of ournoise band N(x). For fixed j, the function reaches a maximum at` = 1/2. The rate of falloff of the maximum contribution can therefore be seen by plotting C( j, 1/2) as function of j, as is done in
Figure 5. The falloff is rapid, so we can safely truncate the summation at j = 0 with minimal loss of detail. (In practice, to ensure thatthe contribution completely vanishes at j = 0 even for non-integer`, we fade out the noise as we approach j = 0.)

3.6

Noise in multiple dimensions

The two-dimensional noise function N(x, y) is constructed from a
two-dimensional array nix ,iy of coefficients that are blended together
to form a uniform quadratic tensor-product B-spline:(29)
N(x, y) = ∑ nix ,iy B(2x − ix )B(2y − iy )ix ,iyThe coefficient array N = (..., nix ,iy , ...) is obtained from an arbitrary array R = (..., rix ,iy , ...) by applying the downsampling filter
and then the upsampling filter to each row of R to obtain an arrayR↓↑ . Next, the downsampling filter followed by the upsampling filter is applied to each of the columns of R↓↑ to obtain an array ↑↓ R↓↑ .
This array is then subtracted from the initial array R to produce thenoise coefficient array N. Generation of coefficients for noise inthree and higher dimensions is accomplished in a similar fashion,
as shown in the code in Appendix 1.
In analogy to the orthogonality condition for one-dimensional
noise, the two-dimensional noise function N(x, y) satisfies the twodimensional orthogonality condition:
Z

N(2 j x, 2 j y)B(x − ix )B(y − iy )dx dy = 0,

j≥0

(30)

meaning that bands of noise for j ≥ 0 can be safely truncated.

N(2 j x, 2 j y, 2 j z)B(x − ix )B(y − iy )B(z − iz )dx dy dz = 0(31)

for j ≥ 0, meaning that bands of noise for j ≥ 0 can be safely truncated. This integral is relevant when 3D noise is used to texture a
volume.

3.7

slice is obtained by integrating the Fourier transform along the direction perpendicular to the slice. For example, the 1D Fouriertransform of the slice shown in Figure 6(a) is obtained by integrating the 2D Fourier transform shown in Figure 6(b) along lines
perpendicular to each point on the white line. Notice that the slice
can contain arbitrarily low frequencies because the low-frequency
line integrals (i.e., near the origin of the Fourier transform) pick up
contributions from non-zero high-frequency parts of the 2D Fouriertransform. For example, the integral along the dotted line, which
determines a low frequency, is in general non-zero since it passes
through the higher frequencies in the band-pass.
Fortunately, we can use the wavelet orthogonality machinery to address this issue. Consider first the case of a noise-axis-aligned planar slice such as the plane 2 j z = z0 . Rather than texturing the planeat (x, y) with the 3D function N(2 j x, 2 j y, z0 ), we instead use the 2D
function Nz0 (2 j x, 2 j y), whereZ
Nz0 (2 j x, 2 j y) :=N(2 j x, 2 j y, 2 j z)B(z − z0 )dz

Projected noise

Texture values on a 2D surface are frequently determined by sampling a 3D noise function. This approach runs into a fundamental
problem, however: a 2D slice through a 3D band-limited function
is in general not band-limited. As far as we are aware, this fact does
not seem to have been previously recognized in the texture synthesis literature.3 As an example, Figures 8(a) and (b) show how 3D
Perlin noise is somewhat band-limited, but a planar slice through it
is not band-limited at all.
To see why this is so, we use the Fourier Slice Theorem (cf.
[Malzbender 1993]), which states that the Fourier transform of a
3 [Lewis 1989] observed this phenomenon in a 1D slice through 2D Per-lin noise but attributed it to Perlin’s interpolation method rather than theslicing process itself.

(32)

This shows how to project wavelet noise onto a surface: instead ofsimply point sampling the texture at the intersection point, we perform a line integral orthogonal to the surface, where the integrandis the 3D noise weighted by a 1D filter kernel with its center at the
point of intersection. Using Equation 31, this definition guarantees
the orthogonality
condition we desire:Z
Nz0 (2 j x, 2 j y)B(x − ix )B(y − iy )dx dy = 0

Similarly, in three dimensions, N(x, y, z) satisfies
Z

(a)
(b)
Figure 6: A band-limited function (a) and its Fourier transform (b).
The Fourier transform of the white slice in (a) is obtained by integrating the 2D Fourier transform along the direction perpendicular
to the white slice shown in (b).

(33)

for all integer z0 and j ≥ 0.
Moreover, Equation 32 is a projection operator from the tri-variate
wavelet space W0 (x, y, z) to the bi-variate wavelet space W0 (x, y).
(This can be established by expanding N(x, y, z) in a wavelet basis for W0 (x, y, z).) As a consequence, the other important property, band orthogonality, is also preserved: i.e., Nz0 (2i x, 2i y) andNz0 (2 j x, 2 j y) are orthogonal whenever i 6= j.
In some cases, this axis-aligned projection is just what we want.
For example, when projecting 4D noise (indexed by xyz and time)
onto 3D images (xyz) in an animated sequence, the projection isalong the time dimension so the integration is axis-aligned. But inmany cases the projection is not axis-aligned, so the integration is
not axis-aligned. In our numerical experiments that consider noninteger j and non-axis alignment, we have found that the deviationfrom orthogonality is small in the non-axis-aligned case and falls
off rapidly in j in a fashion similar to that shown in Figure 5.

4

Implementation

In this section we examine a variety of issues that arise in practicalimplementations of wavelet noise.

4.1

Noise evaluation

In practice, renderers approximate the integral of Equation 2 using
a quadrature formula. In the case where S(x) = N(x) this becomesQ

Pixel(i) ≈ ∑ Bq N(xq )

(34)

q=1

where the quadrature weights Bq are typically just set to B(xq ), butthey could also be weights that provide higher order accuracy, using for instance Gaussian quadrature. If the renderer’s quadratureformula is sufficiently accurate, we can use N(xq ) directly. N(xq ) is
calculated by interpolating the coefficients ni around x = xq using
B(x − xq ), which is 3 noise coefficients wide.
However, if the number of samples Q is too small, a more accurate
value for the pixel can be obtained by replacing N(xq ) by N̄(xq ), aweighted average of the value of N(x) in the vicinity of xq . In thelimit of one sample (Q = 1), we are setting Pixel(i) = N̄(i), so N̄(i)would ideally be equal to the value of integral in Equation 2. The
optimal computation of N̄(xq ) depends on the details of the quadrature calculation of the renderer, but a reasonable and convenient
approximation is to use a widened version of B(x) for the finest
bands. For example, for Q = 1, we have had good success interpolating the noise coefficients using B((1 + 2 j )(x − xq )), which is 6
noise coefficients wide at j = 0 and 4.5 wide at j = −1, and thengradually unwidening the basis function in the region −2 < j ≤ −1
so that for j ≤ −2, we use the regular, unwidened B(x − xq ).
The projected noise integral in Equation 32 can be expressed interms of the noise coefficients n as (taking j = 0 for simplicity):Nz0 (x, y) = ∑ nix iy iz B(2x − ix )B(2y − iy )β (2z0 − iz )(35)

Figure 7: The quintic projection function β (z) (in blue) can beclosely approximated by a quadratic basis function (in red).
Thus if the ri are Gaussian, the ni will be too. There are various
standard ways to generate such values ri ; we use the polar formof the Box-Muller transformation (cf. [Knuth 1997]) which gives
us a Gaussian distribution with a mean of 0 and a variance of 1.
The variance σn2 of the ni can then be determined analytically using
the ri weights that are equivalent to the upsample, downsample,
subtract process; however, it is simpler to calculate σn2 numericallyfrom the ni values.
N(x) is computed as a linear combination of the ni using a set ofbasis functions evaluated at x. As a result, the variance σN2 of N(x)differs from σn2 by the sum of squares of evaluated basis functions.
We could use a different variance for every x, but we have foundit preferable to use the average variance. The former allows a bit
more control over the variance; the latter preserves orthogonality.
The average variance can be determined analytically, but once again
it is simpler to evaluate it numerically. For quadratic B-splines, thisaverage σN2 is 0.265 for 2D noise, 0.210 for 3D noise, and 0.296 for
3D noise projected onto a 2D surface.
As bands of noise are added, the variance is further changed by theband weights wb . The variance of the final multiresolution noise
M(x) is:ix iy iz

2
σM
= σN2

where

Z

β (z) =

B(z0 )B(2z0 − z)dz0

4.2

Noise variance

Another issue with Perlin noise is that there is no model for the
expected statistical distribution of the intensity values or how that
distribution changes as a function of the number of bands and theirweights. By contrast, as we show in this section, the statistical
distribution of wavelet noise can be calculated and controlled.
Every step leading from the random variables ri to the noise coefficients ni is linear, meaning that the noise coefficients ni are linear
combinations of the ri . A standard result from statistics says that a
linear combination of independent, identically distributed Gaussian
random variables is again a Gaussian random variable. Specifically,
let x1 , ..., xn be Gaussian random variables with variance σx2 . The
variableny = ∑ w i xi

(37)

i=1

has variancenσy2 = σx2 ∑ w2i
i=1

(38)

(39)

j= jmin(36)

is a piecewise quintic whose support covers 9 noise coefficients.
z+3/2
This quintic can be closely approximated by β (z) ≈ B( 2 ),which is a double-width quadratic B-spline, whose support covers
6 coefficients as shown in Figure 7. This allows us to combine projection and evaluation using a single asymmetric tri-variate basis
function. Note that if the normal is not parallel to one of the axesof the noise coefficient grid, the computation will be less efficientbecause the support of the basis function will not be aligned withthe grid. Code for this is included in Appendix 2.

0

∑ w2j

Because the values in M(x) have a Gaussian distribution with a predictable variance, we can change the values to match any desired
distribution. To match a Gaussian distribution with a different variance, we simply scale the results; for a white distribution, we use
the erf function (cf. [Abramowitz and Stegun 1970]); otherwise webuild a lookup table. This control is useful, but making the distribution non-Gaussian does violate the orthogonality assumption. Like
any extreme transformations applied to the noise, abrupt transitions
within the desired distribution can lead to aliasing during animation, but we have found that smoothly varying distributions do not
introduce any noticeable artifacts.

4.3

Noise tiles

The algorithm as described so far is not very efficient because thewidth of the downsampling filter makes the generation of the noise
coefficients expensive. Instead of incurring this expense for every
noise sample during rendering, we can use a small pre-computed
volume of noise coefficients and then tile our space with that volume.4 The tile can optionally include a border the width of thebasis function so that tile boundary checks can be avoided in theevaluation loop.
The downside of using noise tiles is that they introduce lowfrequency repeating patterns into the noise. Fortunately, because
the noise band contains only small-scale features, this large-scale
repetition is often not objectionable and can just be ignored. For
some applications, it may be appropriate to use a larger tile to make4 Although it isn’t immediately obvious, the tile size must be even; thisis because the period of R↓↑ (x) is equal to the period of R(x) if the size iseven, but twice the period of R(x) if it is odd.

(a) 2d Perlin noise(d) 2d wavelet noise(b) 2D slice through 3D Perlin noise(e) 2D slice through 3D wavelet noise(c) 2D white noise(f) 3D wavelet noise projected onto 2D

Figure 8: Noise patterns (left) with their Fourier transforms (right). For Perlin noise, we use the RenderMan implementation of [Perlin 2002].
the repetition less obvious. For other applications, the repeating
pattern can be eliminated by using a hash function to select a different noise tile to use in each tile location. This strategy avoids
repetition but creates discontinuities at the tile boundaries. These
discontinuities can be avoided by overlapping the tiles slightly andblending from one tile to the other, but this creates another artifact:
blending reduces the variance of the noise as noted in Section 4.2.
Using Equation 38, if the tiles are blended using α of one tile and1 − α of another, the variance of the blended value is reduced by afactor of α 2 + (1 − α)2 . This appears as areas of reduced contrast
along the tile borders. We can correct this simply by dividing theblended tile values by the variance reduction. Since the α values
are a function of the location within the tile, we can pre-multiply
the noise coefficients by the variance-corrected α values at tile creation time, eliminating the run-time cost of the α multiplicationand variance correction. To reduce storage, we can store a single
tile pattern and create 48 different variations by permuting the x, y,
and z tile indices (6 possibilities) and by stepping through the coefficients in reverse order in x, y, or z (2x2x2 possibilities). Thisrequires storing only 4 numbers per variation: step sizes in x, y, andz and a pointer to the (0,0,0) location.
Finally, we note that the upsampling and downsampling filters treat
the even and odd indexed coefficients differently, which introduces

a subtle variation in variance. This variation is easily eliminatedby adding two noise tiles offset by an odd number of coefficientsin each dimension. This can be done at tile creation time so that
there is no additional run-time cost. This was done prior to thecomputation of the variances listed in Section 4.2.

5

Results

5.1

Fourier transforms

A comparison of Perlin and wavelet noise patterns and their Fouriertransforms is shown in Figure 8. The center of each Fourier transform corresponds to the DC term, with x frequencies increasing to
the right and y frequencies increasing vertically. A bright spot at
x = 0.25, y = 0, for example, would indicate a signal that repeats
every 4 pixels in the x direction. The Fourier transform of a noiseband should be dark where |x| and |y| are both < 0.5, with no otherpatterns. Wavelet noise (Figure 8(d)) has these characteristics. Perlin noise (Figure 8(a)), by contrast, contains significant energy inthe low frequency areas near the center of the Fourier transform.
When we look at the Fourier transform of a 2D slice through 3Dnoise, Perlin noise (Figure 8(b)) becomes as non-band-limited as
2D white noise (Figure 8(c)). Wavelet noise (Figure 8(e)) fairs

(a)
(b)
Figure 9: Fourier transforms of 3 bands of 3D noise projected onto
2D with (a) Perlin noise and (b) wavelet noise. In both images,
the middle band (magenta) should nestle exactly between the next
higher band (cyan) and the next lower band (yellow).
better, but low frequencies still leak through. When we project
wavelet noise onto a 2D surface using the technique described in
Section 3.7, however, the band limits are preserved (Figure 8(f)).
Figure 9 shows the Fourier transforms of three adjacent bands ofnoise in different colors and illustrates the band-limiting characterof wavelet noise vs. Perlin noise.
The square nature of the wavelet noise Fourier transform reflects
the asymmetry of our separable quadratic B-spline basis function
and is less evident for non-axis-aligned views. We have not foundthis asymmetry to be objectionable, but it can be reduced by either(1) averaging multiple tiles with different orientations or (2) using a higher B-spline, thus more closely approximating a Gaussian,which is symmetric.

5.2

Run-time cost

The run-time cost of noise depends on many factors of the specificimplementation, but a rough guide to the performance is the number of noise coefficients involved in the interpolation calculation.
The run-time cost of the basic wavelet noise algorithm involves 27(3x3x3) noise coefficients compared to the 24 (2x2x2 triples) usedin Perlin noise. So we would expect the cost of the two techniques
to be similar.
In timing tests on a 1.4 GHz Macintosh G4, however, our moderately optimized implementation of wavelet noise was about 30%
faster than the highly optimized RenderMan implementation of Perlin noise (0.81 µsec vs. 1.06 µsec per noise evaluation). When we
include tile meshing (Section 4.3) with a tile overlap of 1/16 of atile, the cost of wavelet noise increased to 0.97 µsec, which is still
about 10% faster than Perlin noise.
Wavelet noise uses more memory than Perlin noise, but even with arelatively large tile the memory requirements are small by today’s
standards. For example, a 128x128x128 tile uses 8 Mb compared
to about 16 Kb for the RenderMan implementation of Perlin noise.
Using the optional parts of our technique adds additional cost. Projecting 3D noise onto a 2D surface (Sections 3.7 and 4.1) approximately doubles the cost in the axis-aligned case; the non-axisaligned case is more expensive because the interpolation is less efficient. Basis function widening (Section 4.1) approximately triples
the cost of the finest band, doubles the cost of the second band, anddoesn’t affect coarser bands. The decision of which of these options
to use will depend on the application.

5.3

ExamplesWe have embedded our implementation of wavelet noise in a 2D
testbed with controls for shaping the weights of the bands and for(a)(b)
Figure 10: 2D noise patterns with (a) 12 bands with a Gaussian
distribution and (b) 8 bands with a white distribution. The blue bars
are the band weights.
specifying the desired distribution of the result. Figure 10 shows
two patterns produced by this program. Because of the independence of our noise bands, the spectral shape is accurately controlled
by the band weights so that we have not found it necessary to provide a lacunarity control. The disadvantage of values of lacunarity
less than 2 is that they make the bands overlap and lose their mutual
orthogonality.
We have also embedded our implementation of wavelet noise in a
RenderMan shader DSO. Figures 1(a) and 1(b) are stills from two
animations that compare Perlin and wavelet noise. The geometry is
the same in both animations; they differ only in the version of noise
used. The level of noise used in Figure 1(a) shows the amount ofdetail obtainable using the RenderMan version of Perlin noise. Figure 1(b) shows the same scene rendered using wavelet noise using
the same noise band weights.

5.4

Infinite noiseBand weights usually correspond to fixed scales in world space.
Another option we have found useful is to make the weights relative to each location’s scale in screen space. This produces infinite
noise, the ability to zoom into a scene indefinitely, with low frequencies gradually fading out as higher frequencies gradually fadein. We have found that this trick works remarkably well and thatthe fading out is not noticeable if it’s done over a few bands (asin Figure 10(a), for example). We speculate that this is because it
echoes the way human perception adapts to the average intensity ofa scene.

6

Summary

The wavelet noise preprocessing calculations are:
• Create a tile of noise coefficients by filling the tile with random noise, then downsampling, upsampling, and subtracting(Section 3.4 and Appendix 1).
• (Optional) Add two such tiles together to correct for even vs.
odd variance (Section 4.3 and Appendix 1).
• (Optional) Adjust the tile borders to accommodate tile meshing (Section 4.3).
The run-time calculations for each band b, given a location and an
object scale s, are:• Determine the noise resolution j from s and b (Section 3.1 andAppendix 2).
• (Optional) Use a hash function to determine which tile pattern(s) to use for this location (Section 4.3).

• (Optional) If the renderer undersamples, widen the basis function in the tangent directions in the last bands (Section 4.1).
• (Optional) Project the 3D noise onto the 2D surface by doubling the support of the basis function in the direction normal
to the surface (Section 3.7 and Appendix 2).
• Evaluate the noise function by multiplying the noise coefficients by the basis function (Section 3.4 and Appendix 2).
• (Optional) Correct the noise for the desired distribution (Section 4.2 and Appendix 2).
In conclusion, we have shown that wavelets are well suited to producing noise for use in procedural textures. Because the bands
are orthogonal, they provide a set of independent controls over the
shape of the spectrum. The distribution of the final result can bepredicted and controlled. Most importantly, the noise is truly bandlimited, so that virtually all of the detail can be rendered with minimal aliasing, even when projecting 3D noise onto a 2D surface.

Acknowledgements
Thanks to Chris Bernardi and Eben Ostby for pointing out the problems with Perlin noise and getting us to look in this area in the firstplace. Loren Carpenter and Michael Kass helped with the Fourier
analysis.

References
A BRAMOWITZ , M., AND S TEGUN , I. A. 1970. Error function andfresnel integrals. In Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables. Dover, ch. 7.1.1,
297.
C HUI , C. K. 1992. An introduction to wavelets. Academic Press
Professional, Inc., San Diego, CA, USA.
FARIN , G. 2002. Curves and surfaces for CAGD: a practical guide,
5th ed. Morgan Kaufmann Publishers Inc.
K NUTH , D. E. 1997. The Art of Computer Programming, third ed.,vol. 2. Addison-Wesley.
L EWIS , J. P. 1989. Algorithms for solid noise synthesis. In
SIGGRAPH ’89: Proceedings of the 16th annual conference on
Computer graphics and interactive techniques, ACM Press, New
York, NY, USA, 263–270.
M ALZBENDER , T. 1993. Fourier volume rendering. ACM Transaction on Graphics 12, 3 (July), 233–250.
P EACHEY, D. 2003. Building procedural textures. In Texturing andModeling: A Procedural Approach, third ed. Morgan Kaufmann
Publishers Inc., ch. 2.
P ERLIN , K., AND V ELHO , L. 1995. Live paint: painting with procedural multiscale textures. In SIGGRAPH ’95: Proceedings ofthe 22nd annual conference on Computer graphics and interactive techniques, ACM Press, New York, NY, USA, 153–160.
P ERLIN , K. 1985. An image synthesizer. In SIGGRAPH ’85: Proceedings of the 12th annual conference on Computer graphics
and interactive techniques, ACM Press, New York, NY, USA,
287–296.
P ERLIN , K. 2002. Improving noise. In SIGGRAPH ’02: Proceedings of the 29th annual conference on Computer graphics
and interactive techniques, ACM Press, New York, NY, USA,
681–682.
S TOLLNITZ , E., D E ROSE , T., AND S ALESIN , D. 1996. Wavelets
for Computer Graphics. Morgan Kaufmann Publishers.

Appendix 1. Noise tile generation/* Note: this code is designed for brevity, not efficiency; many operations can be hoisted,* precomputed, or vectorized. Some of the straightforward details, such as tile meshing,* decorrelating bands and fading out the last band, are omitted in the interest of space.*/static float *noiseTileData; static int noiseTileSize;
int Mod(int x, int n) {int m=x%n; return (m<0) ? m+n : m;}#define ARAD 16
void Downsample (float *from, float *to, int n, int stride ) {float *a, aCoeffs[2*ARAD] = {0.000334,-0.001528, 0.000410, 0.003545,-0.000938,-0.008233, 0.002172, 0.019120,-0.005040,-0.044412, 0.011655, 0.103311,-0.025936,-0.243780, 0.033979, 0.655340,
0.655340, 0.033979,-0.243780,-0.025936, 0.103311, 0.011655,-0.044412,-0.005040,
0.019120, 0.002172,-0.008233,-0.000938, 0.003546, 0.000410,-0.001528, 0.000334};
a = &aCoeffs[ARAD];for (int i=0; i<n/2; i++) {to[i*stride] = 0;
for (int k=2*i-ARAD; k<=2*i+ARAD; k++)to[i*stride] += a[k-2*i] * from[Mod(k,n)*stride];}}void Upsample( float *from, float *to, int n, int stride) {float *p, pCoeffs[4] = { 0.25, 0.75, 0.75, 0.25 };p = &pCoeffs[2];for (int i=0; i<n; i++) {to[i*stride] = 0;
for (int k=i/2; k<=i/2+1; k++)to[i*stride] += p[i-2*k] * from[Mod(k,n/2)*stride];}}void GenerateNoiseTile( int n, int olap) {if (n%2) n++;/* tile size must be even */int ix, iy, iz, i, sz=n*n*n*sizeof(float);
float *temp1=(float *)malloc(sz),*temp2=(float *)malloc(sz),*noise=(float *)malloc(sz);/* Step 1. Fill the tile with random numbers in the range -1 to 1. */for (i=0; i<n*n*n; i++) noise[i] = gaussianNoise();/* Steps 2 and 3. Downsample and upsample the tile */for (iy=0; iy<n; iy++) for (iz=0; iz<n; iz++) {/* each x row */i = iy*n + iz*n*n; Downsample( &noise[i], &temp1[i], n, 1 );Upsample(&temp1[i], &temp2[i], n, 1 );}for (ix=0; ix<n; ix++) for (iz=0; iz<n; iz++) {/* each y row */i = ix+ iz*n*n; Downsample( &temp2[i], &temp1[i], n, n );Upsample(&temp1[i], &temp2[i], n, n );}for (ix=0; ix<n; ix++) for (iy=0; iy<n; iy++) {/* each z row */i = ix+ iy*n;Downsample( &temp2[i], &temp1[i], n, n*n );Upsample(&temp1[i], &temp2[i], n, n*n );}/* Step 4. Subtract out the coarse-scale contribution */for (i=0; i<n*n*n; i++) {noise[i]-=temp2[i];}/* Avoid even/odd variance difference by adding odd-offset version of noise to itself.*/
int offset=n/2; if (offset%2==0) offset++;for (i=0,ix=0; ix<n; ix++) for (iy=0; iy<n; iy++) for (iz=0; iz<n; iz++)
temp1[i++] = noise[ Mod(ix+offset,n) + Mod(iy+offset,n)*n + Mod(iz+offset,n)*n*n ];for (i=0; i<n*n*n; i++) {noise[i]+=temp1[i];}noiseTileData=noise; noiseTileSize=n; free(temp1); free(temp2);}

Appendix 2. Noise evaluation
float WNoise( float p[3]) {/* Non-projected 3D noise */int i, f[3], c[3], mid[3], n=noiseTileSize;/* f, c = filter, noise coeff indices */float w[3][3], t, result =0;/* Evaluate quadratic B-spline basis functions */for (i=0; i<3; i++) {mid[i]=ceil(p[i]-0.5); t=mid[i]-(p[i]-0.5);w[i][0]=t*t/2; w[i][2]=(1-t)*(1-t)/2; w[i][1]=1-w[i][0]-w[i][2]; }/* Evaluate noise by weighting noise coefficients by basis function values */for(f[2]=-1;f[2]<=1;f[2]++) for(f[1]=-1;f[1]<=1;f[1]++) for(f[0]=-1;f[0]<=1;f[0]++) {float weight=1;
for (i=0; i<3; i++) {c[i]=Mod(mid[i]+f[i],n); weight*=w[i][f[i]+1];}
result += weight * noiseTileData[c[2]*n*n+c[1]*n+c[0]]; }return result;}float WProjectedNoise( float p[3], float normal[3]) {/* 3D noise projected onto 2D */int i, c[3], min[3], max[3], n=noiseTileSize;/* c = noise coeff location */float support, result=0;/* Bound the support of the basis functions for this projection direction */for (i=0; i<3; i++) {support = 3*abs(normal[i]) + 3*sqrt((1-normal[i]*normal[i])/2);
min[i] = ceil( p[i] - (3*abs(normal[i]) + 3*sqrt((1-normal[i]*normal[i])/2)) );max[i] = floor( p[i] + (3*abs(normal[i]) + 3*sqrt((1-normal[i]*normal[i])/2)) ); }/* Loop over the noise coefficients within the bound. */for(c[2]=min[2];c[2]<=max[2];c[2]++) {for(c[1]=min[1];c[1]<=max[1];c[1]++) {for(c[0]=min[0];c[0]<=max[0];c[0]++) {float t, t1, t2, t3, dot=0, weight=1;/* Dot the normal with the vector from c to p */for (i=0; i<3; i++) {dot+=normal[i]*(p[i]-c[i]);}/* Evaluate the basis function at c moved halfway to p along the normal. */for (i=0; i<3; i++) {t = (c[i]+normal[i]*dot/2)-(p[i]-1.5); t1=t-1; t2=2-t; t3=3-t;
weight*=(t<=0||t>=3)? 0 : (t<1) ? t*t/2 : (t<2)? 1-(t1*t1+t2*t2)/2 : t3*t3/2;}/* Evaluate noise by weighting noise coefficients by basis function values. */result += weight * noiseTileData[Mod(c[2],n)*n*n+Mod(c[1],n)*n+Mod(c[0],n)];}}}return result;}float WMultibandNoise( float p[3],float s,float *normal,int firstBand,int nbands,float *w) {float q[3], result=0, variance=0; int i, b;for (b=0; b<nbands && s+firstBand+b<0; b++) {for (i=0; i<=2; i++) {q[i]=2*p[i]*pow(2,firstBand+b);}
result += (normal) ? w[b] * WProjectedNoise(q,normal) : w[b] * WNoise(q);}for (b=0; b<nbands; b++) {variance+=w[b]*w[b];}/* Adjust the noise so it has a variance of 1. */if (variance) result /= sqrt(variance * ((normal) ? 0.296 : 0.210));return result;}