Sparse PDF Maps for Non-Linear Multi-Resolution Image Operations
Markus Hadwiger
KAUST

Ronell Sicat
KAUST

Johanna Beyer
KAUST

Jens Krüger
IVDA, DFKI, Intel VCITorsten Möller
Simon Fraser University

Figure 1: sPDF-maps are a compact multi-resolution image pyramid data structure that sparsely encodes pre-computed pixel neighborhood
probability density functions (pdfs) for all pixels in the pyramid. They enable the accurate, anti-aliased evaluation of non-linear image
operators directly at any output resolution. A variety of operators can be computed at run time from the same pre-computed data structure in a
way that scales to gigapixel images, such as local Laplacian filters for (b,d) detail enhancement or (c,e) smoothing, (f) median filters, (g)dominant mode filters, (h) maximum mode filters, (i) bilateral filters. The original image (a) has resolution 16, 898 × 14, 824 (250 Mpixels).

Abstract

1

We introduce a new type of multi-resolution image pyramid forhigh-resolution images called sparse pdf maps (sPDF-maps). Eachpyramid level consists of a sparse encoding of continuous probability
density functions (pdfs) of pixel neighborhoods in the original image.
The encoded pdfs enable the accurate computation of non-linear image operations directly in any pyramid level with proper pre-filtering
for anti-aliasing, without accessing higher or lower resolutions. The
sparsity of sPDF-maps makes them feasible for gigapixel images,
while enabling direct evaluation of a variety of non-linear operatorsfrom the same representation. We illustrate this versatility for antialiased color mapping, O(n) local Laplacian filters, smoothed local
histogram filters (e.g., median or mode filters), and bilateral filters.

The recent increase in the resolution of acquired and computed image
data has resulted in a need for novel multi-resolution techniques, e.g.,
representing, processing, and rendering gigapixel images [Kopf et al.
2007; Summa et al. 2010]. One of the most prevalent types of multiresolution image hierarchies are pyramid representations, such as
mipmaps [Williams 1983], or Gaussian pyramids [Burt and Adelson
1983]. Image pyramids store pre-filtered and downsampled versionsof the original image, where the pre-filtering is crucial for avoiding
aliasing. Because standard pre-filters are linear operators, further
linear operators (e.g., further smoothing) can be applied accurately
in any pre-computed pyramid level, as they commute with the prefilter. However, this does not apply to non-linear operators, suchas color mapping or edge-preserving filters. In this case, one has
to process the original image and re-compute the pyramid. This is
impractical for gigapixel images, especially when the goal is theinteractive display of processed images at a lower output resolution,
which is becoming more important as the gap between the acquired
image and display resolutions widens. Therefore, in such a scenario
all operators are usually applied directly at the output resolution.
Nevertheless, downsampling the image first using a standard prefilter, e.g., bicubic, followed by a non-linear operation, can introduce
false colors. These artifacts due to resampling are a form of aliasing,though not of the same kind as the well-known moiré patterns. A
naive solution would be to remove the pre-filter and resort to nearestneighbor downsampling, but then moiré patterns would appear.

CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism—Color, shading, shadowing, and texture;
Keywords: image pyramids, mipmapping, anti-aliasing, multiresolution filtering, display-aware filtering, smoothed local histogram filtering, bilateral filtering, local Laplacian filtering
Links:DL

PDF

Introduction

This paper introduces sparse pdf maps (sPDF-maps), which areprobabilistic image pyramids that sparsely encode probability density functions (pdfs) of local pixel neighborhoods. Our main goals
are (1) the accurate evaluation of non-linear operators at any output
resolution, and (2) scalability to gigapixel images. The sPDF-map
representation is computed in a pre-computation stage. However, at
run time different operators with arbitrary parameter settings can beevaluated interactively. sPDF-maps enable the accurate evaluationof non-linear operators directly at any resolution without accessing

range rr
computesPDF-map
coefficientspyramid levelsrrr

convolve
selectedlevelwith W

computeexpectedvaluesraspace s

b

s

s

s

c

s

d

s

Figure 2: Non-linear operator evaluation from pdfs encoded as sPDF-map coefficients. The original signal (a) is represented in coarser
pyramid levels by coefficients in the (space × range) domain (b). A non-linear operator, such as bilateral filtering, is evaluated via (c) spatialconvolution of the coefficients, and (d) computing the expectation of each pixel (blue) as the sum of look-ups in pre-computed 1D functions.

higher or lower resolutions. This property makes our approach fully
scalable to gigapixel images. Efficient out-of-core processing andrendering are greatly facilitated by only having to load the visible image tiles of the desired output resolution. In this context, a drawbackof other multi-scale representations such as wavelets, e.g., [Fattal
2009], is that they always require accessing multiple scales.
Each pyramid level of an sPDF-map comprises a set of coefficientsthat sparsely represent the 2D image domain combined with its 1Drange (i.e., intensities). Such a combined 3D domain was previously
exploited for fast bilateral filtering [Paris and Durand 2006; Chenet al. 2007]. Non-linear operators in the 2D image domain can thenbe evaluated from linear convolutions in the 3D sPDF-map domain.
We compute the sPDF-map coefficients via Matching Pursuit [Mallat
and Zhang 1993], which is common in sparse signal representation.
In fact, our representation is so sparse that we usually require no
more coefficients than image pixels. Despite this compactness, ourrepresentation implicitly represents the entire continuous 1D pdf
associated with each image pixel. See Fig. 2 for an overview.

preserving filters, e.g., for shape and detail enhancement [Fattal
et al. 2007]. A very powerful family of multi-scale representations are wavelets [Mallat 2009], which can be adapted to better
observe edges, e.g., edge-avoiding wavelets [Fattal 2009], or à-trous
wavelets [Hanika et al. 2011]. Our goal, however, is the direct evaluation of non-linear operators at any desired output resolution, whichenables their scalable evaluation for gigapixel images [Kopf et al.
2007], whose display, filtering, and processing is an active area ofcurrent research [Kazhdan and Hoppe 2008; Summa et al. 2010].
Anti-aliasing. Our main goal is preserving the non-linearity of image operations in a pyramid of images that have been pre-filtered foranti-aliasing [Crow 1977]. Our approach is conceptually equivalent
to knowing subpixel coverage and weighting subpixel contributions
accordingly [Carpenter 1983]. For shadow mapping, this is equivalent to percentage-closer filtering [Reeves et al. 1987]. We represent
the entire range of pixel values, although we “forget” the exact pixel
locations like locally orderless images [Koenderink and Van Doorn1999]. In contrast, many methods represent only Gaussian distributions faithfully [Donnelly and Lauritzen 2006; Younesy et al. 2006].

We present a general method that uses the sPDF-map coefficientsfor the efficient and accurate evaluation of a variety of non-linear
operators in the combined (space × range) domain, with the range
being continuous (Fig. 2). We introduce the sPDF-map data structure, which compactly stores the coefficients in a manner similar tostandard image pyramids. This facilitates efficient parallel implementation on GPUs. We demonstrate the versatility of sPDF-maps
for color mapping, edge-aware processing with local Laplacian filters [Paris et al. 2011], smoothed local histogram filters [Kass andSolomon 2010], and bilateral filters [Chen et al. 2007]. Our method
unifies all of these operations, enabling their computation with arbitrary parameters from the same pre-computed data structure.

