EUROGRAPHICS 2010 / M. Chen and O. Deussen(Guest Editors)

Volume 30 (2011), Number 2

Goal-based Caustics
Marios Papas1

Wojciech Jarosz1

Wenzel Jakob1

1 Disney Research Zürich

Szymon Rusinkiewicz2

2 Princeton University

Wojciech Matusik1

Tim Weyrich3

2 University College London

Abstract
We propose a novel system for designing and manufacturing surfaces that produce desired caustic images whenilluminated by a light source. Our system is based on a nonnegative image decomposition using a set of possibly
overlapping anisotropic Gaussian kernels. We utilize this decomposition to construct an array of continuous surface patches, each of which focuses light onto one of the Gaussian kernels, either through refraction or reflection.
We show how to derive the shape of each continuous patch and arrange them by performing a discrete assignmentof patches to kernels in the desired caustic. Our decomposition provides for high fidelity reconstruction of natural
images using a small collection of patches. We demonstrate our approach on a wide variety of caustic images by
manufacturing physical surfaces with a small number of patches.
Categories and Subject Descriptors (according to ACM CCS):
Generation—Line and curve generation

I.3.3 [Computer Graphics]: Picture/Image

1. Introduction
Recently, increased attention has been paid to the designof physical materials with desired reflectance or refraction properties [WDB∗ 07, WPMR09, MAG∗ 09, DWP∗ 10,
HFM∗ 10, FDL10]. In this paper, we propose a novel systemfor designing and manufacturing surfaces that produce desired caustic images when illuminated by a light source. This
allows, for instance, to create a translucent object that, when
lit by direct sunlight, projects a custom image on a projection plane behind it. (This is different from projecting lightthrough a slide, as light is redirected rather than absorbed,
maintaining higher light levels and contrast.)
A traditional application of reflective caustic design is luminaire design: the creation of lamp reflectors that create adesired radiance distribution for a given light source position. Other applications could be for instance in architecture(pattern glass that casts attractive caustics into a room) or
security features, creating a hard-to-forge artefact that canbe “read out” under normal daylight.
Previous work on physical surface reflectance or refraction design that would be capable of creating caustic patterns
either uses a bottom-up approach, constructing geometry
from discrete building blocks [WPMR09], or uses a global
analysis-by-synthesis approach [PP05, PPV07, FDL10]. The
bottom-up approach suffers from quantization artifacts dueto their choice of primitives, as we will show, while global,
c 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.
Published by Blackwell Publishing, 9600 Garsington Road, Oxford OX4 2DQ, UK and
350 Main Street, Malden, MA 02148, USA.

Figure 1: We design and mill piecewise-smooth arrays of
micropatches which reshape incident light onto a collectionof Gaussian caustics, to approximate arbitrary images at a
projection plane. Above, we show a photograph of a physicalsample projecting the Lena image.
continuous optimizations cannot guarantee convergence.
In this paper, we propose a novel, discrete (bottom-up)
approach for surface generation that band-limits the constructed appearance and thus avoids quantization artifacts
while guaranteeing convergence. In addition, we show that
our surfaces are more amenable to milling than the current
state of the art in caustic generation [FDL10], at the samePapas et al. / Goal-based Caustics

time providing higher contrast.

3. Overview

At the core of our system is a nonnegative image decomposition into (possibly overlapping) anisotropic Gaussian
kernels, the superposition of which yields an aliasing-freereconstruction of an input image. From this decomposition,
we create a surface consisting of an array of continuous surface elements, each of which focuses light onto one of thekernels, forming a Gaussian caustic through either refraction
or reflection. In analogy to the concept of microfacets, wecall these smooth surface elements micropatches. Together,
these micropatches create a joint caustic that approximatesthe input image.

The final output of oursystem is a regular arrayat surfaceof (reflective or refractive)
micropatches, where, under fixed illumination conGaussian
milled surface
ditions, each micropatch caustic
of micropatches
generates a Gaussian caustic and where all elemenprojection
tary caustics sum up toplaneform a desired input image
in a known projection plane. Figure 1 shows a photograph ofthis and we illustrate this setup in the inset above.

