Volume 35 (2016), Number 7

Pacific Graphics 2016
E. Grinspun, B. Bickel, and Y. Dobashi(Guest Editors)

An Error Estimation Framework for Many-Light Rendering
K. Nabata1 , K. Iwasaki1,3 , Y. Dobashi2,3 , and T. Nishita3,4
1 Wakayama University, Japan
2 Hokkaido University, Japan

3

4

UEI Research, Japan
Hiroshima Shudo University, Japan

5%+

(a)

(b)

(f)

(d)

0%

(c)

(g)

(e)

Figure 1: Rendering results of (a) our method, relative error images of (b) our method and (c) Lightcuts [WFA∗ 05] in Sponza scene. Relative
error threshold ε = 2% and confidence level α = 95% are specified. 92.96% pixels satisfy that the relative error is within 2% in our method,while only 43.67% pixels satisfy the condition in Lightcuts. In (d) Kitchen scene and (e) San Miguel scene, Cook-Torrance BRDFs andAshikhmin-Shirely BRDFs are used, which cannot be used in Lightcuts. (f) and (g) show relative error images of (d) and (e), respectively.
Abstract
The popularity of many-light rendering, which converts complex global illumination computations into a simple sum of theillumination from virtual point lights (VPLs), for predictive rendering has increased in recent years. A huge number of VPLs
are usually required for predictive rendering at the cost of extensive computational time. While previous methods can achievesignificant speedup by clustering VPLs, none of these previous methods can estimate the total errors due to clustering. Thisdrawback imposes on users tedious trial and error processes to obtain rendered images with reliable accuracy. In this paper,
we propose an error estimation framework for many-light rendering. Our method transforms VPL clustering into stratified
sampling combined with confidence intervals, which enables the user to estimate the error due to clustering without the costly
computing required to sum the illumination from all the VPLs. Our estimation framework is capable of handling arbitrary
BRDFs and is accelerated by using visibility caching, both of which make our method more practical. The experimental results
demonstrate that our method can estimate the error much more accurately than the previous clustering method.
Categories and Subject Descriptors (according to ACM CCS): I.3.3 [Computer Graphics]: Picture/Image Generation—Display
algorithms I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Raytracingc 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John
Wiley & Sons Ltd. Published by John Wiley & Sons Ltd.

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering

1. Introduction
Efficient rendering methods to predict visual appearances are necessary for many applications including design software (e.g. Autodesk 360 [AUT13]) and in the film industry. For such applications, designers and manufacturers need to rely on rendered images
to design, create, and manipulate products, and thus the rendered
images are expected to predict the visual appearance faithfully.
In recent years, many-light rendering [DKH∗ 13] has attracted
attention due to its simplicity, efficiency, and wide applicability.
Many-light rendering simplifies complex global illumination computations by using a large number of virtual point lights (VPLs).
The incident illumination onto points to be shaded (we refer tothese points as shading points) is calculated by using VPLs. The
outgoing radiance from a shading point is the sum of the contributions from all VPLs. For many-light rendering, the accuracy of
the rendering increases as the number of VPLs increase. Therefore,
an exact solution that sums up the contributions from all the VPLs
is computationally expensive. To address this problem, clustering
methods [WFA∗ 05,WABG06,WKB12,HPB07,OP11,BMB15] for
VPLs have been proposed with which efficient, scalable renderingwith a huge number of VPLs can be provided. These methods cluster similar VPLs and the representative VPL for each cluster is used
to calculate the outgoing radiance at each shading point. Unfortunately, none of these methods proposes an error estimation framework for many-light rendering, which results in tedious trial anderror processes to produce images whose errors from the reference
solution are guaranteed to be less than user-specified thresholds forpredictive rendering.
We propose an error estimation framework for many-light rendering by using confidence intervals. Our method builds upon
scalable clustering approaches (e.g. Lightcuts [WFA∗ 05]) to estimate the error of the outgoing radiance directly, while the previous methods [WFA∗ 05, BMB15] can only estimate the maximum error of each cluster instead of that of the outgoing radiance.
By specifying the relative error threshold and a confidence interval parameter, our method automatically partitions a huge number
of VPLs into clusters so that the outgoing radiance estimated bysampling two representative VPLs from each cluster is within thestochastic bound of the true value. Contrary to the previous methods [WFA∗ 05, BMB15] that limit BRDFs whose upper bounds
within the VPL cluster can be computed cheaply and tightly, ourmethod can handle arbitrary BRDFs whose upper bounds are estimated by using BRDF importance sampling. Furthermore, ourmethod introduces a visibility caching method to accelerate VPL
rendering. Unlike previous caching methods [GKPS12, YNI∗ 15]
that record the visibilities between sparsely sampled points (calledcache points) and all the VPLs, our method records the average visibilities between the cache points and the clusters of VPLs, whose
number is much smaller than that of the VPLs, which reduces boththe computational time and memory footprint required for the visibilities at the cache points.
The contributions of our method are as follows.
• An error estimation framework for VPLs is proposed. Our
framework can be easily incorporated into Lightcuts [WFA∗ 05].
• Our method removes the limitation of BRDFs applicable toLightcuts, and thus widens the applicability of VPL rendering
further.
• A visibility caching method is introduced to accelerate rendering
without noticeable decrease in the error estimation accuracy.