Non-linearity of shading. An important non-linear operation ingraphics is shading. Computing pre-filtered normal maps is a hard
problem. Approaches in this area that are similar to ours are fitting
the underlying distribution of normals [Han et al. 2007], or the
distribution of reflectance (BRDFs) [Tan et al. 2008]. Both of theseapproaches use expectation maximization (EM) for fitting, whereas
we obtain a sparse representation via Matching Pursuit [Mallat andZhang 1993]. A major difference of sPDF-maps is that they exploit
coherence between neighboring pixels in the 3D (space × range)
domain, which leads to very few coefficients, whereas all previous
approaches are fitting the distribution of each pixel individually.

Main contribution. Our work is the first to consistently and compactly store a pre-computed representation of the continuous 1D
probability distribution describing the neighborhood of each image
pixel in each level of an image pyramid, while providing efficient andeasy means for directly using this pdf representation for non-linear
operator evaluation anywhere in this multi-resolution hierarchy.

Local Laplacian filters enable powerful edge-preserving filtering
via simple non-linear operators [Paris et al. 2011]. sPDF-maps allow
these filters to be evaluated for gigapixel images at near-interactive
rates, by computing each Laplacian pyramid coefficient directly.
This results in O(n) complexity, instead of the original O(n log(n)).
Recent work by Aubry et al. [2011] also presents an O(n) approach,however at the cost of storing multiple pre-computed pyramids.

2

Bilateral filtering [Tomasi and Manduchi 1998] is a non-linear
weighting of the neighborhood N (p)1 of a pixel p, with intensity
Ip , with Gaussians Gσs (spatial weight) and Gσr (range weight):

Related Work

Mipmaps [Williams 1983] are pre-computed texture pyramidsthat enable efficient pre-filtering for anti-aliasing in texture mapping [Heckbert 1989]. In contrast to using multi-scale decompositions for image processing, mipmaps focus on using any desired
resolution directly. This is also an important property of sPDF-maps,
but the encoding of probability distributions enables us to faithfully
apply non-linear operators after the pre-computed pre-filtering.
Image pyramids such as Gaussian and Laplacian pyramids [Burtand Adelson 1983] target multi-scale analysis and filtering. Multiscale image decompositions can also be computed via edge-

Ipbilat. =

1
wp

X

Gσs (kp − qk) Gσr (|Ip − Iq |) Iq ,(1)q∈N (p)normalizing with the sum of all weights wp . This equation is nonlinear, because the intensity Iq appears in the argument of Gσr (·).
It can be computed from linear convolutions by going from the
1 Technically, the sum is over the whole image, but we use a neighborhood

N (p) in 2D (N3 (p, Ip ) in 3D) to be closer to the actual implementation.

~I

H0 =~I ∗ K

range filter K

dPDF map H jsparse V jcoefficients

sPDF map
data structure

spatialpre-filter W
pyramid levels

spatialconvolution
expectation computation
histogramslicing

non-linear
filtering

input image
1~
Map Input Image to I*K Volume2

Compute Dense PDF Map

3

Compute Sparse PDF Map

4Image Reconstruction

Figure 3: Overview of the sPDF-map pipeline. The input image is lifted from the 2D spatial domain into the 3D (space × range) domain. We
then apply range smoothing via the filter K. Next, we compute an intermediate dPDF-map Hj , by iteratively applying the spatial pre-filter W
and downsampling. The dPDF-map is then transformed into an sPDF-map using Matching Pursuit. Each level of an sPDF-map comprises
sPDF-map coefficients stored as a pair of index and coefficient images. At run time, non-linear operators are evaluated using these coefficients.