We show how to derive the shape of each micropatch
and how to near-optimally arrange such patches so that they
project onto the kernels in the desired caustic. Our decomposition provides for high fidelity reconstruction of natural images even for small collections of patches. We compare theresults of our algorithm to previous work and demonstrate
our approach on a wide variety of caustic images by manufacturing physical surfaces with a small number of patches.

2. Related Work
The process of finding the underlying structure that matches
the high-level reflection or refraction specification havebeen explored by a number of researchers. The computerbased optimization approaches for designing custom reflectors have been investigated since the early 1990s (see, forexample, Neubauer [Neu94]). The most common approach
uses analysis-by-synthesis, where the simulation framework
evaluates the light distribution for a given state and compares it to the desired distribution. Typically the reflector
surface is assumed to be continuous and therefore continuous optimization approaches can be used. In this context, thework by Patow et al. [PP05, PPV07] use SAN optimization
to obtain the reflector surface for a given desired distribution and a fixed light position. As the search space is inherently large and difficult to optimize researchers have used
GPUs to decrease the simulation time [MMP09, FDL10],
have used more efficient surface representations (e.g., BSplines [FDL10], NURBS [MMP09,ASG08]), have reduced
the search space to radially symmetric reflectors [ASG08],
and they have used more efficient optimization strategies (forexample, SPSA [FDL10]).
Our system is the most similar to Weyrich et al.’s approachfor the fabrication of microgeometry [WPMR09], in that
we use a regular assembly of discrete surface elements that
redirect light. We also use a similar optimization framework
to find a near-continuous arrangement of these elements.
In contrast, however, we use higher-order elements (reflectors/refractors with a Gaussian caustic), which leads to better
tonal reproduction, and our optimization has to consider local light sources and local target points when orienting theseelements, as we will discuss.

In the remainder, we present the components of this system. Section 4 describes our nonnegative image decomposition into Gaussian kernels. Section 5 derives the construction of an individual micropatch that generates a desired
Gaussian caustic. Finally, in Section 6 we present a global
optimization scheme to implement the Gaussian image decomposition with micropatches, while considering physicalconstraints of a fabrication process.
The system described here works equally well for reflective and refractive surfaces. For simplicity, we will focus on
the more complex effect of refraction, while the reflective
case follows by analogy.
4. Gaussian Image Decomposition
We approximate images using nonnegative linear combinations of anisotropic Gaussian kernel functions. In the context
of normalized statistical distributions, such representations
are usually referred to as Gaussian mixture models (GMMs).
The parameters to be computed in our application consistof the weights wi , means µi and covariance matrices Σi of a2-dimensional m-term GMM of the following form:m
wi
f (x) = ∑ √exp − 12 (x − µi )T Σ−1(1)i (x − µi )i=0 2π det Σi

Rather than finding the parameters of this model using general nonlinear optimization techniques (as done by Tsaiand Shih [TS06] for approximation environment maps using spherical radial basis functions), we use expectation
maximization (EM) [DLR∗ 77] specialized to GMMs, whichrobustly converges to high-quality approximations using acomparably small number of Gaussians terms.
Each of the summands in Equation (1) eventually leads toa micropatch on the milled surface. To facilitate the manufacturing process, it will be useful to introduce the requirement that all of the patches collect an equal amount of energy, allowing them to be manufactured with identical sizes.
In the above expression, this constraint translates into setting
wi := m1 (i = 1, . . . , m).
Note that instead of directly operating on the input image,the algorithms in this section instead make use of a densec 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

Papas et al. / Goal-based Caustics
Input image

Sampled representation

CCVT parameter estimate

E Step

Converged result

M Step(as discs)(as discs)

Figure 2: We approximate images using a linear combination of anisotropic Gaussian kernels. An initial estimate of theirparameters is computed and iteratively improved using a sparse version of expectation maximization.

(a) m = 256(b) m = 1024(c) m = 4096(d) m = 16384(e) Input image