2. Previous Work2.1. Many-light Rendering Methods
Recent work on many-light methods [DKH∗ 13, KGKC13] has revealed that global illumination effects can be well approximated byusing many virtual light sources. Keller introduced the instant radiosity method, which calculates the indirect illumination by using
virtual point lights (VPLs) [Kel97]. Following this seminal work, anumber of many-light methods that can handle highly glossy materials [HKWB09, DKH∗ 10, SHD15] and can achieve interactive
rendering [RGK∗ 08, Tok15], have been proposed.
Clustering VPLs: Walter et al. proposed a scalable solution to
many-light methods using a hierarchical representation of VPLs,
called Lightcuts [WFA∗ 05, WABG06, WKB12]. Hasan et al. represent the many-light method as a large matrix, and explore thematrix structure by using row-column sampling [HPB07]. Ou andPellacini [OP11] cluster the shading points into groups calledslices and perform row-column sampling for each slice. Wanget al. [WHY∗ 13] proposed an out-of-core many-light rendering
method. Huo et al. [HWJ∗ 15] proposed a matrix sampling-andrecovery approach for many-light rendering. Bus et al. [BMB15]
cluster the product-space of all shading point-VPL pairs andachieve speedup compared to Lightcuts.
Sampling VPLs: Wang and Akerlund presented a BRDF importance sampling method for VPLs [WA09]. Georgiev et al. proposed an importance sampling method for VPLs by recording thecontributions of VPLs at each cache point [GKPS12]. Wu andChuang [WC13] proposed the VisibilityCluster, which approximates the visibilities between each cluster of VPLs and those of
the shading points by estimating the average visibility. Yoshida etal. [YNI∗ 15] improved the importance caching [GKPS12] by adaptively inserting cache points.
Although many methods have been proposed for many-lightrendering, an error estimation framework has not been proposed.
Lightcuts attempts to control the error due to clustering by estimating the upper bound error of each cluster, but cannot control the
total error from all clusters. This results in tedious trial and error
processes to obtain rendered images whose errors from the reference solution are bounded.
Visibility Caching: Clarberg et al. proposed a visibility caching
method for direct illumination by storing 2D visibility maps at eachcache point [CAM08]. Popov et al. quantize the visibility functions between two points by those between two clusters [PGSD13].
Ulbrich et al. presented a progressive visibility caching approach
using shadow maps [UNRD13]. Although these methods can accelerate rendering significantly by using visibility caching, they
introduce errors due to quantization [PGSD13] and discretization [CAM08, UNRD13], which makes the error estimation more
complex.
c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering

2.2. Rendering with Confidence Interval
Purgathopher [Pur87] introduced the confidence intervals to control the number of samples for distributed ray tracing. Tamstorf andJensen [TJ97] proposed an adaptive sampling method for path tracing using confidence intervals. Rigau et al. [RFS03] presented refinement criteria using convex functions. Hachisuka et al. [HJJ10]
proposed an error estimation framework for progressive photon
mapping. Dammertz et al. [DHKL10] introduced a hierarchicalautomatic stopping criterion for Monte-Carlo methods. Moon etal. [MJL∗ 13] used confidence intervals to detect homogeneous pixels for image denoising. These methods, however, do not estimate
the errors of many-light rendering.

Algorithm 1 Our Algorithm (Input : ε and α, Output : L̂). Procedures of our error estimation framework are highlighted with underlines, and the rest of procedures are identical to Lightcuts.
1: Priority Queue : Q ← C1 , Cut : C ← C1
2: Estimate L̂, ∆L, and σ1√
3: while ∆L > εL̂ ∨ 2σi > εL̂ (∀Ci ∈ C) do
4:
Ck ← Q.pop()5:
Replace Ck in C with two children CL and CR
6:
Estimate L̂CL , L̂CR , σL , and σR
7:Update L̂ and ∆L
8:
Q ← CL , CR
9: end while

2.3. Preliminary of Lightcuts
Since our method builds upon the scalable clustering approaches,
we briefly review Lightcuts. In many-light rendering, the outgoingradiance L(x, xv ) at the shading point x towards the viewpoint xv is
calculated from the following equation.
Nvpl

L(x, xv ) = ∑ I(yi ) f (xv , x, yi )G(x, yi )V (x, yi ),

(1)

i=1

where Nvpl is the number of VPLs, I(yi ) is the intensity of the ith VPL yi , f (xv , x, yi ) is the BRDF, G(x, yi ) is the geometry term,and V (x, yi ) is the visibility term that returns 1 if x and yi are mutually visible and returns 0 otherwise. To render realistic images, alarge number of VPLs are usually required. In this case, summing
the VPL contribution I(yi ) f (xv , x, yi )G(x, yi )V (x, yi ) for all VPLs
is computationally expensive. Rather than summing all VPL contributions, Lightcuts partitions similar VPLs into N clusters, andthe outgoing radiance LCi (x, xv ) illuminated by VPLs in cluster Ci
is calculated from
LCi (x, xv ) = ∑ I(y) f (xv , x, y)G(x, y)V (x, y).

(2)

y∈Ci

Instead of summing all the VPLs in cluster Ci , LCi is estimated bysampling NS VPLs with probability p using the following equation.
L̂Ci (x, xv ) =1 NS I(ys ) f (xv , x, ys )G(x, ys )V (x, ys ).
∑NS s=1p(ys )(3)