2D spatial domain to the 3D (space × range) domain, and using
homogeneous notation for normalization [Paris and Durand 2006]:
hi X
wp Ipbilat.
, wp =G(σs ,σr ) (p, Ip ) − (q, Ib ) I˜ (q, Ib ), (2)(q,Ib )∈N3 (p,Ip )with quantized intensity bins b (1 ≤ b ≤ B) with intensities Ib ,˜ Ib ) = [Iq , 1], if Iq maps to bin b, and [0, 0] otherwise.
and I(q,G(σs ,σr ) = (Gσs ⊗ Gσs ) ⊗ Gσr , where ⊗ is the tensor product.
sPDF-maps also operate in the (space × range) domain. However,they represent a continuous range axis, which we denote as r (insteadof Ib ), whereas a bilateral grid [Chen et al. 2007] quantizes the range
with B bins. In fact, each pyramid level of an sPDF-map can be
dynamically sampled into a bilateral grid of any range quantization.
Also, the bilateral grid leverages coarse-resolution computations forhigher performance, whereas we are interested directly in coarser
output images, like display-aware approaches [Jeong et al. 2011].
Histograms are often employed as a compact, quantized representation of distributions. For example, Thompson et al. [2011] use pixel
histograms (hixels) in the context of representing data uncertainty.
Smoothed local histograms estimate the local distribution of pixel
values, which enables a variety of non-linear filters [Kass andSolomon 2010]. Channel smoothing [Felsberg et al. 2006] is a similar robust filter. We compute normalized smoothed local histogramsfor the initial estimation of sampled pdfs, but our sPDF-maps representation then only stores sparse coefficients that represent continuous pdfs. The smoothed local histogram of a pixel p is defined as:
hp (Ib ) =X

W (p − q) K (Ib − Iq ),(3)q∈N (p)where the kernel K sums to one and smoothes the histogram, andthe kernel W performs smoothing in the spatial domain [Kass andSolomon 2010]. We will exploit that Eq. 3 can also be considered in
the continuous range r, instead of quantized histogram bins b.

Basic setup. We consider a given pixel p with neighborhood N (p).
We take a probabilistic viewpoint, and consider the pixel value at p
to be a continuous random variable Xp . We describe the distributionof Xp by its probability density function pdfp (r), informed by the
pixel values in the neighborhood N (p). We then have two basic
goals: First, we want to come up with a compact representation forthe pdfp (r) of every p that, despite its compactness, allows for easy
access. Second, given this representation, we want to use it directly
for the evaluation of a variety of non-linear image operations.
Non-linear image operations. Given pdfp (r), we will show that
we can compute the result of many non-linear operations as the
expectation of a suitably chosen function tp applied to the randomvariable Xp , using an appropriate normalization factor wp :E [tp (Xp )] =Method Overview

Basic intuition. Our goal is to apply non-linear image operators atcoarse image levels without sacrificing accuracy. Similarly to thelinearity of Eq. 2 vs. the non-linearity of Eq. 1, this is possible bypre-filtering and downsampling the 3D (space × range) domain of animage, instead of its 2D spatial domain. We therefore substitute eachpixel p by a 1D function pdfp (r) that describes the entire range r.

Z 1
tp (r)pdfp (r) dr.

(4)0With tp (r) = r and wp = 1, Eq. 4 simply computes E [Xp ], theexpected value of the random variable Xp . In order to gain intuition,
we consider evaluating a bilateral filter. We can estimate pdfp (r) as
the normalized smoothed local histogram hp (r) (Eq. 3), with thespatial filter W = (Gσs ⊗ Gσs ) (see Eq. 1), and K as the Diracdelta δ0 (r). We can then evaluate a bilateral filter using Eq. 4 with:Z 1
tp (r) = r · Gσr (|Ip − r|) , and wp =0

tp (r)
pdfp (r) dr,r(5)

using the range weight Gσr from Eq. 1. In this way, the spatialsmoothing via Gσs is contained in pdfp (r), while the range weighting of the bilateral filter is contained in the particular choice of tp (r).
Smoothed histogram volumes. In order to prepare for carrying out
operations in 3D instead of in 2D, we define a smoothed histogram
volume H(p, r) as:
H(p, r) =

3

1
wp

1
hp (r), with wp =wp

Z 1
hp (r) dr,(6)0with hp (r) defined by Eq. 3, in the continuous range r, and wp used
to guarantee proper normalization.
˜ Ib )
sPDF-map coefficient volumes. In analogy to the volume I(p,in Eq. 2, we define a sparse volume V (p, r) that is non-zero only at
specific positions (pn , rn ), i.e., V (pn , rn ) = cn with cn 6= 0. We
then represent V solely by the set of tuples (pn , rn , cn ), which wecall sPDF-map coefficients (see Fig. 2). We compute a V (p, r) thatis as sparse as possible, while approximating H(p, r) well via:
H (p, r) ≈X(W ⊗ K) (p, r) − (q n , rn ) V (q n , rn )(q n ,rn ,cn )q n ∈N (p)= V ∗ (W ⊗ K),(7)

In the next step, an intermediate dense pdf map is computed by
iteratively applying a spatial pre-filter and downsampling (Sec. 4).
The pre-filter size is determined to match the desired footprint ofeach pixel (Fig. 4). These two steps together result in smoothed
histograms that estimate the pdfp (r) that we require. The dense pdf
map is then converted into an sPDF-map level by level, by iteratively
computing sPDF-map coefficients (pn , rn , cn ) (Sec. 5). After thesPDF-map has been computed, the dense pdf map is discarded.

where ∗ denotes a convolution, and (W ⊗ K) is a 3D kernel.
A crucial property of our approach is that coefficients influence aneighborhood of pixels N (p) as determined by the spatial kernel W,
instead of individual pixels. This reduces the number of required
coefficients by exploiting coherence in the (space × range) domain.
Sec. 5 describes how we compute V (p, r). For now, we assume that
we know the set of sPDF-map coefficients (pn , rn , cn ).
Non-linear operations with sPDF-map coefficients. Instead ofcomputing the pdfp (r) = H(p, r) explicitly from V (p, r), we plugthe convolution from Eq. 7 directly into Eq. 4. Also exploiting theseparability of (W ⊗ K) then allows us to rewrite Eq. 4 as:
Z 1
1t̃p (r) (V ∗ W ) dr, with t̃p = tp ∗ K.
E [tp (Xp )] =wp 0
(8)This formulation has moved the range convolution with K into
the new function t̃p . Since V (p, r) only consists of sparse coefficients (pn , rn , cn ), Eq. 8 can be computed exactly as the sum:E [tp (Xp )] =

1
wp

X

t̃p (rn )W (p − q n ) V (q n , rn ).

In order to evaluate a non-linear operator, we therefore first define asuitable function tp and convolve it with the range smoothing kernel K to obtain t̃p . Then, we evaluate Eq. 9, which is a simple sum
over all sPDF-map coefficients, multiplying each coefficient withthe spatial smoothing weight W, times the look-up result from t̃p .
Image pyramids store a pre-computed hierarchy of images
smoothed and downsampled in the spatial domain. sPDF-maps
are similar pyramids, but additionally smooth and represent therange r of the image. For each pixel p, the function pdfp (r) represents the pixel values in the neighborhood of p corresponding to its
footprint [Greene and Heckbert 1986] in the original image (Fig. 4).
sPDF-maps pipeline overview. Our pipeline for computing andusing sPDF-maps is illustrated in Fig. 3. The smoothing in the range
and in the spatial domain, respectively, are computed in two different
stages. The initial step computes a smoothed histogram for eachpixel in the original image, with a spatial neighborhood of 1 × 1
pixel and without spatial smoothing. The required quantization rateis determined by the amount of range smoothing applied.

convolution

Dense PDF Maps

We call a pyramid of smoothed histogram volumes Hj (p, r), consisting of pyramid levels j, a dense pdf map (dPDF-map). As in
most image pyramids, from level j to j + 1 we decrease the spatial
resolution by a factor of two in each axis. However, the range axis r
is always sampled with B samples, such that every Hj comprises B
image layers. Each layer b corresponds to histogram bin b for all
pixels in level j. The choice of B does not have to correspond to
the quantization of the intensity axis of the original image. Since
we perform range smoothing, a sampling rate that guarantees proper
signal representation can be chosen depending on K (Eq. 3).
The dPDF-map is the first step to building an sPDF-map. It is adense representation of the normalized local histograms. Later, we
compute a sparse approximation of it, which is the sPDF-map that
we actually use at run time.
We note that each Hj is similar to the homogeneous channel of abilateral grid [Chen et al. 2007].

(9)(q n ,rn ,cn )q n ∈N (p)histogram smoothing estimated pdf

4sparse pdf
approximation

Figure 4: Pixel neighborhoods in sPDF-maps and dPDF-maps.
The pdf of each pixel represents the neighborhood (footprint) of the
pixel in the original image. Continuous pdfs are initially estimatedas discrete smoothed local histograms, but then represented by
sPDF-map coefficients in the combined (space × range) domain.

4.1

Pixel neighborhoods in dPDF-maps

The spatial filter W in a regular smoothed local histogram (Eq. 3)
corresponds to the desired spatial neighborhood size and smoothing,as discussed by Kass and Solomon [2010]. However, for computing
a dPDF-map, we define W to act as the spatial pre-filter for coarser
pyramid levels. Thus, W defines the neighborhood of each pixel p
as the footprint of p in the original image. This is illustrated in Fig. 4.
In principle, any spatial pre-filter W can be chosen for dPDF-map
computation. Due to their straightforward iterative evaluation, we
use the 5 × 5 Gaussian-like filter of Burt and Adelson [1983], or the
2 × 2 box filter commonly used in mipmapping [Williams 1983].

4.2

Iterative dPDF-map computation

In principle, every pyramid level j of a dPDF-map can be computed
directly from the original image. However, it is more efficient tocompute the pyramid iteratively, i.e., computing each level j withj > 0 from (j − 1), where j = 0 corresponds to the original
image resolution. This is analogous to the iterative computation of
mipmaps and Gaussian pyramids. However, instead of smoothing
and downsampling a regular image in each level j, the computationof a dPDF-map requires doing this for each image layer b.
Separable smoothing in space and range. Instead of computing
the spatial smoothing (W ) and the range smoothing (K) together as
in Eq. 3, we perform these two steps separately, using the approach
depicted in Fig. 3 1 and 2 . We first compute the smoothed histogram volume for level j = 0, i.e., H0 (p, r), using Eq. 3 for eachpixel p. However, we do not perform spatial smoothing in this step,i.e., we set W to a 1 × 1 box filter. We perform range smoothing
using a Gaussian kernel, i.e., K = Gσr . In the next step, we compute all downsampled levels j with j > 0 iteratively. Each Hj (p, r)
is computed from Hj−1 (p, r), by applying the spatial smoothingkernel W to each of the B image layers, followed by downsamplingby a factor of two in each spatial axis.

˜ r) withFormally, we lift the input image I(p) into a volume I(p,˜ r) = 1 when r = Ip , and zero otherwise, and compute Hj as:
I(p,
H0
Hj

=
=

I˜ ∗ K,
2↓ Hj−1 ∗ W for all j > 0,

(10)
(11)

where 2↓ performs spatial downsampling by a factor of two.
a5

Storing gigapixel images as dPDF-maps is not feasible, because
each Hj consists of B image layers. Therefore, our goal is tocompute a sparse pdf map (sPDF-map) from the dPDF-map, anddiscard the dPDF-map afterward. The sPDF-map also consists ofpyramid levels j, comprising the sparse volumes Vj (p, r).
sPDF-map pyramid levels and coefficients. In contrast to a dPDFmap, each pyramid level of an sPDF-map solely consists of thesparse set of sPDF-map coefficients (pn , rn , cn ) introduced inSec. 3. The set of coefficients of level j can be used directly with
Eq. 9 in order to evaluate non-linear image operations in level j
without accessing other pyramid levels. We compute each level j ofthe desired sPDF-map directly from its corresponding level j in the
input dPDF-map. Each level is completely independent of all otherpyramid levels, which facilitates out-of-core approaches.
sPDF-map computation. For each pyramid level j, our goal is tocompute the sparse volume Vj (p, r) represented by the smallest
set of sPDF-map coefficients (pn , rn , cn ) that well approximate apre-scribed dPDF-map level Hj (p, r) in the L2 sense, when Hj iscomputed using Eq. 7. We thus quantify similarity using the L2
distance between each dPDF-map level Hj , and its approximation
based on the corresponding coefficients (pn , rn , cn ):E(Vj ) = kHj − Vj ∗ (Wj ⊗ K)k,(12)

where (Wj ⊗ K) acts like (W ⊗ K) in Eq. 7. However, in thiscontext we will call the spatial kernel Wj the reconstruction filter.
(Wj ⊗ K) is the basis function used for fitting, which is then usedlater for image reconstruction. Note that Wj does not have to match
the pre-filter W used in the dPDF-map computation (Eq. 11). In
principle, Wj can also be chosen independently for each sPDF-map
level j, but we currently use the same Gaussian kernel for all levels.

5.1

b

c

Sparse PDF Maps

Computation of sPDF-map coefficients

In order to compute the sparse volume Vj (p, r) of each sPDFmap level j, we use the Matching Pursuit algorithm of Mallat andAlgorithm 1 Data fitting via Matching Pursuit
Input: dense histogram volume Hj
Output: sparse volume Vj given by a list of triplets (pn , rn , cn )1: E = kHj k;
2: n = 0; chunk = 1;
3: while E > ε and chunk ≤ max chunks do4:
Vj ← ∅
5:
while n < chunk · mj do
6:find (pn , rn ) that maximizes the inner product of Eq. 13
cn = hHj (s, r), (Wj ⊗ K) ((s, r) − (pn , rn ))i7:
Hj ← Hj − cn (Wj ⊗ K) (pn , rn );8:
E ← kHj k;
9:
Vj ← Vj ∪ (pn , rn , cn )10:n ← n + 1;
11:
end while12:
Vj (chunk) ← Vj ; chunk ← chunk + 1;
13: end while

Figure 5: sPDF-maps vs. image pyramids. An image pyramid (a)consists of a single image per level. Each level of an sPDF-map(b) comprises pairs of index and coefficient images that store thesPDF-map coefficients describing the (space × range) domain (c).

Zhang [1993]. Matching Pursuit is a greedy iterative algorithm that,in our framework (Algorithm 1), finds a position (pn , rn ) at eachstep of the iteration, such that the inner product of the basis function(Wj ⊗ K) centered at (pn , rn ) with Hj is maximized:(pn , rn ) = argmax hHj (s, r), (Wj ⊗ K) ((s, r) − (q, i))i.
(q,i)∈D(Hj )(13)
Finding each (pn , rn ) requires an exhaustive search over all possible(q, i) in the domain of Hj (p, r). The corresponding cn is thencomputed as given in Algorithm 1. In practice, we only try a discrete
set of range positions i for each q, and sample the L2 errors ateven fewer range positions r. Denoting the input quantization by B
(e.g., B = 256 for 8-bit images), our implementation considers B
different i, and computes the L2 errors at B/16 different r.

5.2 sPDF-maps data structure
In order to store the Vj (p, r) efficiently and facilitate GPU implementation, we define the data structure illustrated in Fig. 5b, mimicking the basic geometry of a standard image pyramid (Fig. 5a). In
contrast to standard image pyramids, each level j of an sPDF-map
consists of one or multiple coefficient chunks, each of which is stored
as a pair of images: an index image, and a coefficient image. Fig. 5b
depicts a single coefficient chunk, i.e., one such image pair per level.
Coefficient chunks. We define the concept of a coefficient chunk inorder to be able to store all data comprising an sPDF-map in images
and facilitate parallel image reconstruction on GPUs. Each coefficient chunk is simply a sequence of mj coefficients (pn , rn , cn ),where mj is the number of image pixels in level j. Therefore, instead of computing a single sequence of coefficients, Algorithm 1
computes coefficients in multiples of coefficient chunks. When onlya single chunk is computed, this implies that there is exactly one
coefficient per image pixel on average. That is, some pixels might
not have any coefficients, whereas other pixels might have multiple
coefficients, but the total number of coefficients equals the total number of pixels. When two coefficient chunks are computed instead,
there will be two coefficients per pixel on average, and so on.
In principle, coefficients could be represented as a stream of coefficients, but our choice of grouping coefficients into chunks enables
storing all coefficients in images. This guarantees a constant memory
overhead per pixel and facilitates using 2D textures on GPUs.
Index and coefficient images. The sequence of coefficients computed by Algorithm 1 is sorted in order of descending inner product. However, this order does not facilitate efficient parallel reconstruction. In each chunk, we therefore reorder the coefficients bygathering all coefficients of each pixel p into consecutive memory
locations. We then tightly pack (in scanline order) all coefficients

Figure 6: Fitting quality of sPDF-maps. (a) Flower image (512 × 512). We evaluate sPDF-map level 1 (256 × 256). We show the distributionof coefficients for (b) the first coefficient chunk, and (c) the second chunk (64K coefficients each). White maps to 9 coefficients. We compare
E[Xp ] computed from the sPDF-map with the one from the dPDF-map. (d) shows the scaled difference between these two images with one (lefthalf) and two (right half) coefficient chunks. For one chunk, the RMSE between the two images is 2.39, PSNR 40.57. For two chunks the RMSE
is 2.22, PSNR 41.20. (e) the green curve shows how the initial RMSE of 93.75 between the two expected value images decreases as coefficientsare added, the blue curve shows how the initial RMSE of 0.0113 between the dPDF-map and its approximation by the sPDF-map decreases.
These errors are for the downsampled luminance channel of (a) using a 5 × 5 Gaussian for downsampling (W ) and reconstruction (Wj ).

of all pixels into the coefficient image C(·) of this chunk. In orderto remember where the coefficients of each pixel p are stored, we
also construct an index image X(p) that stores a pair of values(index, count)p for each p: the position of the first coefficient of pin the coefficient image: indexp , and the number of coefficients of p:countp . Each entry in C(·) then simply consists of the pair (rn , cn ),because pn from the tuple (pn , rn , cn ) is implicit in X(p).
The coefficients (rn , cn ) of pixel p are therefore enumerated asfollows: C (X(p).index) .. C (X(p).index + X(p).count − 1).

6Image Reconstruction

Given the set of sPDF-map coefficients (pn , rn , cn ), which comprise the sparse volume Vj (p, r), we can reconstruct the output
image corresponding to various non-linear operations using Eq. 9.
The non-linear operator is determined by the choice of function
tp (r) and the normalization wp . However, exploiting the propertiesof different classes of operators, and the corresponding definitionsof tp (r), enables different strategies and conceptual optimizations.

6.1

Global functions

The simplest non-linear operator applies the same function t to every
pixel p, i.e., tp (r) = t(r) for all p. A global t is sufficient to apply
arbitrary global re-mapping functions, e.g., color maps. In this case,the spatial convolution in Eq. 9 does not need to be computed forevery coefficient individually. We can rewrite Eq. 9 as:
E [t(Xp )] =

1
wp

X

Wj (p − q)q∈N (p)Xt̃(rn )V (q n , rn ). (14)(q n ,rn ,cn )q=q nRemember that V (q n , rn ) = cn . Since t is the same everywhere,
the second sum can be computed once per pixel p, and then used inall spatial convolutions. This enables evaluating Eq. 14 as follows(setting wp = 1). We first compute two images T (p) and N (p) as:XX
T (p) =t̃(rn )cn , and N (p) =K̄(rn )cn ,(15)(pn ,rn ,cn )
p=pn

R1

(pn ,rn ,cn )
p=pn

where K̄(r) = 0 K(x − r) dx is the integral of K centered at r
and clamped to the range [0, 1]. We compute N (p) in order to enableproper normalization, since the fitting operation of Algorithm 1 does
not guarantee pdfs that sum exactly to one. Then, both images areindividually convolved with Wj , and the output image is obtainedby dividing each pixel in T by the corresponding weight in N :T ∗ Wj (p).
(16)
E [t(Xp )] =N ∗ Wj (p)Color mapping. We can compute the result of anti-aliased color
mapping by simply using the color map as t(r) (handling eachchannel individually). Both t̃(r) = t(r) ∗ K and K̄(r) can be precomputed, because they are the same for all pixels. Computationally,
we only have to perform a look-up in t̃(r) and K̄(r) per coefficient,sum the weighted look-ups for each pixel, perform a 2D convolution
for T and N each, followed by one division per pixel. Despite its
simplicity, this approach accurately evaluates Eq. 8 for any global t.

6.2

Local functions

If the function tp (r) is different for every pixel p, the simplificationof Eq. 9 described in the previous section cannot be used directly,
because the spatial convolution must not mix different tp .
Local Laplacian filters can be computed via Eq. 9 by defining thefunction tp (r) to be the re-mapping function for each pixel p (Pariset al. [2011], Sec. 5.2, Fig. 6). These functions perform local contrast
adjustment around each pixel p, depending on its intensity Ip , which
we compute as E[Xp ] (Sec. 3). Apart from the inaccuracies causedby Algorithm 1, no normalization is necessary, i.e., conceptually
wp = 1 (Eq. 8). However, in practice a normalization factor wp 6= 1
must be computed and divided by, with K̄(r) as defined above:X
wp =K̄(rn )Wj (p − q n ) cn .
(17)
(q n ,rn ,cn )q n ∈N (p)We note that the computation of each Laplacian coefficient depends
on an upsampling step of Gaussian pyramid coefficients, because
it is obtained as the difference between Ip from level j and thecorresponding upsampled intensity from (j + 1). It is important thatfor each upsampled pixel p in pyramid level j, the corresponding
neighborhood in level (j + 1) is re-mapped with the same tp (r).
Bilateral filters can be computed via Eq. 9 and the definitions oftp (r) and wp given in Sec. 3. However, the spatial kernel W must
include both the reconstruction kernel Wj and the spatial neighborhood of the bilateral filter (Wbl ). Both can either be pre-convolved(W = Wbl ∗ Wj ) or be evaluated in two successive steps. Alternatively, instead of evaluating the bilateral filter directly, it is possible

ImageNight Scene
BelliniNASA Bathymetry
Machu Picchu
Rock
Beach
FlowerBarbaraSalt and Pepper Noise

Resolution
47, 908 × 7, 531
16, 898 × 14, 824
21, 601 × 10, 801
9, 984 × 3, 328
876 × 584
800 × 533
800 × 533
512 × 512
320 × 428

Size [Mpix.]
361
250
233
33
0.51
0.43
0.43
0.26
0.14

Pyramid Levels
17
16
16
15
11
11
11
10
10

Overhead [Mpix.]
120
83
77
11
0.17
0.14
0.14
0.09
0.05

Pre-Computation [s/Mpix.]
127
133
128
132
139
135
142
138
143

Table 1: Images in this paper with corresponding sPDF-map properties and pre-computation times. The overhead column gives the total
number of pyramid pixels (in Mpixels) in addition to the original image size. This is the same number of additional pixels required by a regularmipmap or Gaussian pyramid. For each of these additional pyramid pixels, every sPDF-map coefficient chunk requires 64 bits of storage per
pixel in our current implementation. All our result images use a single coefficient chunk, except where indicated in Fig. 12 and Fig. 13e-m.

to dynamically compute a bilateral grid [Chen et al. 2007] from an
sPDF-map, similar to the histogram slicing process described below.
This enables all performance optimizations of the bilateral grid.

6.3

Histogram slicing

Smoothed local histogram filters [Kass and Solomon 2010] arecomputed from histograms hp (Ib ) as defined by Eq. 3. sPDF-maps
support this kind of filter by simply sampling the continuous pdfp (r)from the sPDF-map into the hp (Ib ) for each pixel p. This can bedone with an arbitrary sampling rate B without changing the sPDFmap. We call this process slicing the sPDF-map into histogram bin
images Hb (p), each of which contains a bin b, i.e., hp (Ib ), for all p.
Each bin image can be computed identically for all p by evaluating
Eq. 8 using the Dirac-delta tp (r) = t(r) = δIb (r) and wp = 1.
Because the required function t(r) is the same for all p, each bin image can be computed as follows: First, we compute an intermediateslice coefficient image Sb (p) corresponding to bin b, as:XXSb (p) =δIb ∗ K (rn )V (pn , rn ) =K(Ib − rn )cn ,(pn ,rn ,cn )
p=pn(pn ,rn ,cn )
p=pn

(18)
which sums up the contributions of all coefficients of pixel p,weighted by a look-up into the range kernel K corresponding to Ib .
Each Sb (p) must then be convolved by W :Sb∗ (p) = Sb ∗ W (p),(19)
where the desired spatial histogram neighborhood (Wsh ) can either
be computed at the same time (W = Wsh ∗ Wj ), or in two successive convolution passes. In order to obtain a normalized smoothed
histogram that represents a pdf, each intermediate bin image mustthen be normalized by the sum of all bin images. This yields the
final histogram bin image Hb (p) of each bin b as:
B

Hb (p) =X ∗Sb∗ (p)with N (p) =Sb (p).
N (p)

(20)

7.1

Pre-computation

dPDF-map level j = 0. The first step is reading the input image,conceptually computing the corresponding volume I˜ (Sec. 4.2), andapplying the range kernel K to compute sPDF-map level 0, i.e., H0 .
In order to make this scalable, we read the input image one scanline
at a time. For each pixel p, we directly compute its contribution
to H0 by applying the range filter K. This results in a plane of
H0 in (space × range) for each scanline. Instead of writing eachplane to disk immediately, we accelerate disk access by avoiding
interleaving reads from the input with writes to the output. We
forward each plane to a simple caching system that only writes to
disk when a user-defined memory limit is exceeded. This guarantees
zero overhead for images that can easily be processed in-core.
Reformatting into tiles. We want to store each dPDF-map level asa collection of tiles, with a tile size of (n + o) × (n + o) pixels,
where n denotes the inner area of a tile, and o denotes an overlap
region that is replicated amongst neighboring tiles. In order to do
this, we collect the corresponding planes of H0 computed above,and for each tile re-arrange them into a stack of B slice images thatsample the range. In our implementation, we are using an innertile size of n = 256, and an overlap of o = 8, which facilitatesstraightforward upsampling in Gaussian pyramids without accessing
neighboring tiles. The overlap depends on the size of the spatial
filters that must be evaluated. Our relatively small tile size ensuresthat the tiles fit into the L1 cache of the processor and can also be
used directly by the tiled rendering system for image reconstruction.
dPDF-map levels j with j > 0. To compute one new level in thedPDF-map, 2 × 2 tiles are loaded at a time, and downsampled toa single tile by applying the kernel W (Sec. 4.1). This is repeateduntil all tiles comprising level j have been processed. Instead ofwriting each downsampled tile directly to disk, it is passed to ourcache which is flushed to disk once it exceeds a pre-defined size. In
this way, the full pyramid is computed incrementally by applying
the same operation to the downsampled tiles until only a single tileis left. Finally, this tile is iteratively downsampled to a single entry.

b=1

Using different look-up tables for the range kernel K, including derivative and integration kernels as used by Kass andSolomon [2010], using the appropriate normalization, enables all
smoothed local histogram filters to be evaluated via an sPDF-map.

7

Implementation and Performance

In order to support the computation of sPDF-maps of gigapixel
images using a limited memory footprint, we have implemented theentire pipeline depicted in Fig. 3 in an out-of-core fashion.

sPDF-map computation. Next, we compute the sPDF-map fromthe dPDF-map tile by tile. Because we process each tile independently, our approach naturally scales to images of arbitrary resolution. In every iteration of Algorithm 1, the best position (pn , rn )with its corresponding coefficient cn must be computed by finding
the (pn , rn , cn ) that reduces the residue E the most. A straightforward but inefficient way of doing this is just going through all
possible (pn , rn ), computing the coefficients and residue decreases,
then choosing the position and coefficient with the highest decrease.
Another alternative is to keep a priority queue of the residue decreases for all possible positions. However, this queue can quickly

Figure 7: Color mapping. (Top row) Gray-scale 21, 601 × 10, 801 (233 MPixels) bathymetry image from the NASA Blue Marble collection [NASA 2005]. (Center row) Anti-aliased color mapping computed from the sPDF-map; (Bottom row) Standard pre-filtering anddownsampling followed by color mapping: coarser resolutions introduce wrong colors, and whole structures are changing or disappearing.
become too large to be practical and can be slow to maintain. Instead, we subdivide each tile of the dPDF-map into smaller cubicalblocks in (space × range) and track the optimal positions and residue
decreases for each block in a small list. For each iteration, we thenperform a parallel reduce step on the list to determine the maximum
residue decrease and the corresponding position. We maintain thislist continuously, and update all affected entries whenever the residue
is modified. To improve efficiency further, we have implemented
this approach on the GPU using CUDA.

performance. In our implementation, each pixel has a constantmemory consumption of 64 bits per coefficient chunk. We pack the(index, count)p pair of each pixel in the index image X(p) into asingle 32-bit integer (24 bits for indexp , 8 bits for countp ), and use
an integer texture for X(p). The (rn , cn ) pair of each pixel in thecoefficient image C(·) is stored as two 16-bit floats in a 16-bit float
texture. In addition, our implementation incurs the typical overhead
for the tiling data structure that enables us to properly deal withgigapixel images. However, this is independent of the sPDF-maps
representation and inherent to the use of a tiled out-of-core approach.

7.2

The image reconstruction stage of E[t(Xp )] for a global function t
is extremely efficient. On a GTX 580 GPU, we compute anti-aliased
color mapping (Sec. 6.1) in 12 ms per Mpixel and coefficient chunk,with a Wj of 5 × 5. A histogram slice (Sec. 6.3) can be obtainedat the same rate. This performance is roughly constant for all ourimages, and scales linearly with the number of coefficient chunks.
For Laplacian filtering, the performance of computing Laplacianpyramid coefficients is independent of the pyramid level that iscomputed. This is not the case for the original approach of Pariset al. [2011], where the time for computing each coefficient must
increase with the original image resolution. The approach by Aubryet al. [2011] is much faster, but depends on pre-computing fullresolution pyramids, which makes its scalability to gigapixel images
unclear. Using an sPDF-map, all Laplacian pyramid coefficients fora 1 MPixel image can be computed in around 320 ms, independent of
the original image resolution, scaling linearly. Overall performancemust also take into account the Laplacian pyramid collapse, which
depends on the number of pyramid levels. For a 1024 × 1024 view
showing a zoom-in of level 0 of the 47, 908 × 7, 531 Night Scene
image (Fig. 8), the pyramid collapse takes around 5 ms, and overall
computation time for the entire 1024 × 1024 view is below 400 ms.

Image reconstruction

Our implementation computes the output image E[tp (Xp )], given
various tp (r), from the sPDF-map representation on the GPU using
CUDA. We create a thread for every output pixel p, and computethe necessary intermediate images, as well as the final image, by
evaluating the equations given in Sec. 6. Each thread accesses all
coefficients that are required to compute E[tp (Xp )] for its output
pixel p by first accessing the corresponding pixel in the index image
X(p), and then enumerating and fetching all pairs (rn , cn ) from thecoefficient image C(·) (Sec. 5.2). The required functions t̃p (r) =tp (r) ∗ K are pre-computed and stored in 1D arrays.

7.3

Performance and memory consumption

We have found that in general our fitting procedure (Algorithm 1)
converges quite quickly (see Fig. 6e). It is usually sufficient tocompute just one chunk of sPDF-map coefficients. Using morechunks does not dramatically improve the fitting quality for mostimages. Exceptions are images such as Fig. 12a and Fig. 13e.
The pre-computation time for sPDF-maps, as well as memory consumption, are summarized in Table 1 for various images. On an
NVIDIA GTX 580 GPU, we currently require about two minutes tocompute an sPDF-map with one coefficient chunk for one Mpixel.
This time scales linearly both with image size and number of chunks.
While we consider this pre-processing time feasible for large images, we believe using an optimized fitting algorithm could improve8

Applications of sPDF-maps

We view the versatility of the sPDF-maps representation as one ofits biggest strengths. In this section, we illustrate several non-linear
image operators that are all computed from the same data structure.

Figure 8: Local Laplacian filtering with an sPDF-map in O(n) time. (Top row) Night scene of resolution 47, 908 × 7, 531 (361 Mpixels).
The top third of the image is shown with detail enhancement (σr = 0.2, α = 0.25), the center third is the original image, and the bottom thirdis shown with smoothing (σr = 0.2, α = 3.0). (Bottom row) Images used by Paris et al. [2011]: left-hand image of each pair with detailenhancement (σr = 0.4, α = 0.25), right-hand image with smoothing (σr = 0.2, α = 2.0). RGB color channels were computed separately.

8.1

Anti-aliased color mapping

Fig. 7 illustrates anti-aliased color mapping vs. first downsampling
and then color mapping for a high-resolution gray-scale image fromthe NASA Blue Marble collection [NASA 2005]. The color map is
a standard false color coding often used by scientists to visualize theunderwater depth in bathymetry images. The bottom row illustrates
the two common problems of color-mapping gigapixel images: (1) Incoarser resolution levels, wrong color values are introduced, because
the color map is applied after downsampling. This is especiallyvisible around the islands of the Philippines depicted in the left zoomin. (2) Structures change their shape and/or topology, disappear, orappear from pyramid level to pyramid level. In the left zoom-in,
entire islands disappear, whereas in the right zoom-in structures ofvery large depth in the Mariana Trench successively disappear.
Both of these problems can be avoided by using the same color
map as the function t(r) with the corresponding sPDF-map. Thisis illustrated in the center row of the figure. The output image iscomputed as described in Sec. 6.1. Each channel of the RGB output
is computed separately, but from the same sPDF-map input, using
the corresponding channel of the color map as the function t(r).

non-linear filter directly for a downsampled image cannot recover
the edges already smoothed away by the linear pre-filter. In contrast,
if the non-linear filter uses the sPDF-map, it is conceptually applied
to the original image followed by pre-filtering and downsampling,which correctly preserves the non-linearity of the histogram filter.
Fig. 12 illustrates this in more detail for median filtering. The
original image (a) contains salt and pepper noise, which is easily
removed by median filtering (b). However, when the image is first
downsampled (f) and then median-filtered (g), artifacts from thedownsampling are visible. In contrast, sPDF-map median filtering
is able to remove the noise without artifacts. For this example, we
illustrate two different reconstruction kernels Wj , because matching(a) exactly is hard for sPDF-maps. Our standard Wj of 5 × 5smoothes the noise (j); 31 dB compared to (b). Using a Wj of
3 × 3 produces a better match, however at the price of using two
coefficient chunks (e, bottom half); 33 dB. If only a single chunkis used, there are not enough coefficients to match the noise (e, tophalf); 27 dB. However, median filtering quality is good for a Wj of
3 × 3 and two chunks, as well as for a Wj of 5 × 5 and one chunk.

8.4

Bilateral filtering

Fig. 8 illustrates different examples of detail enhancement andsmoothing, respectively, with the local Laplacian filtering approachof Paris et al. [2011], but evaluated using sPDF-maps. We have
implemented an interactive application for gigapixel viewing andfiltering that only computes the pixels of image tiles that are currently visible on screen. This is demonstrated in the video. The
Laplacian pyramid is only computed and collapsed for the parts
actually visible on screen. This approach is greatly facilitated by
the fact that sPDF-maps enable direct computations in each pyramid
level. We compute the filtering as described in Sec. 6.2. Fig. 9
and Fig. 10 compare our results with the publicly available Matlab
implementation of Paris et al. [2011]: Fig. 9 for a zoom-in of the
Bellini image (Fig. 1), and Fig. 10 for the flower shown in Fig. 6a.

Fig. 13 illustrates bilateral filtering using the sPDF-map representation vs. standard downsampling followed by bilateral filtering.
Applying the bilateral filter to a downsampled image loses many of
its edge-preserving qualities, because many edges have already beensmoothed away by the linear pre-filter. This problem is reduced
considerably by applying the bilateral filter to the sPDF-map as
described in Sec. 6.2. The top half of Fig. 13 illustrates a zoom-inof the Bellini image (level 4), where high-frequency detail is lostin the naive approach. The bottom half of Fig. 13 shows zoom-insof level 1 and level 2 of the Barbara image. This illustrates that
applying the bilateral filter to a downsampled image is also problematic for already anti-aliased edges, where using sPDF-maps achievesmuch better results. Note that due to the many very thin stripes inFig. 13e, we have computed the corresponding sPDF-map with a
Wj of size 3 × 3 and two coefficient chunks to obtain a good fit.

8.3

9

8.2

Fast local Laplacian filtering

Smoothed local histogram filtering

Fig. 11 illustrates examples of computing smoothed local histogramfilters in image pyramids, and compares the results of standard
downsampling and then filtering with sPDF-maps. Evaluating theDiscussion and Limitations

We envision sPDF-maps as a powerful alternative to regular image
pyramids, however enabling filter evaluation to exploit direct accessto local distribution functions with a minimal storage overhead.

Figure 9: Local Laplacian smoothing. sPDF-map results vs. the original implementation of Paris et al. (σr = 0.2, α = 2.0). (a,b) level 0(1, 744 × 1, 160); (c,d) level 3 (218 × 145). (a,c) Paris et al.; (b,d) sPDF-map. Luminance PSNR [dB] between (a,b) 35, (c,d) 36.

Approximation quality. sPDF-maps offer surprisingly good approximation quality with very few coefficients, e.g., the same numberof coefficients as image pixels, by exploiting the coherence between
neighboring pdfs in the 3D (space × range) domain. This property
also enables accurate approximation with Gaussians of constant size.
However, it would be interesting to experiment with adaptive (possibly anisotropic) kernels in the future. For most images, a Gaussianreconstruction kernel Wj of size 5 × 5 with just one coefficientchunk achieves very good quality, both visually and numerically.
The only exceptions that we have encountered are images with very
high frequency content, e.g., Figs. 12a/13e, for which a smaller Wj
of 3 × 3 and more coefficient chunks can be used. The results of a3 × 3 kernel with one chunk are numerically on a par with the naive
approach (Fig. 12), and two chunks are clearly better. Furthermore,in our experience the visual quality is often better than the PSNR
values suggest, i.e., visually closer to the ground truth than the naive
results, even for small PSNR differences. In practice, we use a 5 × 5
kernel Wj with one chunk, or a 3 × 3 kernel with two chunks.
Classes of filters. The sPDF-maps representation enables the efficient and accurate evaluation of a wide variety of filters. Essentially,sPDF-maps support any filter that can be evaluated from local 1D
pdfs, which includes many important filters used for photo editing.
However, filters that require more information cannot be evaluateddirectly. Examples would be filters that require gradients, such as
anisotropic diffusion, or filters that compare exact neighborhoods,
such as non-local means filtering. Extending sPDF-maps to higherdimensional pdfs could enable these kinds of filters in the future.

Figure 10: Local Laplacian detail enhancement. sPDF-map results vs. the original implementation of Paris et al. (σr = 0.2,
α = 0.5). Flower (800 × 533) level 0: (a) Paris et al. (b) sPDFmaps. Luminance PSNR between (a,b) 37 dB. Original in Fig. 6a.

References
AUBRY, M., PARIS , S., H ASINOFF , S. W., K AUTZ , J., AND D U RAND , F. 2011. Fast and Robust Pyramid-based Image Processing. Tech. rep., MIT. MIT-CSAIL-TR-2011-049.
B URT, P., AND A DELSON , T. 1983. The Laplacian pyramid as acompact image code. IEEE Transactions on Communications 9,
4, 532–540.
C ARPENTER , L. 1983. The A-buffer, an antialiased hidden surface
method. In Proceedings of SIGGRAPH 1983, ACM Press / ACM
SIGGRAPH, 1–11.

Apart from these limitations, we think that the biggest current limitation of our method are the pre-computation times required tocompute the sPDF-map coefficients. We would like to explore different optimization strategies for better performance in future work.

C HEN , J., PARIS , S., AND D URAND , F. 2007. Real-time edgeaware image processing with the bilateral grid. ACM Transactions
on Graphics (Proceedings of SIGGRAPH 2007) 26, 3, 103:1–
103:9.

10

C ROW, F. C. 1977. The aliasing problem in computer-generated
shaded images. In Proceedings of SIGGRAPH 1977, ACM Press/ ACM SIGGRAPH, 799–805.

Conclusions

We have introduced the new sPDF-maps representation and data
structure to compactly represent pre-computed pdfs of pixel neighborhoods in multi-resolution image pyramids. In order to use sPDFmaps in practice, we have presented an efficient unified method
for evaluating a variety of non-linear image operations from thesame pre-computed representation, which illustrates the versatility of sPDF-maps. However, we imagine that in the future more
operations can make use of the information represented by pixel
neighborhood pdfs, using the same sPDF-maps data structure.

Acknowledgements
We would like to thank the anonymous reviewers for helping us improve the paper; Peter Rautek, Thomas Theußl, Thomas Höllt, andWito Engelke for their help; Justin Solomon for a smoothed local histogram filter implementation, and Sylvain Paris et al. for their local
Laplacian filter implementation, which we used for comparisons.

D ONNELLY, W., AND L AURITZEN , A. 2006. Variance shadow
maps. In Proceedings of 2006 Symposium On Interactive 3D
Graphics and Games, ACM Press, 161–165.
FATTAL , R., AGRAWALA , M., AND RUSINKIEWICZ , S. 2007.
Multiscale shape and detail enhancement from multi-light image
collections. ACM Transactions on Graphics (Proceedings of
SIGGRAPH 2007) 26, 3, 51:1–51:9.
FATTAL , R. 2009. Edge-avoiding wavelets and their applications.
ACM Transactions on Graphics (Proceedings of SIGGRAPH
2009) 28, 3, 22:1–22:10.
F ELSBERG , M., F ORSSEN , P.-E., AND S CHARR , H. 2006. Channel smoothing: Efficient robust smoothing of low-level signal
features. IEEE Transactions on Pattern Analysis and MachineIntelligence 28, 2, 209–222.

Figure 11: Smoothed local histogram filtering. (a,b) Beach image in (a) dominant mode-filtered (luminance only) in (b). (f,g) Rock image in
HSV color model: Standard downsampling introduces strong haloes of the wrong color around the rock (f), whereas dominant mode-filtering
the H and V channels correctly preserves the circular domain of the hue channel (g). (c,d,e) and (h,i,j): Median filtering (luminance only) using
sPDF-maps vs. downsampling and then filtering. (c,h) Original zoom-ins of the Night Scene (c) and the Machu Picchu (h) images. (d,i) Medianfiltering with sPDF-maps prevents over-smoothing by properly preserving the non-linearity of the image operation. (e,j) Median filtering after
downsampling introduces strong over-smoothing that cannot be reversed by the median filter applied directly at the coarser resolution.

G REENE , N., AND H ECKBERT, P. 1986. Creating raster omnimax images from multiple perspective views using the ellipticalweighted average filter. IEEE Computer Graphics and Applications 6, 6, 21–27.
H AN , C., S UN , B., R AMAMOORTHI , R., AND G RINSPUN , E. 2007.
Frequency domain normal map filtering. ACM Transactions on
Graphics (Proceedings of SIGGRAPH 2007) 26, 3, 28:1–28:12.
H ANIKA , J., DAMMERTZ , H., AND L ENSCH , H. 2011. Edgeoptimized À-trous wavelets for local contrast enhancement withrobust denoising. Computer Graphics Forum 30, 7, 1879–1886.

NASA, 2005. Blue marble: Next generation. NASA Earth Observatory. http://earthobservatory.nasa.gov/Features/BlueMarble/.
PARIS , S., AND D URAND , F. 2006. A fast approximation of thebilateral filter using a signal processing approach. In Proceedingsof European Conference on Computer Vision (ECCV) 2006, 568–
580.
PARIS , S., H ASINOFF , S. W., AND K AUTZ , J. 2011. Local Laplacian filters: Edge-aware image processing with a Laplacian pyramid. ACM Transactions on Graphics (Proceedings of SIGGRAPH
2011) 30, 4, 68:1–68:12.

H ECKBERT, P. 1989. Fundamentals of Texture Mapping and Image
Warping. Master’s thesis, U.C. Berkeley.

R EEVES , W., S ALESIN , D., AND C OOK , R. 1987. Rendering antialiased shadows with depth maps. In Proceedings of SIGGRAPH
1987, ACM Press / ACM SIGGRAPH, 283–291.

J EONG , W.-K., J OHNSON , M. K., Y U , I., K AUTZ , J., P FISTER ,H., AND PARIS , S. 2011. Display-aware image editing. In
Proceedings of IEEE International Conference on Computational
Photography.

S UMMA , B., S CORZELLI , G., J IANG , M., B REMER , P.-T., AND
PASCUCCI , V. 2010. Interactive editing of massive imagery
made simple: Turning Atlanta into Atlantis. ACM Transactions
on Graphics 30, 2, 7:1–7:13.

K ASS , M., AND S OLOMON , J. 2010. Smoothed local histogram filters. ACM Transactions on Graphics (Proceedings of SIGGRAPH
2010) 29, 4, 100:1–100:10.

TAN , P., L IN , S., Q UAN , L., G UO , B., AND S HUM , H. 2008. Filtering and rendering of resolution-dependent reflectance models.
IEEE Transactions on Visualization and Computer Graphics 14,
2, 412–425.

K AZHDAN , M., AND H OPPE , H. 2008. Streaming multigrid
for gradient-domain operations on large images. ACM Transactions on Graphics (Proceedings of SIGGRAPH 2008) 27, 3,
21:1–21:10.
KOENDERINK , J. J., AND VAN D OORN , A. J. 1999. The structureof locally orderless images. International Journal of Computer
Vision 31, 2-3 (Apr.), 159–168.
KOPF, J., U YTTENDAELE , M., D EUSSEN , O., AND C OHEN , M. F.
2007. Capturing and viewing gigapixel images. ACM Transactions on Graphics (Proceedings of SIGGRAPH 2007) 26, 3,
93:1–93:10.
M ALLAT, S., AND Z HANG , Z. 1993. Matching pursuits with timefrequency dictionaries. IEEE Transactions on Signal Processing
41, 12 (Dec), 3397 –3415.
M ALLAT, S. 2009. A Wavelet Tour of Signal Processing: The
Sparse Way, 3rd ed. Academic Press.

T HOMPSON , D., L EVINE , J., B ENNETT, J., B REEMER , P.-T.,
G YULASSY, A., P EBAY, P., AND PASCUCCI , V. 2011. Analysis of large-scale scalar data using hixels. In Proceedings ofIEEE Symposium on Large-Scale Data Analysis and Visualization(LDAV), 23 – 30.
T OMASI , C., AND M ANDUCHI , R. 1998. Bilateral filtering for grayand color images. In Proceedings of International Conference on
Computer Vision (ICCV) ’98, 839–846.
W ILLIAMS , L. 1983. Pyramidal parametrics. In Proceedings of
SIGGRAPH 1983, ACM Press / ACM SIGGRAPH, 1–11.
YOUNESY, H., M ÖLLER , T., AND C ARR , H. 2006. Improving thequality of multi-resolution volume rendering. In Proceedings of
Eurovis 2006, 251–258.

Figure 12: Median filtering. (a) Original image (320 × 428) with salt and pepper noise. (b) Ground truth 5 × 5 median applied to level 0,then downsampled to level 1. (g) Naive equivalent median computed in level 1. sPDF-map where Wj is a Gaussian of size (c,d,e) 3 × 3,
(h,i,j) 5 × 5. Median from sPDF-map with (c,h) 1 coefficient chunk, (d,i) 2 chunks. (f) Gaussian pyramid level 1 (160 × 214). (e,j) E[Xp ] ofsPDF-map level 1 (top half: 1 coefficient chunk, bottom half: 2 chunks). PSNR [dB] between (g,b): 30, (c,b): 30, (d,b): 36, (h,b): 38, (i,b): 39.

Figure 13: Bilateral filtering. (Top row) Zoom-ins of original 16, 898 × 14, 824 image (a) at level 4: (b) Ground truth bilateral, (c) sPDF-map
bilateral, (d) Naive bilateral. Luminance PSNR [dB] between (c,b) 43, (d,b) 41. (Bottom row) (e) Original image (512 × 512). The sPDF-mapof (e) uses a Wj of size 3 × 3 and two coefficient chunks. Zoom-ins from (f,g,h,i) level 1, (j,k,l,m) level 2. (f,j) Downsampled image, no bilateralfiltering. (g,k) Ground truth bilateral. (h,l) sPDF-map bilateral. (i,m) Naive bilateral. PSNR [dB] between (h,g) 37, (i,g) 35, (l,k) 38, (m,k) 37.