Figure 3: Our image representation produces high-quality approximations using comparatively few Gaussian terms. The obtained quality continuously improves as more terms are added. In this example, the number of samples xi is held constant, sothat sampling noise becomes visible for large numbers of Gaussians (d).
point representation, which is obtained by sampling a large
number of points proportional to the pixel luminances (weuse n = 4, 000, 000 in our examples, in order for it to beorders of magnitudes larger than m, the number of patches).
In the following, these will be referred to as xi (i = 1 . . . n).
4.1. Initialization using CCVT
It is important to note that while EM is guaranteed to converge, it will generally only find a locally optimal solution
to Equation (1). Hence it is important to supply it with areasonable initial parameter estimate. We experimented withseveral strategies for computing such an estimate, and obtained the most reliable results using regular point sets that
adapt to the image intensity.
Such a point set can easily be created by computing a Capacity-Constrained Voronoi Tesselation (CCVT,[BSD09]). Following that, the centroid of each Voronoi region is chosen as the center of a Gaussian in the initial parameter set. (This bears similarity to the stippling approachby Weyrich et al. [WPMR09]. Note, however, that in our
approach, stipple point locations are only a starting guess toguide the optimization of a more complex image approximation.)
Initially, Gaussian kernels are centered at the stipple points and initialized to an isotropic covarianceΣi = diag(σ 2 , σ 2 ), with σ set to the radius of a k-nearest
neighbor query across sample points xi near its center, andc 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

(a) Initial estimate(b) Isotropic(c) Anisotropic

Figure 4: Starting from the initial estimate, our optimization
can compute either isotropic or anisotropic approximations.
Due to the significant qualitative improvement, we always
use anisotropic Gaussians for our surfaces.
with k equal to the capacity of the associated Voronoi region.
This leads to a similar weighted contribution of samples toeach Gaussian’s region.
4.2. Refinement using EM
The EM algorithm improves the initial estimate using a sequence of alternating E and M steps; the E step computes theprobabilities of each sample conditioned on the current set ofparameters, whereas the M step maximizes the expected log
likelihood of the parameters subject to the previously computed probabilities.
The standard formulation of EM considers all pairs ofGaussians and sample points. In our context, this would leadto prohibitive running times and storage requirements, given
the large numbers of Gaussians and samples involved. In-Papas et al. / Goal-based Caustics

stead, we disregard pairs whose distance exceeds a certain
threshold, on the grounds that our Gaussians are comparatively narrow and approach zero for larger distances. Thiseffectively turns each E and M step into a sparse problem.
By using a kd-tree to look up relevant pairs, we speed up the
EM convergence from O(n2 ) to O(n log n), if we assume mand n to be proportional to each other. In addition, we use asparse matrix representation to pass probabilities into the M
step, to avoid O(n2 ) memory usage.
For the number of sample points and Gaussians used inthis paper, each pair of EM-iterations took approximately
17 seconds on a on a 3Ghz Intel Core 2 Duo machine. In
all cases, we ran between 50 and 100 iterations. The aforementioned sparsity optimizations were crucial to make thefitting process practical, since the O(n2 ) memory costs of
the unmodified method quickly exhausted our machine’s 8
GiB of memory.

(a)Incident Light
Micropatch
Domain

z=0 










Gaussian Intensity
Distribution

Projection Plane











D

z = −D

Compute Normals
Normal Field(b)

Micropatch
DomainIntegrate Normal Field

5. Micropatch Design
We wish to compute the shape of a square surface patch
which produces a 2D Gaussian irradiance footprint on a
projection surface when illuminated with incident lighting.
We call such a surface patch a micropatch. In the reflective
case, this will be a minute concave reflector, in the refractive
setting a plano-convex microlens. We start out by determining the normal field of the desired patch, obtaining the final
shape through subsequent integration. As stated before, weare confining the description to the refractive case.
Assumptions: We assume the configuration illustrated in
Section 3, showing incident lighting, the milled surface, andthe projection plane. Without loss of generality we describe
the case where the projection surface is a plane parallel to themicropatch (at z = −D), but the approach is the same given
an appropriate 2D domain on an arbitrary curved projection
surface. We orient our surfaces to have the milled side facing the projection plane and the flat side facing the light.
For simplicity of exposition, we assume distant illumination
coming from the z+ direction and can therefore ignore thefirst refraction (on the un-milled side of the surface). This is
not a restriction of our method, and we can in fact accountfor local light sources and refractions at both interfaces (bysolving for the first refraction point as done by Walter et.
al [WZHB09]), with the remaining procedure largely unchanged.
Approach: To compute the surface normal of the facet wefirst need a mapping between points on the micropatch, g,
and the corresponding points, p, on the projection plane. To
enforce that the resulting irradiance profile is a 2D Gaussian, the determinant of the Jacobian of the mapping needs
to be a 2D Gaussian. Given this continuous correspondence,
we analytically solve for the normal direction producing therequired refraction direction at any point of the micropatch,
and we then integrate a discretized version of this normal(c)Micropatch
Surface