Lightcuts samples the VPLs based on the probability proportional
to the intensity (i.e. p(ys ) = I(ys )/ ∑y∈Ci I(y)) as in the following
equation.
1 NS I(ys ) f (xv , x, ys )G(x, ys )V (x, ys )L̂Ci (x, xv ) =∑NS s=1I(ys )/ ∑y∈Ci I(y)=Ii NS
∑ f (xv , x, ys )G(x, ys )V (x, ys ),NS s=1(4)

where Ii = ∑y∈Ci I(y) is the sum of intensity, and Lightcuts sampled
one representative VPL (i.e. NS = 1 is used). The standard variation
σi for Ci is calculated using the following equation (arguments areomitted for brevity and p = I/Ii is substituted).
ppσi = Var[I f GV /p] = Ii Var[ f GV ].
(5)Lightcuts clusters VPLs for each shading point and represents
VPLs and clusters with a binary tree, referred to as light tree. In
c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

the light tree, each leaf node corresponds to each VPL and each inner node corresponds to the cluster that contains all the VPLs of thedescendant nodes. A set of clusters is represented by cut, where every path from the root node to each leaf node contains one node of
the cut. From the root node, cut C is calculated by subdividing eachcluster until each cluster Ci in C satisfies ECi < εL̂(x, xv ), where ε
is a threshold specified by the user, and ECi is the upper bound of
the clustering error calculated from the following equation.
ECi = Ii fub (xv , x, Ci )Gub (x, Ci )Vub ,(6)

where fub , Gub , and Vub represent the upper bound of each function
within the cluster Ci , respectively. The upper bound of visibility Vub
is set to 1. The estimated outgoing radiance L̂(x, xv ) is calculated byL̂(x, xv ) = ∑N
i=1 L̂Ci (x, xv ).
Although Lightcuts can render realistic images efficiently, it has
some limitations as follows. Firstly, Lightcuts do not directly estimate the total error ||L(x, xv ) − L̂(x, xv )|| (i.e. the sum of errors ofeach cluster Ci ). Although it would be possible to estimate the total
error by improving Lightcuts to sum the upper bound of clustering errors as ∑N
i=1 ECi < εL̂, this overestimates the total error andmay lead to a numerous number of clusters and expensive computational time. Secondly, applicable materials in Lightcuts are limited to Lambertian, Blinn-Phong, and isotropic Ward BRDFs only,
since the calculation method of the upper bound fub is currently
developed only for these three BRDF models.
3. Error Estimation Framework for VPLs
3.1. Overview
Fig. 2 and Algorithm 1 show an overview and the pseudo codeof our method, respectively. Our goal is to estimate the error kL̂(x, xv ) − L(x, xv )k by clustering and sampling VPLs without evaluating L. Our method represents VPLs with a lighttree [WFA∗ 05] and the cluster in the light tree is referred to as a
VPL cluster. From the root node of the light tree, our method estimates L̂ by sampling a small number of VPLs from each VPL
cluster Ci in cut C, and determines whether the following criterion
is satisfied.
kL̂ − Lk < εL̂,

(7)

where ε is the relative error threshold specified by the user, where
the arguments are omitted for the simplicity. The criterion, how-

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering
light source(a)light source

(b)

VPL

light source

(c)

(d)

Figure 2: Overview of our method. (a) Root node of VPL cluster that encompasses all VPLs. (b) Light tree. Each leaf node corresponds to
VPL yi and each inner node corresponds to VPL cluster Ci . From the root node C1 , two VPLs are sampled and L̂ and ∆L are calculated. (c)If L̂ does not satisfy ∆L < εL̂, C1 is replaced with two child nodes C2 and C3 . Since y3 and y5 are evaluated in the parent node C1 , they arereused and only y2 and y8 are sampled for C2 and C3 , respectively. Then L̂ is estimated by summing L̂C2 and L̂C3 . (d) If L̂ does not satisfy
the criterion, standard variations σ2 and σ3 are calculated and VPL cluster with the largest standard variation, C2 in this case, is replacedwith C4 and C5 . These processes are repeated until L̂ satisfies the criteria.

ever, contains the unknown (or computationally expensive to evaluate) value of L. To address this problem, our method employs aconfidence interval [L̂ − ∆L, L̂ + ∆L] that contains the true value Lwith probability α which is specified by the user. That is, the probability Pr that the error is less than a tolerance ∆L is α as specified
by the following equation.
Pr (kL̂ − Lk ≤ ∆L) = α,(8)where the tolerance ∆L is calculated by using the tdistribution [Tho12]. By using the confidence interval, thecriterion is calculated by ∆L < εL̂ with probability α. To use theconfidence interval, our method imposes an additional criterion
described in Sec. 3.2. Until the criteria are satisfied, the VPL
cluster with the maximum standard variation is selected from thecut nodes and is replaced with its two child nodes.
3.2. Clustering VPLs
To achieve both scalable rendering and efficient sampling for alarge number of VPLs, our method employs a stratified sampling,which partitions VPLs into clusters (or strata) of similar VPLs andsamples representative VPLs from each cluster. In stratified sampling, if the samples in each cluster follow a normal distribution
and the number of samples for each cluster is allocated by using
the Neyman allocation [Tho12], statistics T calculated by the following equation follows a t-distribution with (n − N) degrees of
freedom.
rn−NL̂ − L
q
T=
,
(9)
n
2 /n
s
∑N
i
i=1 i
where N is the number of clusters, ni is the number of samples forcluster Ci , n is the total number of samples (i.e. n = ∑Ni=1 ni ), ands2i is the sample variance of Ci . Based on Eq. (9), the tolerance ∆L
is calculated using the following equation.
v
uu n N s2
∆L = tα t(10)∑ nii ,
n − N i=1
where tα is the α quantile of theR t-distribution t(x) with (n − N)tαdegrees of freedom that satisfies −tt(x)dx = α.
αThe rendering efficiency depends on how the VPL clusters are
partitioned and how the samples ni for each VPL cluster are allocated when the total number of samples n is fixed. To this end,
we propose an efficient solution for clustering VPLs and allocating
samples for each VPL cluster. In stratified sampling, it is known
that partitioning into many strata and estimating with a small number of samples is more efficient than a few strata with many samples [PH10]. The variance of estimator in stratified sampling is thesum of variances of each stratum. Therefore, selecting the stratumwith maximum standard variation and subdividing the stratum isefficient to reduce the variance. As shown in Eq. (10), the sample
variance s2i of each VPL cluster is required to compute the confidence intervals. Since the minimum number of samples to computethe sample variance is two, ni should be larger than two. As mentioned before, the tolerance ∆L is calculated by using the Neyman
allocation, which is the optimal method for allocating samples toeach cluster under a fixed number of samples n, as shown in thefollowing equation.
ni =

σi
∑N
k=1 σkn,(11)

where σi is the standard variation of VPL cluster Ci . By subdividing the VPL cluster with the maximum standard variation repeatedly, the standard variations of VPL clusters approach uniform andthe numbers of samples for each VPL cluster also approach uniform as shown in Eq. (11). Then two samples (ni = 2) are allocatedfor each VPL cluster and n is calculated as 2N. By substituting
ni = 2 and n = 2N, the tolerance ∆L simplifies to the following
equation.
vuNu∆L = tα t ∑ s2i .
(12)
i=1

Standard Variation Estimation: To select the cluster to subdivide, the standard variation σi for Ci is required. Although σi canbe estimated by using the sampled variance s2i , the estimated standard variation seems poor since only two samples are used. Instead,our method estimates σi using Eq. (5). We estimate the variance
Var[ f GV ] by using the estimated upper bound of BRDF withincluster Ci , fˆub (xv , x, Ci ), and the upper bound of the geometry term
c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering
0.5

(a)

(b)

(c)
normal

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0
-5

Figure 3: Estimate the upper bound of BRDF fˆub by using (a)f (xv , x, y1 ) and f (xv , x, y2 ) which are evaluated at two VPL samples y1 and y2 , (b) maximum value f max based on BRDF importance sampling, and (c) BRDF fωCi at direction ωCi that bounds
minimum angle between the normal and the VPL cluster.

Gub as in the following equation.
p
σi ≈ Ii fˆub (xv , x, Ci )Gub (x, Ci ) Var[V ].

(13)

The variance Var[V ] for Ci is represented by the average of thevisibilities V̄i within Ci as given by the following equation.
V̄i = ∑ V (x, y)p(y),y∈Ci
2Var[V ] = V̄i − V̄i = −(V̄i − 1/2)2 + 1/4,

(14)

where p is the probability to sample
VPL as p(y) = I(y)/Ii . Ourpmethod conservatively bounds Var[V ] with the maximum value
0.5, while it can tightly estimate Var[V ] by using the visibility
caching described in Sec. 4.
Estimation of the Upper Bound of BRDF: To calculate the standard variation σi for each VPL cluster Ci , the upper bound of
BRDF within cluster is required. Some BRDF models such as Lambertian, Blinn-Phong, and isotropic Ward BRDF models have asimple and tight way to compute upper bounds as used in Lightcuts [WFA∗ 05]. This limitation is lifted in our method by using
BRDF importance sampling and exploiting VPL samples as shownin Fig. 3. BRDF importance sampling [WA09] is performed to estimate fˆub for glossy BRDFs. Based on BRDF f at shading point x,
N f rays are sampled. Then intersection tests between these rays andthe bounding box of the VPL cluster are performed. The maximum
value of BRDF f max is calculated from among the intersected rays.
As mentioned before, since our method samples two VPLs y1 , y2 ,BRDF values f (xv , x, y1 ) and f (xv , x, y2 ) are evaluated and reused
for fˆub . In addition, we calculate the minimum angle between thenormal at the shading point x and the direction ωCi that intersects
the bounding box of the cluster, which is also required in the calculation of Gub . Our method evaluates f by using ωCi as fωCi . The
maximum value of f max , f (xv , x, y1 ), f (xv , x, y2 ), and fωCi is used
for fˆub .
Additional Clustering Criterion: As mentioned before, ∆L is derived from the assumption that samples in each cluster follow a normal distribution and the Neyman allocation with uniform standard
variations is used. However, it is quite difficult to check whether
this assumption is valid or not from samples. Instead, t-distribution
is robust to violations of this assumption and Eq. (10) is valid when
the number of samples n is sufficiently large and the standard variations σi are uniform. We consider that the number of samples issufficient and the standard variations σi are uniform when the cric 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

0.5(a) n=256

0.4

0

5

0

(b) n=512

-5

0

5

Figure 4: Histograms of T for n = 256 and 512. Horizontal axis
represents statistics T and vertical axis represents the probability density. The line shows the t-distribution. Histograms of T aremeasured in the center pixel of the inset image. By increasing n, T
approaches the t-distribution.