Figure 5: We (a) define a bijective mapping between points
in the micropatch domain and points on the projection plane,such that the irradiance distribution follows a Gaussian profile. We then (b) analytically compute the surface normals
that refract/reflect the light in this way. We integrate thisnormal field to arrive at the required micropatch surface (c).
field to obtain the elevated micropatch surface. We illustrate
this process in Figure 5.
5.1. Irradiance Warp
We first find a mapping from any point in the micropatch
domain to a zero-mean unit Gaussian at the projection planeq. We subsequently explain how we modify this mapping toaccount for anisotropy and arbitrary means.
Since our micropatch is square, we parametrize positions
as g(s, t) = [gx , gy , gz ]T where (s, t) ∈ [0, 1]2 . We therefore
need a transformation from the uniform square to a standard
2D normal distribution. We could use the Box-Muller transformation [BM58] to map points on the patch to points on
the projection plane. However, we wish this mapping to be“well behaved” (i.e. continuous, low-distortion, and without singularities) to ensure the resulting normal field is integrable. Unfortunately, the original Box-Muller transformdoes not satisfy these requirements. In its standard form it
transforms points on a disk to points on a Gaussian, andmoreover this mapping contains a singularity in the polarradius (by mapping the boundary of the disk to the origin,and the center of the disk to infinity).
To overcome this we use a two-stage transformation.
We first apply Shirley’s low-distortion concentric mapping [SC97] to map points (s, t) in the square micropatchc 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

Papas et al. / Goal-based Caustics
Concentric
MappinguModied BoxMuller TransformsvtMicropatch Surface

Unit Disk

Unit Gaussian

Figure 6: Our multi-stage process of mapping points in the
micropatch domain (s, t) to points on a Gaussian caustic.
domain to points on a unit disk (u, v). We then apply a modified Box-Muller transform, which non-uniformly scales theradius without singularities. Given the point (u, v), we compute the new radius as:r−2 log(1 − u2 − v2 )R=,(2)
u2 + v2
and the resulting point on a unit Gaussian is q = [uR, vR]T .
We illustrate this two-stage mapping in Figure 6.
5.1.1. General Gaussians

We now need to transform our unit Gaussian position q toobey a specific mean µ and covariance matrix Σ. Standard
approaches for doing this are typically only concerned withproperly transforming the distribution, by e.g. interpretingthe coordinates as offsets along the scaled eigenvectors of Σ.
This results in a sequence of scale, rotation and translation
operations. Unfortunately, using this approach can introduce
large distortions in the mapping from patch coordinates tothe final projection plane positions, affecting the integrability of the normal field. To avoid this, we eliminate rotationfrom our transformation, and instead directly scale the point
q along the directions of the eigenvectors of Σ and then translate the point by the mean µ . This is equivalent to transforming q by the square root of the covariance matrix (as definedby the matrix logarithm):√(3)p = − Σ q + µ.
We complete our mapping by placing this 2D point onto theprojection plane, which involves appending a z-componentof pz = −D to the above result.
In Equation (3), we actually scale by the negative squareroot, which introduces a reflection about the origin prior totranslation. This results in micropatches which always focus
light (between the micropatch and the projection plane) regardless of the size of the Gaussian caustic. This gives ourrefractive micropatches more curvature, which is better for
milling.
5.2. Normal Computation
Given this mapping, we need to compute the surface normal at each location g which would refract light towards p.
c 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

Since the refracted direction, surface normal, and incident direction must lie in the same plane,
we view this as a 2D problem inthis coordinate system, which we
illustrate on the right.

−"n
θi








g
"n

θ
The angle between the outgo- Dθb oing direction and the incident di- 





rection is θb = arctan (β ) wherepd
β = D . By using the fact that% &' (dθo = θi − θb and combining with
Snell’s Law, we can solve for the incident angle asarctan(ηo sin(θb ))
θi =
,
(4)
ηo cos(θb ) − ηiwhere ηi and ηo are the incident and outgoing indices of refraction. We obtain the surface normal vector by rotating theincident vector towards the outgoing vector by θi radians.