terion ECi < εL̂ used in Lightcuts is satisfied. Since σi has the same
factors fub and Gub as in ECi as shown in Eqs. (6) and (13), ourmethod utilizes σi instead of ECi with the criterion√
2σi < εL̂ (∀Ci ∈ C),(15)√
where 2 is the scaling factor to equate σi p
and ECi , which stems
from the difference of maximum values of Var[V ] and Vub , andthe different number of samples. Fig. 4 shows the histograms of T
in Eq. (9) with different n. As shown in Fig. 4, T approaches thet-distribution by increasing n.
In summary, our method partitions VPL clusters until ∆L < εL̂
and Eq. (15) are satisfied. If these two criteria are not satisfied,
the VPL cluster Ck with the maximum standard variation σk is selected, and Ck from cut C is replaced with the left child CL and theright child CR . L̂CL , L̂CR , σL , and σR are estimated by sampling two
VPLs and reusing from Ck . Then L̂ and ∆L are updated by using
L̂CL and L̂CR .
4. Visibility Caching
4.1. Overview
Our method exploits the spatial and directional coherence of visibilities between nearby shading points and VPLs using visibility
caching. Shading points are partitioned into clusters referred to as
shading clusters. Shading clusters are initially partitioned based on
the proximities of the normals and positions of the shading points(Fig. 5(a)) in a similar way [WC13]. Several cache points are randomly sampled from the shading points in each shading cluster,and the average of the visibilities between each cache point andeach VPL cluster is estimated and recorded as visibility caches(Fig. 5(b)). The visibility caches are averaged over all the cache
points in each shading cluster (Fig. 5(c)). If the average visibility caches differ from those stored at the cache points, the shadingcluster is refined (Fig. 5(d)).
4.2. Calculating the average visibilities and refining theshading clusters
At each cache point c j , the outgoing radiance L̂ is estimated by partitioning the VPL clusters as described in Sec. 3. The average visibility V̂i between Ci and c j is estimated by

V (c j ,y1i )+V (c j ,y2i ), where2

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering

(c)

(b)

(a)

(d)

shading cluster

Figure 5: Visibility caching. (a) The average visibilities between Ci and cache point c1 are estimated and recorded at each VPL cluster
Ci . (b) The average visibilities between Ci and c2 . (c) The averaged visibilities at ancestor VPL clusters are estimated, then V1 and V2 areconstructed. In V2 , average visibilities at descendant nodes C6 and C7 are borrowed from that at C3 . V̄ is calculated by averaging V1 andV2 . δ1 and δ2 are calculated. (d) If δ1 or δ2 is larger than the threshold δVC , the shading cluster is subdivided into two clusters.

y1i and y2i are the sampled VPLs in Ci , respectively. The averagevisibility V̂a for the ancestor node Ca of cut C is estimated by
!
V̂a =

!

∑ p(y) V̂L + ∑ p(y) V̂R ,
y∈CL

(16)

y∈CR

where CL and CR are the left and right child nodes of Ca , V̂L andV̂R are the estimated average visibilities at CL and CR , and p is
the probability of sampling the VPLs. By recursively estimating
the average visibilities based on the bottom-up approach, the average visibilities at ancestor nodes of VPL clusters are estimated.
The visibility cache stores the set of estimated average visibilities for cache point c j as V j = {· · · , V̂i , · · · }. Please note that our
method estimates and records the average visibilities for the ancestor nodes of the VPL clusters only, while previous caching methods [GKPS12, YNI∗ 15] evaluate and store contributions from all
the VPLs at each cache point, resulting in expensive computationaltime and requiring an extensive memory footprint.
Then our method calculates the average V̄ of all the visibility
caches V in the shading cluster. We measure the similarity between
V j and V̄ by using the following error δ j .
N

1 j
δj =
∑ Ṽi −
N j i=1

V (c j , y1i ) +V (c j , y2i )
2

!!2
,

(17)

where N j is the number of VPL clusters in the cut C for c j , Ṽi isthe average value stored in V̄. If δ j is greater than a threshold δVC ,the caching data V̄ is not considered to be adequate for the shadingcluster and thus the shading cluster is subdivided into two clusters.
Allp
the shading points in the shading cluster use Ṽi for V̄i to estimate Var[V ] in Eq. (13). If Ṽi = 0 or Ṽi = 1, Var[V ] and σi become
zero, although VPL cluster Ci has variance in f and G. This leads to
Ci being unpartitioned and introduces errors. To address this problem, our method clamps Ṽi between 0.01 and 0.9 in order that σicannot to be zero.

5. Results
Figs. 1, 6 to 8 show the rendering results and visualize the relative
errors estimated by using our method. The computation times were
measured on a PC with an Intel Xeon E5-2697 v2 2.70 GHz CPU.
All computations were performed in parallel using multithreading.
The image resolution of Figs. 1, 7, and 8 is 1280 × 720, and thatof Fig. 6 is 10242 . In our experiments, the relative error threshold ε = 2% based on Weber’s law [WFA∗ 05], α = 95%, numberof rays N f = 256, δVC = 0.1 are used otherwise stated. In visibility caching, ten cache points are used in each shading cluster.
We have experimented with various δVC and the numbers of cache
points, but these parameters do not affect the rendering performance and the estimation accuracy so much. Our method measures
the rate Rε of pixels, which satisfy kL̂ − Lk < εL, in the image.
Each VPL cluster stores multiple representative VPLs in a similar
way [WABG06]. Our method generates VPLs in the similar way
to the traditional many-light rendering method [DKH∗ 13], and thecomputational times to generate VPLs and light trees are omitted
in the following reported computational times since those are common to Lightcuts and our method.
Fig. 1 shows the comparisons of relative error images between(b) our method (282s) and (c) Lightcuts (60s) when ε = 2% is specified. Our method achieves Rε = 92.96% while Lightcuts achieves
43.67% in (c). In order to exceed R2% = 92.96% , ε for the condition ECi < εL̂ of Lightcuts is set to 0.3% with trial and error processes. The rendering time for ε = 0.3% is 206s in Lightcuts that
estimates fˆub by using analytical solution, while our method estimates fˆub by using BRDF importance sampling. To make the comparison fair, our method employs the analytical solution to calculate fˆub , and the computational time Tc is 209s. That is, our method
can compete with Lightcuts in rendering performance, while ourmethod can directly render images with reliable accuracy without
tedious trial and error processes. We have tried to estimate the sumof the upper bounds of the clustering errors ∑N
i=1 ECi < 0.02L̂ forLightcuts, but it took more than 16 hours (58,078s) and the relative
mean error is 5.69214e-06. This indicates that Lightcuts with theestimation of total clustering errors overestimates the error.
c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering

Table 2: Computational times Tc , Tu , Rε , and MRE for α = 99%.

(a)

(b)

92

200

90

Sponza388
97.67%
0.0068549
97.68%
0.0068
1.41

Buddha103
95.93%
0.0075126
95.92%
0.0075
1.22

Kitchen89
96.67%
0.0072114
96.69%
0.0072
1.28

San Miguel
3,151
97.94%
0.00683,872
98.06%
0.0068
1.23

150

88
86

100

84

50

82
80
(c)

20%

Scene
Tc
Rε
MRE
Tu
Rε
MRE
Tu /Tc

0
64

128

256

512

1024

(d)

64

128

256

512

1024

0.0

Figure 6: Buddha scene rendered with Cook-Torrance BRDF illuminated by an area light source and an environment map. (a)rendering result, (b) relative error image, graphs of (c) Rε and (d)Tu with different N f . Error bar is shown in the bottom. According
to (c) and (d), N f = 256 is used in our experiments.

Fig. 6 shows the Buddha scene rendered with Cook-Torrance
BRDF. By using our method, Rε = 91.15% pixels satisfy the condition kL̂ − Lk < εL when α = 95% is specified. Please note that
our method bounds the error kL̂ − Lk stochastically, not deterministically. Therefore, some pixels could have higher relative errors
than the specified threshold ε. Figs. 6(c) and (d) show the graphsof Rε and computational times without visibility caching Tu withdifferent N f , respectively. As shown in these graphs, Rε grows by
increasing N f , but it grows slowly when N f exceeds 256, while thecomputational times increase considerably. Therefore, our method
uses N f = 256 in all experiments.
Table 1: Scene information of our method. N f = 256, ε = 2%,
α = 95% are used. Nvpl is the number of VPLs, and Ntri is that oftriangles. Tc and Tu are rendering times in seconds with and without visibility caching, and Tu /Tc indicates the acceleration ratiousing the visibility caching. Rε represents the rate of pixels satisfying kL̂ − Lk < εL in the image. MRE represents mean relative
error.
Scene
Fig
Nvpl
Ntri
Tc
Rε
MRE
Tu
Rε
MRE
Tu /TcSponzaFig. 1
2,061k262k
28292.96%
0.0088401
92.90%
0.0088
1.42

BuddhaFig. 62,075k
1,086k8091.15%
0.009494
91.19%
0.0094
1.18

KitchenFig. 72,082k528k
7091.86%
0.00918692.50%
0.0089
1.23

San Miguel
Fig. 8
1,936k
10,464k
2,230
92.66%
0.00892,860
92.84%
0.00881.29

c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

Figs. 7 and 8 show the relative error images with different relative error thresholds for scenes with complex geometries and
BRDFs. As shown in these images, our method can render images
whose errors from the exact solution are stochastically bounded forvarious relative error thresholds. Our method achieves that Rε exceeds 91% for all experiments when α = 95% is specified.
Table 1 shows the number of VPLs Nvpl , the computational times
Tc with visibility caching and Tu without visibility caching, acceleration ratio Tu /Tc , Rε , and the mean relative error (MRE) of each
image. As shown in Table 1, the visibility caching can speed up
the rendering in 18% to 42%, but the decreases in Rε and the mean
relative error are almost negligible. Table 2 shows the computational times Tu , Tc , the rate Rε , and the mean relative error (MRE)of each scene when α = 99% is specified. Our method achieves
that Rε exceeds about 96% for all experiments when α = 99% is
specified. The computational times Tc for α = 99% are up to 1.41
times longer than those for α = 95% in these examples.
Fig. 9 shows the histogram of relative errors in Buddha scene
shown in the inset. Our method achieves Rε = 93.13% when α =95% and ε = 2% are used. The maximum relative error in Fig. 9
is 20.5%. As shown in Fig. 9, the number of pixels whose relative
errors are larger than 2% decreases exponentially. Histograms of
the relative errors in other scenes are shown in the supplemental
material.
5.1. Discussions
The stochastic evaluation of VPL contributions can introduce noisein the rendered images. Fig. 10 shows the comparison between our
method (Figs. 10(a), (b)), Lightcuts (Fig. 10(c)), and the reference(Fig. 10(d)). Relative error threshold ε is set to 2% in our methodand Lightcuts. As shown in the zoomed-in images, Fig. 10(c) has
perceptible noises since Lightcuts cannot control the relative errorof outgoing radiance, while the noises in Fig. 10(b) are imperceptible since our method can directly control the relative error, and due
to Weber’s law, humans cannot detect the change in less than 2% inpractice.
We have applied our visibilityp
caching method to Lightcuts by
substituting Vub in Eq. (6) with 2 Var[V ]. When ε = 2% is specified, the computational time and Rε are 44s and 37.84%, while
those without visibility caching are 60s and 43.67%, respectively.
Our visibility caching method can accelerate Lightcuts at the costof decreasing Rε , while, as we mentioned before, the decrease in Rε

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering

Figure 7: Kitchen scene rendered with Cook-Torrance, Ashikhmin-Shirely BRDFs. Different relative error thresholds are used. (a) ε = 2%,
Rε = 91.86%, and Tc = 70s, (b) ε = 5%, Rε = 93.65%, and Tc = 40s, (c) ε = 10%, Rε = 95.30%, and Tc = 31s. Our method can estimate
errors well for various relative error thresholds.

Figure 8: San Miguel scene rendered with Cook-Torrance, Ashikhmin-Shirely BRDFs. Different relative error thresholds are used. (a)ε = 2%, Rε = 92.66%, and Tc = 2, 230s, (b) ε = 5%, Rε = 92.79%, and Tc = 757s, (c) ε = 10%, Rε = 92.07%, and Tc = 327s. Our method
can estimate errors well for various relative error thresholds.

0.35

error threshold ε and confidence parameter α, our method can automatically partition VPLs so that the error of the outgoing radiance is stochastically bounded. We demonstrated that our method
can render images with reliable accuracy for complex scenes andvarious materials. In addition, we introduced an efficient visibility
caching method that can accelerate rendering 18% to 42% withouta noticeable decrease in the estimation accuracy.

0.3
0.25
0.2
0.15
0.1
0.05
0
0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

Figure 9: Histogram of relative errors in Buddha scene. Horizontal
axis represents the relative errors, and vertical axis represents theprobability of occurrence. ε = 2%, Rε = 93.13%, and Tc = 86s.

of our method is almost negligible. This indicates that our visibility
caching method is well suited to our error estimation framework.
Since our method is based on many-light rendering, our method
also inherits the drawbacks of many-light rendering, such as the
singularities due to the geometry term and highly glossy BRDFs.
These singularities lead to splotches in the rendered images, but
the splotches are alleviated by increasing the number of VPLs orclamping. Due to the scalability of our method, a large number of
VPLs as shown in Table 1 can be handled very efficiently.
6. Conclusions and Future Work
We have proposed an error estimation framework for many-lightrendering by using confidence intervals. By specifying the relative

In future work, we would like to extend our method to Multidimensional Lightcuts [WABG06]. Our method can be combinedwith Multidimensional Lightcuts by estimating the errors using
clusters of the light tree and the gather tree, but an efficient method
to compute the upper bounds of phase functions and BRDFs forgather points and VPLs is required. Furthermore, we would like toestimate errors for more complex lighting effects (e.g. caustics).

Acknowledgements
This research was partially supported by JSPS KAKENHI, Grantin-Aid for Scientific Research on Innovative Areas, Grant Number
JP15H05924.

References[AUT13] AUTODESK: Autodesk 360 rendering, 2013. URL: http://rendering.360.autodesk.com/. 2[BMB15] B US N., M USTAFA N. H., B IRI V.: Illuminationcut. Computer
Graphics Forum 34, 2 (2015), 561–573. 2[CAM08] C LARBERG P., A KENINE -M OLLER T.: Exploiting visibility correlation in direct illumination. Computer Graphics Forum 27, 4(2008), 1125–1136. 2

c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

K. Nabata & K. Iwasaki & Y. Dobashi & T. Nishita / An Error Estimation Framework for Many-Light Rendering(a)(b) our(c) Lightcuts(d) reference

Figure 10: Comparison between our method, Lightcuts, and the reference. Zoomed-in images show that the noises due to (b) our stochastic
evaluation are imperceptible and indistinguishable from (d) the reference, while (c) those in Lightcuts are perceptible.

[DHKL10] DAMMERTZ H., H ANIKA J., K ELLER A., L ENSCH H.: Ahierarchical automatic stopping condition for Monte Carlo global illumination. In Proceedings of the Winter School of Computer Graphics(2010), pp. 159–164. 3[Pur87] P URGATHOFER W.: A statistical method for adaptive stochastic
sampling. Computer & Graphics 11, 2 (1987), 157–162. 3[DKH∗ 10] DAVIDOVI Č T., K ŘIVÁNEK J., H AŠAN M., S LUSALLEK P.,
BALA K.: Combining global and local virtual lights for detailed glossy
illumination. ACM Transactions on Graphics 26, 6 (2010), 143:1–143:8.
2[RGK∗ 08] R ITSCHEL T., G ROSCH T., K IM M. H., S EIDEL H.-P.,
DACHSBACHER C., K AUTZ J.: Imperfect shadow maps for efficientcomputation of indirect illumination. ACM Transactions on Graphics
27, 5 (2008), 129:1–129:8. 2[DKH∗ 13] DACHSBACHER C., K RIVANEK J., H ASAN M., A RBREE A.,
WALTER B., N OVAK J.: Scalable realistic rendering with many-lightmethods. In Eurographics 2013 - State of the Art Reports (2013). 2, 6