5.3. Normal Field Integration
The remaining step follows the height field optimization of
Weyrich et al. [WPMR09]: similarly, we model the micropatch as a rectilinear grid of elementary microfacets withthe desired normal orientations, solving a Poisson problem
to lift each facet onto a level that leads to a near-continuous
surface. In contrast to the previous work, we do not optimize
for tileability, as only an individual micropatch is modeled.
Note that we also skip the “shuffling” stage, as, by design,
our mapping already leads to a continuous and optimal facet
layout.
In our experiments, we found that for a wide range oftarget Gaussians, the residual of the optimization to be diminishingly small, suggesting that the input normal fields arevery close to integrable. In fact, the error is far below the precision of our manufacturing tools, and subsequent polishingof our samples leads to smooth patches with near-Gaussiancaustics. We show visualizations of individual micropatches
in Figure 7.
6. Micropatch Assembly
Given a Gaussian image decomposition, we aim at finding an
assembly of micropatches that generate the respective Gaussian caustics, while building a surface that is as continuous

Figure 7: Refractive micropatches producing (left) an
isotropic Gaussian caustic, and (right) a highly anisotropic
Gaussian caustic.

Papas et al. / Goal-based Caustics

A naïve approach would be to pre-generate micropatches
for all Gaussian kernels in the image, and to treat their layout
as an independent problem. Unfortunately, since we accountfor local lights and projection distances, the relative positionof each micropatch with respect to its target Gaussian affects
its shape, rendering an independent layout optimization impossible. Conversely, including a recomputation of each micropatch’s shape into a global optimization is costly at best,
and would negatively affect scalability to larger problems.
Instead, we propose an optimization scheme that approximates each patch by a planar proxy geometry that refracts
a central ray toward the center of the target Gaussian. We
found that the orientation of such a proxy is a good approximation of the overall “tilt” of a Gaussian micropatchthat projects toward the same center. We use this planar approximation during optimization to gauge compatibility ofneighboring patches, in the sense of Equation (1) in Weyrich
et al. [WPMR09]. This approximation puts us once again
into the same domain as their microfacet optimization, thistime each microfacet acting as the proxy for one micropatch.
Accordingly, we use a similar optimization scheme, starting with a permutation optimization, followed by a Poisson
integration to solve for facet heights. After the last step, we
replace each microfacet by a micropatch that aims at the corresponding Gaussian kernel.
6.1. Facet Permutation
Due to the similarity of this pipeline, we only point out differences to this previous work, using the same notation asin Weyrich et al. [WPMR09], referring to the previous paper
for more detail.
Permutation Loop. We, too, deploy simulated annealing(SA) as a very simple algorithm to optimize the permutationof facets. Finckh et al. [FDL10] have pointed out the superiority of simulated perturbation stochastic approximation(SPSA) over SA, but we found the simpler SA iteration to
produce satisfactory results nevertheless. The central transition operation is still the swap of two facets, but in contrast
to the case of BRDF design, where only the net effect of a
microfacet assembly under infinite lighting and viewing isof interest, a swap alters the orientation of the two facets:
swapping two facets basically corresponds to swapping theassignment of two kernel centers to their facets, triggering
the recomputation of each facet’s normal.
Energy Function. The previous approach minimizes thefollowing constraints: facet compatibility C, tileability I, andthe avoidance of valleys V . In our setting, tileability is not
required, and valley avoidance is of reduced value, as the(refractive) micropatches create slight ridges at all times. In

order to bias the solution toward more low-frequency solutions, however, we add a smoothness term S = Sx + Sy with
Sx = ∑ ∑yxdz(x, y)
dx

2

,

(5)

with dz defined as in Weyrich et al. [WPMR09], and similarly for Sy . This leads to surfaces that are easier to polishafter the milling process (see Section 7). The final energyfunction is then:
C+S
E=
,
(6)
m2
with m the total number of micropatches. Figure 8 shows theenergy convergence plots for the samples used in this paper.
−5

x 10
3

Venus
EinsteinGreat Wave
Lena

2.5
Energy (E)as possible, in order to facilitate fabrication and to minimize
stray light at boundaries between patches.

2
1.5
1
0.5
0

100

200

300
400
500
Temperature Changes

600

700

Figure 8: Convergence plots of the energy minimization using simulated annealing.
6.2. Height Optimization
The height field integration to determine facet placements
in z direction follows the previous approach, except that weonce more ignore tileability. In our setting, however, changing the height level of a microfacet requires re-aiming it at
its caustic center. We consequently run three iterations ofheight optimization and normal recomputation, which typically leads to reliable convergence after the second iteration.
After this, we replace each facet by the corresponding micropatch, leading to near-continuous surfaces (see Figure 9).
7. Results
We evaluate our refractive designs both through digital renderings (Figures 10 and 14), as well as by milling physicalsamples from PMMA (trade name Plexiglas).
Milling Procedure. We use a Roland DGA EGX-600 Engraver with a maximum cutting area of 61 × 40 × 4.2 cm,
maximum software spatial resolution 0.01 mm, and a repeat accuracy of 0.05 mm. The engraver uses a C2-171010K Carbide Tipped Engraving Cutter (0.13 mm Tip and4.36 mm Shank).
For smoothest surface results, we use one draft and one
fine milling pass, each time along the same axis (x or y).
We run the draft mode with a 0.16 mm path spacing, whichcuts off layers of maximally 2 mm thickness per pass. This
helps to prevent damaging the milling bit. This pass leaves
c 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

Papas et al. / Goal-based Caustics

Lena (refractive)

Einstein (refractive)

Venus de Milo (refractive)

Great Wave (refractive)

Great Wave (reflective)

Figure 9: Visualizations of the micropatch assemblies used
to produce the results in this paper.
0.2 mm of material for the final pass at 0.02 mm path spacing. For our surfaces (10×10 cm), the draft pass takes between 12 and 24 hrs, while the fine mode: requires another24 to 48 hrs, depending on the depth variation of the surface.
Afterwards, we (manually) polish the surface using Novus
Plastic Clean & Shine No. 2, with a final pass with No. 1.
Polishing takes about 1 hr.
Physical Samples. Figure 12 shows two sample datasets
milled according to this procedure. Each of the figures shows
the target image, a decomposition into 1024 Gaussian kernels, a ray-traced simulation of the resulting micropatch assembly, and a photograph of the caustic generated by the real
sample. The samples are 9×9 cm in size, comprising 32×32
micropatches of 2.8 mm edge length. We image the caustics
by using a white-light projector as an approximation of a
distant light source. In the examples, the physical projection
surface was at about 25 cm behind the sample.

the polishing procedure. Great care has to be taken not todamage the material during manual polishing; especially theridges between micropatches are hard to reach and may leadto spurious caustic patterns as particularly visible in Figure 12. But also the milling procedure may lead to imperfections. Systematic error sources, such as stepper quantization,
play a role, and minute shifting of the sample over the long
milling time and local overheating of the material tend to
impair the quality of the results. Given that, we believe that
our results suggest that our approach is robust enough forphysical applications.
7.1. Comparison to Continuous Optimization
Finckh et al. recently presented an algorithm to optimize
a single, continuous surface to exhibit a custom caustic [FDL10]. The authors kindly agreed to run our “Lena”
test specification through their algorithm to allow for a direct comparison to our approach. While their system readily
converged, the result does not span the full contrast range ofthe rather complex target image (see Figure 11). The surface
furthermore turned out extremely flat (0.07 mm depth range
over a 100×100 mm area), which made it very susceptible tostepper quantization of our mill. Figure 11c shows a simulation of the quantization introduced by the mill’s z resolution.
We hence decided to abort the experiment of creating a physical surface from this data set.
While there may be alternative target specifications(shorter throw distance or larger projection) that would stilllead to millable surfaces, we cautiously conclude that ouralgorithm outperforms Finckh et al.’s global optimization interms of physical viability of the results. We hope to be able
to investigate this further to achieve a fair comparison.
8. Discussion and Future Work
Despite other potential applications, our evaluation focuseson image projections onto a plane, which is arguably themost convenient domain to gauge the quality of the method’s
output. It, however, invites comparison with a much easier
method for image generation—that of slide projection. In
contrast to slide projection, our method does not modulate,