[SHD15] S IMON F., H ANIKA J., DACHSBACHER C.: Rich-vpls for improving the versatility of many-light methods. Computer Graphics Forum (Proceedings of Eurographics) 34, 2 (2015), 575–584. 2[GKPS12] G EORGIEV I., K ŘIVÁNEK JAROSLAV J., P OPOV S.,
S LUSALLEK P.: Importance caching for complex illumination. Computer Graphics Forum 31, 2 (2012), 701–710. 2, 6
[HJJ10] H ACHISUKA T., JAROSZ W., J ENSEN H. W.: A progressive
error estimation framework for photon density estimation. ACM Transactions on Graphics 29, 6 (2010), 144:1–144:12. 3[HKWB09] H AŠAN M., K ŘIVÁNEK J., WALTER B., BALA K.: Virtual
spherical lights for many-light rendering of glossy scenes. ACM Transactions on Graphics 28, 5 (2009), 143:1–143:6. 2[HPB07] H AŠAN M., P ELLACINI F., BALA K.: Matrix row-column
sampling for the many-light problem. ACM Transactions on Graphics
26, 3 (2007), 26:1–26:10. 2[HWJ∗ 15] H UO Y., WANG R., J IN S., L IU X., BAO H.: A matrix
sampling-and-recovery approach for many-lights rendering. ACM Transactions on Graphics 34, 6 (2015), 210:1–210:12. 2[RFS03] R IGAU J., F EIXAS M., S BERT M.: Refinement Criteria Basedon f-Divergences. In Eurographics Workshop on Rendering (2003). 3[Tho12]

T HOMPSON S. K.: Sampling, 3rd ed. Wiley, 2012. 4[TJ97] TAMSTORF R., J ENSEN H. W.: Adaptive sampling and bias estimation in path tracing. In Proceedings of the Eurographics Workshop on
Rendering Techniques ’97 (1997), pp. 285–296. 3[Tok15] T OKUYOSHI Y.: Virtual spherical gaussian lights for real-time
glossy indirect illumination. Computer Graphics Forum 34, 7 (2015),
89–98. 2[UNRD13] U LBRICH J., N OVAK J., R EHFELD H., DACHSBACHER C.:
Progressive visibility caching for fast indirect illumination. In Proc. ofInternational Workshop on Vision, Modeling, and Visualization (2013),
p. 8. 2
[WA09] WANG R., A KERLUND O.: Bidirectional importance sampling
for unstructured direct illumination. Computer Graphics Forum 28, 2(2009), 269–278. 2, 5

[Kel97] K ELLER A.: Instant radiosity. In Proc. of SIGGRAPH ’97(1997), pp. 49–56. 2[WABG06] WALTER B., A RBREE A., BALA K., G REENBERG D. P.:
Multidimensional lightcuts. ACM Transactions on Graphics 25, 3 (July
2006), 1081–1088. 2, 6, 8

[KGKC13] K ŘIVÁNEK J., G EORGIEV I., K APLANYAN A., C ANAD J.:
Recent advances in light transport simulation: Theory and practice. In
ACM SIGGRAPH 2013 Courses (2013). 2[WC13] W U Y.-T., C HUANG Y.-Y.: Visibilitycluster: Average directional visibility for many-light rendering. IEEE Transactions on Visualization and Computer Graphics 19, 9 (2013), 1566–1578. 2, 5

[MJL∗ 13] M OON B., J UN J. Y., L EE J., K IM K., H ACHISUKA T.,
YOON S.-E.: Robust Image Denoising Using a Virtual Flash Image forMonte Carlo Ray Tracing. Computer Graphics Forum (2013). 3[WFA∗ 05] WALTER B., F ERNANDEZ S., A RBREE A., BALA K.,
D ONIKIAN M., G REENBERG D. P.: Lightcuts: A scalable approach toillumination. ACM Transactions on Graphics 24, 3 (2005), 1098–1107.
1, 2, 3, 5, 6

[OP11] O U J., P ELLACINI F.: Lightslice: Matrix slice sampling for the
many-lights problem. ACM Transactions on Graphics 30, 6 (2011),
179:1–179:8. 2[PGSD13] P OPOV S., G EORGIEV I., S LUSALLEK P., DACHSBACHER
C.: Adaptive quantization visibility caching. Computer Graphics Forum32, 2 (2013). 2[PH10] P HARR M., H UMPHREYS G.: Physically Based Rendering, Second Edition: From Theory To Implementation, 2nd ed. Morgan Kaufmann Publishers Inc., 2010. 4
c 2016 The Author(s)Computer Graphics Forum c 2016 The Eurographics Association and John Wiley & Sons Ltd.

[WHY∗ 13] WANG R., H UO Y., Y UAN Y., Z HOU K., H UA W., BAO
H.: Gpu-based out-of-core many-lights rendering. ACM Transactions
on Graphics 32, 6 (2013), 210:1–210:10. 2[WKB12] WALTER B., K HUNGURN P., BALA K.: Bidirectional lightcuts. ACM Transactions on Graphics 31, 4 (2012), 59:1–59:11. 2[YNI∗ 15] YOSHIDA H., NABATA K., I WASAKI K., D OBASHI Y.,
N ISHITA T.: Adaptive importance caching for many-light rendering.
Journal of WSCG 23, 1 (2015), 65–72. 2, 6