As the examples show, our technique translates well tophysical constraints, albeit not without artifacts. Differences
between the simulation and the captured results are due toa number of error sources, of which the largest is arguably(a) Surface(a) Input image(b) Decomp.

(c) PBRT sim.

Figure 10: Application to an image of Venus de Milo.
c 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

(b) PBRT sim.

(c) Isocontours

Figure 11: We visualize the surface (a) generated using themethod by Finckh et al. [FDL10], as well as a simulation of
the resulting caustic (b). Unfortunately this surface poses achallenge for physical reproduction due to its small verticalextent, which results in only a small number of quantizedisocontours (c).

Papas et al. / Goal-based Caustics(a) Input image(b) Gaussian decomposition(c) PBRT Simulation(d) Physical capture

Figure 12: Our method applied to the Lena, Einstein, and Venus images. The Gaussian decomposition is shown in (b) The
resulting surfaces are visualized in Figure 9, and we show a PBRT simulation of the caustic in (c). A photograph of the resulting
caustic from a milled reproduction of the surface is shown in (d).
but only redirect light, avoiding light loss due to absorption.
This is similar in spirit to the light reallocating projectionengine of Hoskinson et al. [HSHF10] and meets the hardrequirement in luminaire design, that all energy should be
used for lighting.
One limitation of our current approach is the restriction toequal-sized surface patches, a choice made to simplify thelayout of micropatches. This implies roughly equal energy
for all Gaussians (ignoring Fresnel effects), which leads tothe side effect that brighter regions receive more, and higherfrequency splats than dark regions, which means that brightregions show more detail than darker regions (visible, forinstance, in the Venus dataset). This effect could be avoided
by allowing for micropatches of different sizes, at the costof having to solve the packing problem of seamlessly andoptimally arranging patches of various sizes.
Another limitation is the fact that the generated surface focuses the caustic on a projection plane with a fixed distance.

As shown in Figure 13 the generated caustic is warped anddistorted when the actual distance deviates from the optimal.
It is further conceivable to extend the caustic generation from a single plane to a two-plane target specification. This would enter the realm of light field construction and has been demonstrated as an initial prototype by
Fuchs et al. [FRSL08]. Similarly, it is worth investigating
whether the refractive domain allows for lighting-dependentimage formation, as shown by Alexa et al. for diffuse reflectors [AM10].
9. Conclusion
In this paper we proposed a novel system for designing andmanufacturing surfaces that produce desired caustic images
when illuminated by a light source. The system produces
results superior to previous work on goal-based caustic design and proves to be robust even for complex caustics. Our
nonnegative image decomposition scales well to very largec 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.

Papas et al. / Goal-based Caustics

Closer than Focusing Plane←−−→At Focusing Plane←−−→

Further than Focusing Plane

Figure 13: Positioning the sample at different distances to the projection plane shows the effective depth of field of the caustic.

(a) Input image(b) Decomposition(c) PBRT simulation (refractive)(d) PBRT simulation (reflective)

Figure 14: Application to the woodcut “The Great Wave off Kanagawa” by Hokusai. We show the input image (a), the Gaussian
decompositions (b), and the resulting caustic images simulated using PBRT for (c) a refractive and (d) reflective surface.
numbers of Gaussian kernels, paving the way for highly intricate caustic designs.
At the same time, we show that our anisotropic Gaussianprimitives are flexible enough to describe a wide range ofnatural image content even with only moderate numbers ofkernels, making it possible to fabricate functioning prototypes with off-the-shelf hardware, as we showed with ourrefractive designs.
We hope that our approach of using continuous functions
as discrete building blocks for physical reflectance design
will allow for future developments in creating more complex
lighting and reflectance phenomena.
References[AM10] A LEXA M., M ATUSIK W.: Reliefs as images. ACM
Trans. on Graphics (Proc. SIGGRAPH) 29 (2010), 60:1–7. 8[ASG08] A NSON O., S ERON F. J., G UTIERREZ D.: Nurbs-based
inverse reflector design. In In: Proceedings of CEIG 2008 (2008),
pp. 65–74. 2[BM58] B OX G. E. P., M ULLER M. E.: A note on the generationof random normal deviates. The Annals of Mathematical Statistics 29, 2 (1958), 610–611. 4[BSD09] BALZER M., S CHLÖMER T., D EUSSEN O.: Capacityconstrained point distributions: A variant of Lloyd’s method.
ACM Trans. Graph. (Proc. SIGGRAPH) 28, 3 (2009), 86:1–8. 3[DLR∗ 77] D EMPSTER A., L AIRD N., RUBIN D., ET AL .: Maximum likelihood from incomplete data via the EM algorithm.
Journal of the Royal Statistical Society. Series B (Methodological) 39, 1 (1977), 1–38. 2[DWP∗ 10]

D ONG Y., WANG J., P ELLACINI F., T ONG X., G UO
B.: Fabricating spatially-varying subsurface scattering. ACM
Trans. on Graphics (Proc. SIGGRAPH) 29, 4 (2010), 1–10. 1[HFM∗ 10] H AŠAN M., F UCHS M., M ATUSIK W., P FISTER H.,
RUSINKIEWICZ S.: Physical reproduction of materials withspecified subsurface scattering. ACM Trans. on Graphics (Proc.
SIGGRAPH) 29, 4 (2010), 61:1–10. 1[HSHF10] H OSKINSON R., S TOEBER B., H EIDRICH W., F ELS
S.: Light reallocation for high contrast projection using an analog
micromirror array. To appear in ACM Trans. on Graphics (Proc.
SIGGRAPH Asia) 29, 5 (2010). 8[MAG∗ 09] M ATUSIK W., A JDIN B., G U J., L AWRENCE J.,
L ENSCH H. P. A., P ELLACINI F., RUSINKIEWICZ S.: Printing
spatially-varying reflectance. ACM Trans. on Graphics (Proc.
SIGGRAPH Asia) 28, 5 (2009), 128:1–9. 1[MMP09] M AS A., M ARTÍN I., PATOW G.: Fast inverse reflector design (fird). Computer Graphics Forum 28, 8 (Dec. 2009),
2046–2056. 2[Neu94] N EUBAUER A.: The iterative solution of a nonlinear inverse problem from industry: design of reflectors. In In: Proceedings of the International Conference on Curves and Surfaces inGeometric Design (1994), A. K. Peters, Ltd., pp. 335–343. 2[PP05] PATOW G., P UEYO X.: A survey of inverse surface designfrom light transport behavior specification. Computer Graphics
Forum 24, 4 (Dec. 2005), 773–789. 1, 2
[PPV07] PATOW G., P UEYO X., V INACUA A.: User-guided inverse reflector design. Computers & Graphics 31, 3 (2007), 501–
515. 1, 2[SC97] S HIRLEY P., C HIU K.: A low distortion map between disk
and square. Journal of Graphics Tools 2, 3 (1997), 45–52. 4[TS06] T SAI Y., S HIH Z.: All-frequency precomputed radiancetransfer using spherical radial basis functions and clustered tensor approximation. ACM Trans. Graph. 25, 3 (2006), 967–976. 2[WDB∗ 07] W EYRICH T., D ENG J., BARNES C., RUSINKIEWICZS., F INKELSTEIN A.: Digital bas-relief from 3D scenes. ACM
Trans. Graph. (Proc. ACM SIGGRAPH) 26, 3 (2007), 32:1–7. 1[FDL10] F INCKH M., DAMMERTZ H., L ENSCH H.: Geometry
construction from caustic images. In Proc. of European Conference on Computer Vision (ECCV) (Sep 2010). 1, 2, 6, 7

[WPMR09] W EYRICH T., P EERS P., M ATUSIK W., RUSINKIEW ICZ S.:
Fabricating microgeometry for custom surface reflectance. ACM Trans. on Graphics (Proc. SIGGRAPH) 28, 3(2009), 32:1–6. 1, 2, 3, 5, 6

[FRSL08] F UCHS M., R ASKAR R., S EIDEL H.-P., L ENSCH H.
P. A.: Towards passive 6d reflectance field displays. ACM Transactions on Graphics (Proc. SIGGRAPH) 27, 3 (2008), 58:1–8. 8[WZHB09] WALTER B., Z HAO S., H OLZSCHUCH N., BALA K.:
Single scattering in refractive media with triangle mesh boundaries. ACM Transactions on Graphics 28, 3 (aug 2009). 4

c 2010 The Author(s)Journal compilation c 2010 The Eurographics Association and Blackwell Publishing Ltd.