A Survey on Recent Developments in Global Illumination Techniques
Tomoyuki NISHITA† (Fellow), Hisanari OTSU††† Hiroshima Shudo University / UEI Research, †† The University of Tokyo / UEI Research<Summary> The paper presents recent developments in global illumination techniques for producing photo-realistic images.
In pursuit of the increasing requirement for photo-realism in various industries, research and development into global
illumination techniques has received considerable attention within the graphics community from the early days. In this survey,
we review various global illumination techniques, from the early attempts to recent advances in the two major categories ofglobal illumination techniques: those based on the finite element method (FEM) and those based on the Monte Carlo method.
Recently, in the field of the Monte Carlo based global illumination techniques, some advanced techniques incorporating photon
density estimation or Markov chain Monte Carlo method is actively researched in the graphics community. We attempted topave the way for the reader to understand such advanced topics by describing historical background and motivation behindthese techniques. We hope this survey to help the readers to comprehend the overview of the field and the characteristics ofthese techniques, and to indicate the next direction of the research in the field of global illumination computation.
Keywords: survey, rendering, global illumination, light transport, finite element method, Monte Carlo method

1. Introduction
From the early days of computer graphics, rendering
photo-realistic images has been one of the most important topics inthe field. In order to achieve photo-realism, the rendering
techniques categorized as global illumination techniques, which
simulates the propagation of light energy taking into accountinter-reflections of light between surfaces, has been researched anddeveloped. Research into global illumination techniques started
from the research on the Radiosity method, which is known tooriginate from Nishita and Nakamae’s technique1) (Fig. 1). In someliterature, global illumination techniques are also known as the light

Fig. 1 The scene illuminated with the Radiosity method by Nishita andtransport techniques, since these techniques actually simulate the

Nakamae1) taking into account the diffuse inter-reflection between

transportation of light emitted from light sources to the sensor.

surfaces, which is known as an origin of the global illumination

Global illumination techniques can drastically improve thetechniques5).

appearance of the rendered images, so nowadays their application
covers various industries including the movie production industry,illumination techniques. For example, Robustness is one suchthe game production industry, and the medical industry. In order tocriterion, and it defines the ability to render a scene properly

pursue all of the possible demands from these industries andirrespective of the scene configuration. For instance, a scene withbecause of their increasing importance in the field of photo-realistic

inter-reflections between glossy surfaces or a scene containing

rendering, development and research into global illumination

complex occlusions are known to be difficult cases and even thetechniques will still be an important topic in future research of thelatest techniques often show poor performance. Rendering speed isfield of computer graphics.

another important criterion, which indicates how fast the rendering

However it is a challenging problem to develop efficient

process is completed. Accuracy and precision are also important

global illumination techniques due to the underlying complexities

criteria for numerical techniques. We note that some of theseof materials or the geometries in the rendered scenes. Various

criteria are not always independent and global illumination

criteria can be assumed in order to define the efficiencies of global

techniques often compromise conflicting criterion in some way.

The purpose of this survey is to review the various global

formulation, named the self-adjoint linear operator formulation,

illumination techniques from early attempts to recent developments,

and unifies bidirectional light transports by self-adjoint operators toand to help the readers to comprehend an overview of the field andachieve unifying support for non-symmetric BSFDs. Veach also

the characteristics of these techniques. In this survey, we mainly

proposed the path integral formulation for global illumination

focus on the global illumination effects occurred on surfaces.

computation, which simplifies the hemispherical form of the

Although the theory and techniques with participating media has

rendering equation such that measurements can be obtained frombeen actively researched, we will not review these techniques inthe direct solution rather than the solution of an integral equation.

details. Moreover unfortunately due to the lack of space and time,

Some research has focused on the theory of 2D global illumination,

we omit some types of global illumination techniques includingwhich is useful for the detailed analysis of the global illumination

pre-computation based techniques, virtual point light (VPL) based

techniques38).

techniques, real-time global illumination techniques. The types ofFrom among the above formulations, we review two major

global illumination techniques handled in the survey is summarizes

formulations: the hemispherical formulation and the path integral

in Fig. 2.

formulation. For simplicity, we define the useful terms beforehand

The organization of the paper is as follows. First in section 2,

in Table 1.

we will describe the theory of global illumination computation. In
section 3, we review the finite element method (FEM) based
techniques. In section 4, we review the Monte Carlo method based
techniques. Finally, in section 5, we conclude the survey.
SectionGlobal illumination techniques

Table 1 Mathematical notations
Notation

Meaning

𝑀

Number of pixels in the imageℳ

Union of all surfaces in the scene𝒮

2

Set of outward directions in the unit sphere𝐍(𝐱)FEM based techniques

4

Discretization

4.2

Form factor computation

4.3

Linear system computation

4.4

Monte Carlo method based techniques5

Geometry normal at surface point 𝐱′

Path sampling

5.2

Irradiance / radiance caching

5.3

Photon density estimation

5.4

Markov chain Monte Carlo (MCMC)

5.5

Fig. 2 Classification of global illumination techniques in the survey

𝐱, 𝐱 , 𝐱′′

Points on ℳ𝐱 0 𝐱1 ⋯ 𝐱 𝑛

Light path with vertices 𝐱 0 , 𝐱1 , … , 𝐱𝑛

𝑥Light path in 𝒫

𝑥𝑛

Light path in 𝒫𝑛

𝜍Solid angle measure
Projected solid angle measure𝜍𝐱⊥𝑑𝜍𝐱⊥ (𝜔) = |𝐍(𝐱) ⋅ 𝜔|𝑑𝜍(𝜔)(𝑗)𝑊𝑒 (𝐱, 𝜔),(𝑗)𝑊𝑒 (𝐱 ⟷ 𝐱′)

𝐿𝑖 (𝐱, 𝜔)

Incident radiance function at 𝐱 and direction 𝜔
Outgoing radiance function at 𝐱

𝐿𝑜 (𝐱, 𝜔)and direction 𝜔

2. Theory of Global Illumination𝐱ℳ (𝐱, 𝜔𝑜 )In this section we will briefly describe the theoretical aspectsofglobalillumination

computations.

Globalillumination

Responsively function for 𝑗-th pixel at 𝐱

𝑓𝑠 (𝐱, 𝜔𝑖 → 𝜔𝑜 )computation was initially formulated as rendering equation byFirst point in ℳ from a point 𝐱
to a direction 𝜔𝑜 (ray-casting function)
BSDF(bidirectional

Kajiya2), which is now known as the hemispherical formulation of𝑓𝑠 (𝐱 → 𝐱′ → 𝐱′′ )global illumination. The concept of the rendering equation was

𝐵(𝐱)Radiosity at 𝐱derived from the literature of the radiative heat transfer. With the𝐸(𝐱)Radiant exitance at 𝐱

development of global illumination techniques, the theory of global

𝜌(𝐱)Diffuse reflectance at 𝐱

3)

illumination has become more sophisticated. Arvo

scattering

distributionfunction)
Three-point form BSDF

precisely

formulated radiometry and rendering equation using measure
theory, and proposed the linear operator formulation of global
illumination computation, which is useful for studying the

2.1 Hemispherical formulation
Among

the

available

formulations,theso-calledproperties of global illumination computation using techniques ofhemispherical formulation is one of the most common, and various

functional analysis. Veach4) extended Arvo’s linear operator

literature have adopted this formulation as an introductory theory

for global illumination computation5), 6), 7). Due to a lack of themeasure 𝜇 . The light path space is defined as 𝒫 = ⋃∞𝑛=1 𝒫𝑛space, we will only describe a brief overview of the formulation,

where 𝒫n = *𝐱0 𝐱1 ⋯ 𝐱𝑛 | 𝐱0 , 𝐱1 , … , 𝐱𝑛 ∈ ℳ+ is the set of lightplease refer to above-mentioned literature for details.

paths with n + 1 edges (or length). The light path measure is

The purpose of global illumination computation is to computemeasurements of the light energy that contribute to each pixel. We
let the measurements for 𝑗-th pixels be 𝐼𝑗 (𝑗 = 0, … , 𝑀). Then themeasurement equation is defined as

defined as 𝜇(𝐷) = ∑∞𝑘=1 𝜇𝑘 (𝐷 ∩ 𝒫𝑛 ) where 𝐷 ⊂ 𝒫 and 𝜇𝑘 is
the area-product measure8) on 𝒫 defined as𝜇𝑘 (𝐷) = ∫𝐷 𝑑𝐴(𝐱0 ) ⋯ 𝑑𝐴(𝐱𝑛 ). Then the measurement equation
is written as

𝐼𝑗 = ∫ℳ×𝒮 2(𝑗)
𝑊𝑒 (𝐱, 𝜔)𝐿𝑖 (𝐱, 𝜔)𝑑𝐴(𝐱)𝑑𝜍𝐱⊥ (𝜔).

(1)
𝐼𝑗 = ∫ 𝑓𝑗 (𝑥)𝑑𝜇(𝑥).

defined as an integral equation with respect to the outgoingradiance function 𝐿𝑜 (𝐱, 𝜔) , which is obtained from energy
conservation between the incident and emitted and scatteredradiance at a point, assuming that we observe the equilibrium
radiance:(5)

𝒫

Then we define the light transport equation. This equation is

The function 𝑓𝑗 is the measurement contribution function andis defined for a light path with length 𝑛 as
𝑓𝑗 (𝑥𝑛 ) = 𝐿𝑒 (𝐱0 → 𝐱1 )∞⋅ ∏ 𝐺(𝐱𝑘−1 ⟷ 𝐱𝑘 )𝑓𝑠 (𝐱𝑘−1 → 𝐱𝑘 → 𝐱𝑘+1 )

(6)

𝑘=1

𝐿𝑜 (𝐱, 𝜔𝑜 ) = 𝐿𝑒 (𝐱, 𝜔𝑜 )+ ∫ 𝐿𝑜 (𝐱ℳ (𝐱, 𝜔𝑜 ), −𝜔𝑖 )𝑓𝑠 (𝐱, 𝜔𝑖 → 𝜔𝑜 )𝑑𝜍𝐱⊥ (𝜔𝑖 ).

(𝑗)⋅ 𝐺(𝐱𝑛−1 ⟷ 𝐱𝑛 )𝑊𝑒 (𝐱𝑛−1 ⟷ 𝐱𝑛 ).

(2)

𝒮2

Converting the projected solid angle measure 𝜍𝐱⊥ to the area

3. FEM Based Techniques

measure 𝐴, we obtain the three-point form of the light transport
In this section we review the major global illumination

equation using three surface points 𝐱, 𝐱 ′ , 𝐱′′ instead of angles:

techniques according to our classifications. The first section in this𝐿(𝐱′ → 𝐱′′) = 𝐿𝑒 (𝐱′ → 𝐱′′)classification considers finite element method (FEM) based+ ∫ 𝐿(𝐱 → 𝐱′)𝑓𝑠 (𝐱 → 𝐱 ′ → 𝐱′′)𝐺(𝐱 ↔ 𝐱′)𝑑𝐴(𝐱),(3)ℳtechniques, which achieves global illumination computation byusing FEM. These techniques were originally derived from thewhere 𝐺(𝐱 ↔ 𝐱′) is the geometry term obtained as a result of theRadiosity method 1), 10), 11). Through the development of the Radiosity

measure conversion (Jacobian), which is defined as

method, various techniques have been proposed. These techniques
can be separated into groups according to the focus of the execution|cos 𝜃0 cos 𝜃1 |𝐺(𝐱 ↔ 𝐱′) = 𝑉(𝐱 ↔ 𝐱′),‖𝐱 − 𝐱′‖2(4)

process used in the method: discretization, form factor computation,and linear system computation. In this survey, we adopted thesewhere 𝜃0 and 𝜃1 are the angles between the surface normal at 𝐱′

focuses as a way of classifying the sub-techniques in the Radiosity

and the outgoing vectors from 𝐱′ to 𝐱 and 𝐱′′ , respectively.

method. Although other focuses such as rendering are still𝑉(𝐱 ↔ 𝐱′) is the visibility term defined as 𝑉(𝐱 ↔ 𝐱′) = 1 if 𝐱

important, they are omitted due to a lack of space.

and 𝐱′ are mutually visible, otherwise 𝑉(𝐱 ↔ 𝐱′) = 0.
3.1 Theory of FEM based technique
2.2 Path integral formulation

In the ordinary Radiosity method, all surfaces are assumed to

The path integral formulation of global illumination
computation was proposed by Veach4). Unlike the hemisphericalformulation,

the path

integralformulationbe Lambertian diffuse surfaces. From this assumption, we obtain
the radiosity equation from Equation (3):simplifies themeasurement equation and the light transport equation to use direct

𝐵(𝐱) = 𝐸(𝐱) + 𝜌(𝐱) ∫ 𝐵(𝐱′)𝐺′(𝐱 ↔ 𝐱′)𝑑𝐴(𝐱′),integration rather than an integral equation. Some light pathℳsampling

based

globalillumination

techniques,

such

as

7)

bidirectional path tracing , can be well described by using thisformulation.

(7)

where 𝐺′(𝐱 ↔ 𝐱 ′ ) = 𝐺(𝐱 ↔ 𝐱 ′ )⁄𝜋.
In order to approximate 𝐵(𝐱) following the standard FEM
setup, the domain of 𝐵(𝐱) , that is, the surface mesh ℳ is

The measurement equation in the path integral formulation is

subdivided into small elements and some nodes are assigned on

obtained by recursively expanding Equation (3), which is defined

those elements. Then 𝐵(𝐱) is approximated by a linear sum of theas the integral on the light path space 𝒫 with the light pathbasis functions:𝑛𝑏

𝐵(𝐱) ≈ 𝐵̂(𝐱) = ∑ 𝐵𝑖 𝑁𝑖 (𝐱),

(8)

𝑖=1

where 𝑛𝑏 is the number of basis functions and 𝑁𝑖 (𝐱) is the 𝑖-thbasis function. Various selections for the basis function are possible,
and this has influence on the final result. The simplest of these is
the constant basis function and is defined as 𝑁𝑖 (𝐱) = 1 if 𝐱 is
within the element and otherwise 𝑁𝑖 (𝐱) = 0. Although many

Fig. 3 Scenes with diffuse inter-reflection by Dobashi et al.12)

techniques adopt the constant basis function, some research haswith Radiosity method with point-sampled geometry.

investigated the possibility of using non-constant basis functions 26).
By applying a weighted residual method such as the point
collocation or the Galerkin method, we obtain a way to fit 𝐵̂(𝐱) to𝐵(𝐱) by solving a linear equation. Specifically, when the constantbasis function and the Galerkin method are used, we obtain5) a formof the radiosity equation that can be seen various literatures:al.21) utilize the idea of importance transport, which is a dual
concept of the radiance concept for adaptive mesh subdivision.
More recently, Dobashi et al.12) proposed a Radiosity method forpoint-sampled geometries which does not need explicit meshsubdivision (Fig. 3).

𝑛𝑒

𝐵𝑖 = 𝐸𝑖 + 𝜌𝑖 ∑ 𝐵𝑗 𝐹𝑖𝑗 ,Another important concept related to discretization is to
(9)

𝑗=1

employ a hierarchical structure. Cohen et al.18) introduced atwo-level hierarchical mesh structure and proposed a modified

where 𝑛𝑒 is the number of elements and 𝜌𝑖 is the diffuse
reflectance associated with the 𝑖-th element, and 𝐹𝑖𝑗 is the form
factor defined asversion of the form factor computation. Hanrahan et al.19) extended
and generalized Cohen’s two-level hierarchical radiosity techniquetomulti-level.

Byintroducing

multi-level

hierarchy,

thecomputational complexity is drastically reduced to linear time.
1
𝐹𝑖𝑗 = ∫ ∫ 𝐺′(𝐱 ↔ 𝐱 ′ )𝑑𝐴 (𝐱′)𝑑𝐴(𝐱),𝐴𝑖 ℳ ℳ(10)

Smits et al.22) utilize clustering to accelerate the hierarchicalradiosity. Some methods focus on the hierarchical structure of thebasis functions used in the FEM formulation. Gortler et al20).

where 𝐴𝑖 is the area of the 𝑖-th element.
Some techniques uses the hemispherical form of the Equation

Propose wavelet radiosity which utilizes hierarchical structure for(10), which converts the area measure of the inner integral into

wavelet basis functions converted from the normal finite element

solid angle measure:basis functions, and which also allows us to construct a linear time1
𝐹𝑖𝑗 = ∫ ∫ 𝐺′′(𝐱, 𝜔)𝑑𝜍(𝜔)𝑑𝐴(𝐱),𝐴𝑖 ℳ 𝒮 2technique.
(11)

where 𝐺 ′′ (𝐱, 𝜔) = 𝜋 −1 ⋅ cos 𝜃𝑁𝐱,𝜔 ⋅ 𝑉(𝐱, 𝐱ℳ (𝐱, 𝜔)).

3.3 Form factor computation
One way to compute the form factor is to use an analyticalsolution. Goral et al.10) utilized the contour integral form of

3.2 Discretization

Equation (10) to compute the form factor between rectangular

The first classification involves in discretization of theelements without occlusions. Nishita and Nakame1) resolved theradiosity equation with FEM. We now review some techniques

restriction on occlusions by adding an extra visibility handling

related to this classification.

scheme. Baum et al.23) partially used an analytical approach toFirst we focus on the discretization of the mesh. It is known

solve the inner part of Equation (10).

that various conditions that are introduced in the discretization step

Another way is to use numerical solution techniques. Thison a mesh, such as mesh density or continuity, have a bad influence

approach was initially attempted by Cohen et al.11) and is known as

on the accuracy of the rendered images5). Various techniques havethe hemi-cube technique. This technique utilizes the hemisphericalbeen investigated in an attempt to alleviate this inaccuracy.

form of the form factor (Equation (11)) under the assumption that

15)

17)

proposed a mesh subdivisioneach element has a constant inner integral across the element,

technique for discontinuities introduces by change of illumination.

which means the form factor for an element can be evaluated as a

Lischinski et al.

and Heckbert
16)

proposed an adaptive mesh generation

single integral. The inner integral is computed numerically by

scheme utilizing a binary space partitioning (BSP) tree. Smits etprojecting scene surfaces into hemi-cube, which is defined as aCampbell and Fussellhalf-cube structure whose faces are separated into different

Although various techniques have been developed and researched,elements. Projection onto the hemi-cube can be achieved using

and classification criteria are vast, we adopted the classification

standard rasterization techniques such as the Z-buffer algorithm, orwith 4 categories: path sampling based techniques, irradianceeven with graphics hardware. Various techniques that stem fromcaching based techniques, photon density estimation based

this approach have been proposed. Max24) proposed an optimal

techniques, and MCMC (Markov chain Monte Carlo) based

sampling scheme for the hemi-cube based on some observations ontechniques. Unfortunately, due to a lack of space and time wethe arrangement of the sample directions. Some methods utilize

omitted some important categories of techniques: Pre-computation

other shapes that are different from the hemi-cube. Sillon andbased techniques (e.g. see Fig. 4) or VPL based techniques. Also,

25)

Puech

proposed an alternative technique to the hemi-cube

we will not delve into the topics like participating media (e.g. see

technique that utilizes a plane parallel to the surface of the element

Fig. 5) or adaptive sampling and reconstruction (a.k.a. filtering),with the observation that the angular contribution of the geometry

although they are very fascinating field and have been actively

term is relatively small in the area that is nearly-tangential to theresearched in the latest research. For techniques in these categories,

surface.

please refer to other literature.

Some numerical solution techniques for the form factor
computation utilize the Monte Carlo method. Sbert27), 28) makes use

4.1 Theory of Monte Carlo methodof knowledge of integral geometry (or geometric probability) to

In this section we briefly introduce some selected theoreticalcompute the form factor by using the Monte Carlo method.

concepts related to the Monte Carlo method.

3.4 Linear system computation

1)

Monte Carlo integration

In the Radiosity method, iterative methods known as

In global illumination computation with the Monte Carlorelaxation methods, such as the Jacobi, Gauss-Seidel or Southwellmethod, the central tool is the Monte Carlo integration, which is amethods are often used for linear system computation. However intechniques for numerically integrating arbitrary functions. We letorder to facilitate the potential of these techniques, the execution

the target integral we want to estimate besteps in these methods are often rearranged and interleaved. For
𝐼 = ∫ 𝑓(𝑥)𝑑𝜇(𝑥),instance, the progressive refinement radiosity method proposed by
29)

Cohen et al.

(12)

Ω

relaxes the execution flow of classical linear system
In thiswhere 𝑓 is a real-valued function defined on Ω, and 𝜇 is sometechnique, the form factor computation and the linear systemmeasure on Ω. We let a probability density function (PDF) be

computation are partially interleaved.

𝑝 = 𝑑𝑃 ⁄𝑑𝜇 on Ω with some probability space (Ω, ℱ, 𝑃) where

solution methods to achieve progressive rendering.

technique makes use

𝜇 ≪ 𝑃 ( 𝜇 is absolutely continuous to 𝑃 ) and supp(𝑓) ⊆of the Monte Carlo method to solve the linear system of thesupp(𝑝) (non-zero region of 𝑝 contains non-zero region of 𝑓).

radiosity equation. This technique alleviated the intermediate step

We consider to estimate the integral 𝐼 with 𝑁 number ofof the Jacobi iteration to converge in probability by replacing theindependent and identically distributed (i.i.d.) random variables

matrix-to-vector multiplication with an unbiased Monte Carlo

𝑋𝑖 (𝑖 = 1, … , 𝑁) generated from the distribution 𝑝. Then an

estimator, which significantly reduces the evaluation cost of

estimator 〈𝐼〉𝑁 for the integral 𝐼 can be defined as30) - 32)

The stochastic relaxation radiosity

handling large form factor matrix directly.

𝑁

The discrete random walk radiosity33) - 35) technique is another〈𝐼〉𝑁 =

way of employing the Monte Carlo method to a linear systemsolution. This technique utilizes the discrete random walk process
to estimate the solution of the linear system. In this technique,several counting strategies are considered to obtain an estimatedvalue of the solution with low variance.

1𝑓(𝑋𝑖 )∑.
𝑁𝑝(𝑋𝑖 )(13)

𝑖=1

We note that 〈𝐼〉𝑁 converges to 𝔼𝑝 ,〈𝐼〉𝑁 - = 𝐼 almost surely
as N → ∞ by the strong law of large number (SLLN), where
𝔼𝑝 ,𝑌- denotes the expected value of Y = 𝑓(𝑋) where 𝑋 ∼ 𝑝,which is defined as

4. Monte Carlo Method Based Techniques
Monte Carlo method based techniques is another major way
to achieve global illumination computation, which utilizes theMonte Carlo method to numerically solve the rendering equation.

𝔼𝑝 ,𝑌- = ∫ 𝑓(𝑥)𝑝(𝑥)𝑑𝜇(𝑥).

(14)

Ω

The variance of 𝑌 is another important quantity whichindicates how 𝑌 is dispersed around 𝔼𝑝 ,𝑌- :samples are uniformly distributed over the sample space. An

2𝑉𝑎𝑟𝑝 ,𝑌- = 𝔼𝑝 [(𝑌 − 𝔼𝑝 ,𝑌-) ]advantage of QMC over standard Monte Carlo is faster(15)2= ∫ (𝑓(𝑥) − 𝔼𝑝 ,𝑋-) 𝑝(𝑥)𝑑𝜇(𝑥).

convergence ratio. Although there is little application in global
illumination computation, variance reduction techniques such asΩcontrol variate or Rao-Blackwellization is also known and widely

In Monte Carlo based techniques, the variance is one of theused in the other fields.

most important quantity because the error convergence speed isclosely related to the variance. The asymptotic result of this relation
can be easily confirmed by the central limit theorem (CLT).
According to the theorem,

𝑁 −1⁄2 ⋅ (〈𝐼〉𝑁 − 𝐼)

is converged in

distribution to a normal distribution 𝒩(0, 𝑉𝑎𝑟𝑝 ,〈𝐼〉𝑁 -).
2)

4)

Importance sampling
Selection of the distribution 𝑝 in Equation (13) has

considerable influence on the variance of the estimator. A
techniques of improving the estimator by choosing good
distribution 𝑝 is called the importance sampling. Although it isError and bias

one of the variance reduction techniques, due its importance andWe let 𝐼 be a target quantity to be estimated and its estimatoreffectiveness in global illumination computation we separate the

be 〈𝐼〉𝑁 . Then the error of the estimator 〈𝐼〉𝑁 is defined bytechnique in an independent part. A known result for the〈𝐼〉𝑁 − 𝐼. The bias of the estimator 〈𝐼〉𝑁 is defined as the expected

importance sampling improve efficiency if the distribution 𝑝 isvalue of the error:

chosen so that 𝑝 is similar to the integrant 𝑓. More specifically,
𝐵,〈𝐼〉𝑁 - = 𝔼𝑝 ,〈𝐼〉𝑁 − 𝐼-.

(16)

The estimator 〈𝐼〉𝑁 is called unbiased if 𝐵,〈𝐼〉𝑁 - = 0 for all
𝑁, i.e., 𝔼𝑝 ,〈𝐼〉𝑁 - = 𝐼 for all 𝑁. The estimator 〈𝐼〉𝑁 is called

the optimal choice of the distribution 𝑝∗ that minimizes thevariance is known to be𝑝∗ (𝑥) =consistent if 〈𝐼〉𝑁 is converged to 𝐼 in probability as 𝑁 → ∞,that is, for all 𝜀 > 0
lim 𝑃*|〈𝐼〉𝑁 − 𝐼| > 𝜀+ = 0.

𝑁→∞

(17)

We note that unbiasedness and consistency is orthogonal
concepts, although it is often confused, so we can assume aunbiased and inconsistent estimator, for instance. In global
illumination computation, a techniques is called unbiased orconsistent according to unbiasedness or consistency of theunderlying estimator in the technique.
3)

Variance reduction
As we described above, in order to increase the efficiency of

the Monte Carlo estimator, it is important to reduce variance of theestimator. In this part of the section, we briefly introduce somevariance reduction techniques which is often used in the global
illumination computation. Stratified sampling is one of variancereduction techniques. In this technique the sample space Ω is
subdivided into 𝑚 disjoint sub-spaces Ω𝑖 (𝑖 = 1, … , 𝑚) called
stratum, and for each Ω𝑖 random sample 𝑋𝑖,𝑗 is generated.
However, stratified sampling is not effective on a high dimensional
space because the number of stratum and the number of samples
required for each stratum are increased exponentially. Latin
hypercube sampling addresses the problem by generating samples
so that only one sample is assigned for each row and each column.
Quasi-Monte Carlo (QMC) utilizes a low discrepancy sequence
instead of a sequence of random numbers for estimating an integral.
A low discrepancy sequence is deterministically chosen so that the|𝑓(𝑥)||∫Ω 𝑓(𝑥)𝑑𝜇(𝑥)|

.

(18)

We note that in global illumination computation the integrant isalways positive, so the optimal 𝑝∗ simply becomes the normalizedversion of 𝑓, and in this case the variance is always zero.
Multiple importance sampling proposed by Veach4), 40) is an
importance sampling technique which combines multiple samples
from different distributions (or techniques). These samples are
combined with weights associated with distributions. The technique
is useful for an integrant that is difficult to sample with a single
importance distribution. We let the distributions be 𝑝𝑖 (𝑖 =1, … , 𝑁) and the weighing function associated with the distribution
𝑝𝑖 be 𝑤𝑖 . Then the estimator for multiple importance sampling isdefined as
𝑁

𝑁𝑖

𝑖=1

𝑗=1

𝑓(𝑋𝑖,𝑗 )1〈𝐼〉𝑀𝐼𝑆 = ∑ ∑ 𝑤𝑖 (𝑋𝑖,𝑗 ),𝑁𝑖𝑝𝑖 (𝑋𝑖,𝑗 )(19)

where 𝑁𝑖 is a number of samples taken from the distribution 𝑝𝑖
and 𝑋𝑖,𝑗 ∼ 𝑝𝑖 . The estimator is unbiased and consistent if
following conditions on weights are satisfied: 1) for all 𝑥 ∈supp(𝑓), ∑𝑁
𝑖=1 𝑤𝑖 (𝑥) = 1 , 2) supp(𝑤𝑖 ) ⊆ supp(𝑝𝑖 ). We notethat these conditions imply that we do not need to sample wholedomain of Ω for all techniques. More specifically, for all 𝑥 in
non-zero part of 𝑓, at least one distribution generates 𝑥, that is,transition kernel for the algorithm, a.k.a., Metropolis-Hastings𝑁supp(𝑓) ⊆ ⋃ supp(𝑝𝑖 ).

(20)

𝑖=1

update, is defined as follows. Given the current sample 𝑋𝑖 , wechoose a tentative sample 𝑋′𝑖 ∼ 𝑇(𝑋𝑖 → 𝑋′𝑖 ), where 𝑇(𝑥 → 𝑦)is transition function. Then the next sample 𝑋𝑖+1 is define as

Another focus is selection of the weights. Veach proposedsome simple but effective techniques based on heuristics. One

𝑋𝑖+1 = {technique is balance heuristics:
𝑤𝑖 (𝑥) =𝑁𝑖 𝑝𝑖 (𝑥),∑𝑁
𝑘=1 𝑁𝑘 𝑝𝑘 (𝑥)𝑋𝑖 with probability 𝑎(𝑋𝑖 → 𝑋′𝑖 ),𝑋′𝑖otherwise,(24)

where 𝑎(𝑋 → 𝑋′) is acceptance probability defined as(21)𝑎(𝑋 → 𝑋′) = min (1,which is proved to have good variance bound (Theorem 9.2 of𝑓(𝑋′)𝑇(𝑋 ′ → 𝑋)).
𝑓(𝑋)𝑇(𝑋 → 𝑋′)(25)

Veach4)).

We note that Equation (25) implies 𝑓 need not to be normalized,

5)

and we can observe the transition kernel associated with theMarkov chain Monte Carlo (MCMC)
Although importance sampling is effective for manyalgorithm:

configuration, sampling from the general path space is still𝐾(𝑥 → 𝑦) = 𝑎(𝑥 → 𝑦)𝑇(𝑥 → 𝑦)challenging. In order to resolve the problem, some recent global+ (1 − ∫ 𝑎(𝑥 → 𝑦)𝑇(𝑥 → 𝑦)𝑑𝜇(𝑦)) 𝛿𝑥 (𝑦)illumination techniques begin to use Markov chain Monte Carlo(MCMC) method. In this part of the section, we briefly review thetheoretical aspect of MCMC that is useful for understanding recent
MCMC based global illumination techniques.

(26)

Ω

preserves Equation (23) with 𝑓 as the stationary distribution.
As for the convergence property of a Markov chain, we canshow that if a Markov chain satisfied the following two conditions:

The basic idea of MCMC is to use a correlated sequence of1) 𝑓-irreducible which means that any part of the state space canrandom variables called a Markov chain, which is a sequence ofbe reached from the selection of the initial state, and 2) apediodic

random variable whose samples depend only on samples one

which means the Markov chain does not contain cycles, a sequence

before. Specifically, a Markov chain is defined as a sequence ofof samples is converged to follow the stationary distribution. And

random variables 𝑋1 , 𝑋2 , ⋯ such that 𝑋𝑖+1 depends only on 𝑋𝑖(𝑖 = 1, 2, ⋯), that is, 𝑋𝑖+1 ∼ 𝐾(𝑋𝑖 → 𝑋𝑖+1 ) where 𝐾(𝑥 → 𝑦) is

similar to SLLN for i.i.d. random number case (see Equation (13)),the expected value of the estimator 〈𝐼〉𝑁 converges to 𝐼 evena conditional probability density function called the transition

when a Markov chain satisfying above conditions is used as akernel. The algorithm which generates stationary Markov chain

sequence of samples, which is known as Ergodic Theorem. We also

according to some distribution is called Markov chain Monte Carlonote that MH update satisfies these conditions with a loose(MCMC). A Markov chain is called stationary if there exists a

conditions with 𝑓 and the transition function75).

distribution π such that

Convergenceproperty of Markov chain is important topic but we will not handle

𝜋(𝑥) = ∫ 𝐾(𝑦 → 𝑥)𝜋(𝑦)𝑑𝜇(𝑦).

more it in the survey, for details please refer to Meyn and Tweedy74),(22)ΩIn order to obtain a stationary Markov chain, many MCMC
techniques utilize a stronger condition called detailed balance
condition. A transition kernel 𝐾 satisfied the detailed balance
condition if there exists a distribution 𝜋 for all 𝑥, 𝑦 ∈ Ω,
𝐾(𝑥 → 𝑦)𝜋(𝑥) = 𝐾(𝑦 → 𝑥)𝜋(𝑦).

Robert and Casella75) etc.
4.2 Local sampling based techniques
In this section, we introduce local sampling based techniques,
which focus on solving the rendering equation by Monte Carlo
method with sampling light paths with a sequence of local(23)

sampling on scene surfaces. Here we call the techniques in thisA Markov chain which satisfies this condition is called a reversible

section are local because all sampling processes are based on theMarkov chain, and it can be shown that in this case the distribution

distributions that is defined on a surface point. It implies a light

π is stationary.

path is generated by sequentially sampling next surface points

In order to numerically compute Equation (12) by using

according to some distribution associated with the points.

Monte Carlo integration with MCMC, we want to generate samples

Specifically, we focus on the global illumination techniques withaccording to some distribution proportional to 𝑓 . In order toMonte Carlo method in which the techniques sample the light paths

achieve the goal, now we review one of the most famouswith i.i.d. random variables. This type of the techniques does nottechniques: the Metropolis-Hastings (MH) algorithm. The

facilitate the correlation of samples like techniques described insection 4.5, and many of them are statistically unbiased, so analysis
and implementation of the techniques is relatively simple.
An idea of sampling light path is originated from theWhitted’s ray tracing36) which handles global illumination
restricted to perfect specular surfaces. While the techniques itself isdeterministic, he mentioned the possibility of using random
sampling to simulate specular or glossy surfaces. Initial attempt toutilize the Monte Carlo method for global illumination computation

Fig. 4 Scenes with precomputed radiance transfer (PRT) taking

is the distributed ray tracing proposed by Cook et al.37) which caninto account inter-reflection of light by Iwasaki et al12).

handle graphics effects such as depth-of-filed or motion-blur in theframework of the ray tracing. As we describes in section 2, firstfrequency. For those scene with diffuse surfaces, as we described incomprehensive framework for global illumination computation issection 3.1, the rendering equation can be rewritten with irradiance,developed by Kajiya2), whose technique is now known as the path

which means the transported quantity can be represented as

tracing. Similar to the ray tracing, the path tracing samples a lightirradiance. Irradiance values are associated with positions onpath by tracing a ray from a camera, and when the ray intersects

surface. Although the technique is biased due to the interpolation ofwith the scene surface, next ray direction is sampled from theirradiance values, using the technique, a rendered image without

distribution associated with the BSDF on the intersected point, andhigh frequency noise can be obtained in relatively small running

the next ray is propagated to the direction until the ray hits to a lighttime. Execution steps of this techniques is as follows. If a ray

source. Along with the propagation, the radiance carried along theemitted from a camera is intersected with scene surface, thesampled path is computed according to the equations describes intechnique collects cached irradiance values from the octree storage.

section 2.1. Following Kajiya’s path tracing, various techniques

Only the values which can be well interpolated at the position using

have been proposed. As an reverse approach of the path tracing,a metric with the error estimation based on the split-sphereArvo and Kirk39) proposed the inverse path tracing (a.k.a. the lightheuristics are collected. Reciprocals of estimated error values are

tracing, or the particle tracing). The bidirectional path tracing41), 42)also used as weights for interpolation. If at least a cached value iscombines both light paths traced from a camera and a light source,collected, the interpolation using the samples are dispatched.

which can efficiently generate a light path which is difficult toOtherwise if no cached value is found, a new irradiance value issample by path tracing or light tracing, e.g., sampling from thecomputed with ordinary techniques like path tracing, and computed

scene that a camera and a light source are placed in different rooms

irradiance value is stored into the octree storage.

connected by a small window. According to the combination of thevertex connection between sampled light paths, some light paths

2)

Radiance cachingwith same lengths sampled from different densities are obtained.

Some recent techniques attempt to alleviate the restriction that

These light paths are combined by the multiple importancethe irradiance caching is only applicable to scenes with diffuse

sampling (MIS) described in section 4.1.

surfaces, and to support the material like glossy surfaces. Krivanek
et al.44) proposed a technique called the radiance caching which

4.3 Irradiance / radiance caching based techniques

extends the irradiance caching to support glossy surfaces by

From this section to section 5.4, we will introduce thetechniques based on caching, which stores intermediate data intosome data structure and utilize them for rendering. In this section
we introduce the techniques based of the caching of irradiance orradiance, which is the radiometric quantity handled in the light

representation. In the octree storage, radiance is recorded instead ofirradiance. Extending their technique, Scherzer et al.46) utilizes
improved the efficiency of the radiance caching by alleviating theevaluation cost of SHs by introducing pre-filtered MIP maps. Sometechniques attempt to support glossy surfaces without SHs. The

transport.
1)

precomputing all scene BSDFs as spherical harmonics (SH)radiance-cache splatting proposed by Gautron et al.45) improved

Basicsthe rendering part of the irradiance / radiance caching by
43)

The irradiance caching proposed by Ward et al. is the firsteliminating nearest-neighbor queries to the data structure. In thistechnique of this kind, and the techniques is applicable to scenes

technique, all stored values are splatted to the image plane andwith diffuse surfaces like classical radiosity method. A key

interpolation is computed with the information. The technique canobservation is that on diffuse surfaces the indirect illumination due

be implemented in GPUs except for irradiance / radiance cacheto the diffuse inter-reflection has relatively low variation orrecording parts. The anisotropic radiance-cache splatting proposed

by Herzog et al.46) extended the idea to utilize the orientation of

proposed a technique with a new structure that stores photon ray

splatted cache items in order to reduce errors along illumination

paths instead of photons, which is effective for eliminating bias on

gradient.

density estimation due to topology of the geometry. Herzog et al.54)

3)extended the approach to utilize splatting, projection onto the image

Improving error estimation

plane. Chen et al.55) utilizes clustering of photons according to theAn important focus on the irradiance / radiance cachingintersection history, which is used for generating the polygonaltechniques is how to design the error estimation scheme. Ward etboundaries for expressing complex caustics. Spencer and Jones58)

48)

al.

proposed a new interpolation scheme utilizing the irradiancegradient, which improves the accuracy of error estimation utilizing

improve the final gathering by introducing hierarchical structure to
the photon map.

first order gradient of irradiance. In this technique, additional
information for gradient computation is stored in the octree storage

3)

Improving photon distribution

as well as irradiance values, and later in the interpolation step these

Another focus on improving the technique is improvement of

information is used for gradient computation. Tabellion andphoton distribution in the photon map. Some techniques attempt to49)

proposed an interpolation technique utilizing theachieve this goal by refinement of the photon map withcoherency of the caching which is effective for the surfaces withredistribution of photons after photon tracing step. Spenser et al.56)

highly geometric details such as bump-mapped surfaces. Moreproposed a technique called the photon relaxation that iteratively

recently, Schwarzhaupt et al.50) proposed sophisticated error metric

redistributes photons according to the blue noise distribution. They

utilizing Hessians (a.k.a. second order gradients). This techniquelater extended57) the technique is to support complex illumination

utilizes extended version of an analysis of irradiance Hessian in 2D

structure containing discontinuity by storing photons with their

Lamorlette38)

global illumination . They also utilizes occlusion information in

trajectory information into high-dimensional kd-tree. Anotherorder to improves the accuracy of the Hessian computation.

approach to improve photon distribution is to improve the photon
tracing step itself. Peter and Pietrek62) improved photon distribution

4.4 Photon density estimation based techniquesin photon tracing step by introducing the importance map. Suykens

Next, we review second caching based techniques: photon

and Willems61) controlled photon density by a criterion based on

density estimation based techniques. This type of techniques isimportance. Hachisuka and Jensen59) proposed a photon tracing

originated from the photon mapping proposed by Jensen51).

technique utilizing Markov chain Monte Carlo (MCMC) sampling1)

according to the visibility of photon paths.

BasicsThe photon mapping is executed by two steps: 1) photon

4)

Improving density estimation kernels

tracing step and 2) photon density estimation step. In the first step,

Another way to improve the photon mapping is to selectlight particles a.k.a. photons are traced from the light sources. In theappropriate density estimation kernels. An obvious idea for the

propagation of a photon, intermediate radiance value is computed.

kernel selection problem is to modify the bandwidth of kernels. The

When a photon is intersected with scene surface, the photon is

idea of bandwidth selection has been employed in the context of thestored into the kd-tree storage called the photon map. In the second

photon density estimation. Jensen and Christensen60) initiallystep, the image is reconstructed by ray tracing with photon density

attempt to control the bandwidth by the technique called differential

estimation. When a ray is intersected with diffuse surface throughchecking, which limits the bandwidth by checking the difference inpropagation, the photons near the intersected point is gathered byirradiance. Schregle63) proposed a bandwidth selection techniquenearest-neighbor query to the photon map and used for the density

based on the bias estimation of the reconstructed illumination withestimation. Various techniques originated from the photon mapping

binary search. Schjøth et al.64), 65) proposed a technique based on

is proposed. Havran et al.52) proposed the reverse photon mapping

diffusion filtering, which is an density estimation technique withwhich utilizes the ray tracing for the initial step, and the photon

anisotropic kernel. Hey and Purgathor67) proposed a techniquetracing in the second step. A new kd-tree structure called reverse

utilizing surface geometry information, which is effective for edge

photon map is introduced for preserving information of the rayor corners. More recently, Kaplanyan and Dachsbacher66) proposed

tracing step, and in the second step nearest-neighbor query is issued

an technique with adaptive bandwidth selection for progressiveon the reverse photon map.

photon mapping with plug-in bandwith selection.

2)

5)

Improving structure

Progressive estimation

Some techniques introduces a new structure as a replacement

Although Jensen’s photon mapping is consistent in the sense53)

that if we can take infinite number of samples in the photon map.

of the photon map to improve the efficiency. Havran et al.

However it is virtually impossible due to memory restriction. In
order to resolve the problem and to develop a consistent techniquewith realistic configuration, Hachisuka et al.68) proposed atechnique called progressive photon mapping, which progressively
diminishes kernel radius while increasing number of photons to be
used for the estimation. The technique is extended69) to support
some rendering effects such as depth-of-field or motion blur. Later
Hachisuka et al.71) developed an error estimation framework for thetechnique which can be utilized as a criterion for termination ofrendering. Knaus and Zwicker70) proposed a variant of the
progressive photon mapping which is not dependent on local
statistics and can be utilize arbitrary kernel with a probabilistic
analysis of the technique. More recently, Georgiev et al.72) andHachisuka et al.73) independently proposed a generalized
framework for the progressive photon mapping, which combines
photon density estimation and multiple importance sampling inbidirectional path tracing.

Fig. 5 Scenes with participating media by Yue et al.14). In their technique,free path sampling for inhomogeneous participating media is

4.5 MCMC based techniques

improved by automatic space partitioning scheme.

In this section, we review the techniques based on Markov
chain Monte Carlo (MCMC) method described in section 4.1.

bidirectional path tracing. A sample is assigned one of the space1)and can move beyond the spaces, similarly to serial tempering

Basics
MCMC is initially brought to the global illumination

computation by Metropolis light transport (MLT) proposed by
Veach and Guibas77). MLT utilizes the Metropolis-Hastings
algorithm for sample generation and introduces some light path
specific mutation strategies such as the bidirectional mutation or
perturbations, which correspond to transition function in section
4.1. In these strategies, light paths are modified directly in the lightpath space, for example by slightly changing the next ray directions
or light path lengths, which enables the local exploration in the lightpath space.

which is one of MCMC technique with tempered distributions.
Manifold exploration proposed by Jakob and Marschner 78), 82)
introduces the specular manifold, which is the space defined in thelocal coordinates of the tangent spaces for each intersection points
with the constraints of reflection or refraction vectors for specularsurfaces. In the technique, a new mutation strategy is introduced to
explore on the space with a iterative equation solving technique like
Newton's method. Very recently, the technique is extended andimproved by Kaplanyan et al.83) using a space of the projected half
vectorrepresentation

of

light

paths

with

constraints.
79)

Gradient-domain MLT proposed by Lehtinen et al.

introduced

Changing sample space

image space gradient distribution of the energy contribution

Some variants of MLT have been developed for changing thefunction of the light paths. Instead of directly sampling from energy

sampling space and the target distribution. Primary sample space

distribution, the technique samples a light path in the gradient space.

MLT proposed by Kelemen et al.78) introduces the primary sample

Metropolis-Hastings sampler driven by distribution on the gradientspace, which is a space of high-dimensional unit cube of uniform

space enables for the samples to concentrate the on the portion of

random numbers used for sampling a light path, andthe light path space that energy is changed drastically, e.g., edges on

Metropolis-Hastings algorithm is dispatched in the space. This

the mesh. The final rendered image is reconstructed by Poisson

technique simplifies the mutation strategy and reduces high energy

solver from the gradient image.

region of the light path space thank to the importance sampling.

3)2)

More recently, Hachisuka et al.84) proposed a technique called themultiplexed MLT on Kelemen’s framework. In this technique, aprimary sample space is extended to multiple spaces, and in eachspace corresponding light paths are sampled to follow the
distributions associated with the sampling techniques used in the

Utilizing advanced MCMC techniques
Some technique attempts to use advanced MCMC techniques

other than MH algorithm. Energy preserved path tracing (ERPT)
proposed by Cline et al.80) alleviates the detailed balance andergodicity condition to achieve the techniques which depends only
on Veach's perturbations. Base on the path sampled from theordinary light path sample techniques, some short burns of the

Kaufmann (2010).

perturbations are dispatched followed by the scatter of the energy8) E. Veach and L. J. Guibas, “Bidirectional estimators for light transport”,

over the image plane without introducing any bias. Replica

Proc. of Eurographics Rendering Workshop 1994, pp.147-162 (1995).

85)

exchange light transport proposed by Kitaoka et al.

applies9) P. R. Halmos, Measure Theory, Springer-Verlag (1974).

replica exchange method to the global illumination computation. In10) C. M. Goral, K. E. Torrance, D. P. Greenberg, and B. Battaile,

the technique, a set of paths divided by the feature of the paths such

“Modelling the interaction of light between diffuse surfaces”, Computer

as random paths or caustic paths, which preserves coherency of theGraphics (Proc. of SIGGRAPH 84), Vol.18, No.3, pp.212-222 (1984).

81)

proposed a11) M. Cohen and D. P. Greenberg, “The hemi-cube: A radiosity solution

technique using the population Monte Carlo (PMC) method, which

for complex environments”, Computer Graphics (Proc. of SIGGRAPH

is a variant of MCMC technique. Unlike ordinary MCMC85), Vol.19, No.3, pp.31-40 (1985).

samples in the same type of the paths. Lai et al.

techniques like MH algorithm PMC utilizes a number of samples ina batch (a.k.a. population), which enables to exploit information

12) Y. Dobashi, T. Yamamoto, and T. Nishita, “Radiosity for point-sampledgeometry”, Proc. of Pacific Graphics 2004, pp.152-159 (2004).

collected from the samples. Their technique utilize D-kernel PMC,

13) K. Iwasaki, Y. Dobashi, F. Yoshimoto, and T. Nishita, “Precomputed

which adapts kernel function with populations per iteration. They

radiance transfer for dynamic scenes taking into account light

apply

interreflection”, Proc. of the 18th Eurographics conference on

D-kernel

PMC

to

ERPT

and

iteratively

adaptsenergy-redistributed region.

Rendering Techniques, pp.35-44 (2007).
14) Y. Yue, K. Iwasaki, Bing-Yu Chen, Y. Dobashi, and T. Nishita,“Unbiased,

5. Conclusion
In this survey, we review various types of global illumination
techniques as well as underlying theoretical concepts. As we review
in the previous sections, there are various technique and topic on
the global illumination computation and further researches are
required to solve remaining difficulties. In that means, the researchon the global illumination computation will certainly remain to be
one of the central topics of the computer graphics. We hope thissurvey to help the readers to comprehend the overview of the field
and the characteristics of these techniques, and to indicate the next
direction of the research in the field of global illumination

Adaptive

Stochastic

Sampling

forRendering

Inhomogeneous Participating Media”, ACM Transactions on Graphics(Proc. of SIGGRAPH Asia 2010), Vol.29, No.5, Article 177 (2010).
15) D. Lischinski, F. Tempieri, and D. P. Greenberg, “Discontinuity
meshing for accurate radiosity”, IEEE Computer Graphics andApplications, Vol.12, No.6, pp.25-39 (1992).
16) A. Campbell and D. S. Fussell, “Adaptive mesh generation for global
diffuse illumination”, Computer Graphics (Proc. of SIGGRAPH 90),
Vol.24, No.4, pp.155-164 (1990).
17) P. Heckbert, “Discontinuity meshing for radiosity”, In Third
Eurographics Workshop on Rendering, pp.203-226 (1992).
18) M. Cohen, D. P. Greenberg, D. S. Immel, and P. J. Brock, “An efficientradiosity approach for realistic image synthesis”, IEEE Computer

computation.

Graphics and Applications, Vol.6, No.3, pp.26-35 (1986).
19) P. Hanrahan, D. Salzman, and L. Aupperle, “A rapid hierarchical

References
1) T. Nishita and E. Nakamae, “Continuous tone representation of three
dimensional objects taking account of shadow and interreflection”,
Computer Graphics (Proc. of SIGGRAPH 85), pp.23-30 (1985).
2) J. T. Kajiya, “The Rendering Equation”, Proc. of SIGGRAPH 86,
pp.143-150 (1986).
3) J. Arbo, “Analytic Methods for Simulated Light Transport”, PhD Thesis,
Yale University (1995).
4) E. Veach, “Robust Monte Carlo methods for Light Transport
Simulation”, PhD Thesis, Stanford University (1997).
5) M. F. Cohen and J. R. Wallace, Radiosity and Realistic Image Synthesis,
Academic Press (1993).
6) P. Dutre, K. Bala, and P. Bekaert, Advanced Global Illumination, A.K.
Peters (2006).
7) M. Pharr and G. Humphreys, Physically Based Rendering, Morgan

radiosity algorithm”, Computer Graphics (Proc. of SIGGRAPH 91),Vol.25, No.4, pp.197-206 (1991).
20) S. J. Gortler, P. Schroder, M. F. Cohen, and P. M. Hanrahan, “Wavelet
radiosity”, Computer Graphics (Proc. of SIGGRAPH 93), pp.221-230(1993).
21) B. E. Smits, J. R. Arvo, and D. H. Salesin, “An importance driven
radiosity algorithm”, Computer Graphics (Proc. of SIGGRAPH 92),
Vol.26, No.4, pp. 273-282 (1992).
22) B. Smits, J. Arvo, D. Greenberg, “A clustering algorithm for radiosity incomplex environment”, Proc. of SIGGRAPH 94, pp.435-442 (1994).
23) D. R. Baum, J. R. Rushmeier, and J. M. Winget, “Improving radiosity
solutions through the use of analytically determined form factors”,
Computer Graphics (Proc. of SIGGRAPH 89), Vol.23, No.3
pp.325-334 (1989).
24) N. Max, “Optimal sampling for Hemicubes”, IEEE Transactions on
Visualization and Computer Graphics, Vol.1, No.1, pp.60-76 (1995).

25) F. Sillon and C. Puech, “A general two-pass method integrating
specular and diffuse reflection”, Computer Graphics (Proc. of41) E. Veach and L. Guibas, “Bidirectional estimator for light transport”,
Computer Graphics (Proc. of SIGGRAPH 95), pp.419-428 (1995).
42) E. Lafortune and Y. D. Willems, “Bi-directional path-tracing”,

SIGGRAPH 89), Vol.23, No.3, pp.335-344 (1989).
26) F. Tampieri and D. Lischinski, “The constant radiosity assumption
syndrome”, Proc. of the Second Eurographics Workshop on Rendering,
pp.83-92 (1991).

Compugraphics 93, pp.145-153 (1993).
43) G. J. Ward, F. M. Rubinstein, and R. D. Clear, “A ray tracing solution
for diffuse interreflection”, Computer Graphics (Proc. of SIGGRAPH27) M. Sbert, “An integral geometry based method for fast form-factor
computation”, Computer Graphics Forum, Vol.12, No.3, pp.409-420(1993).

88), Vol.22, No.4, pp.85-92 (1988).
44) J. Krivanek, P. Gautron, S. Pattanaik, K. Bouatouch, “Radiance cachingfor efficient global illumination computation”, IEEE Trans. On28) M. Sbert, “The use of global random directions to compute radiosity –Visualization and Computer Graphics, Vol.11, No.5, pp.550-561 (2005).

global monte carlo techniques”, PhD Thesis, Universitat Politecnica de45) P. Gautron, J. Krivanek, K. Bouatouch, and S. Pattanaik, “Radiance

Catalunya (1996).

cache splatting: a GPU-friendly global illumination algorithm”, Proc. of29) M. Cohen, S. E. Chen, J. R. Wallace, and D. P. Greenberg, “AEurographics Symposium on Rendering 2005 (2005).

progressive refinement approach to fast radiosity image generation”,

46) R. Herzog, K. Myszkowski, and H.-P. Seidel, “Anisotropic

Computer Graphics (Proc. of SIGGRAPH 88), Vol.22, No.4, pp.75-84

radiance-cache splatting for efficiently computing high-quality global(1988).

illumination with lightcuts”, Computer Graphics Forum, Vol.28, No.2,

30) L. Neumann, M. Feda, M. Kopp, and W. Purgathofer, “A new stochastic
radiosity method for highly complex scenes”, Proc. of Fifth

31) L. Neumann, W. Purgathofer, R. Tobler, A. Neumann, P. Elias, M. Feda,and X. Pueyo, “The stochastic ray method for radiosity”, Rendering
Techniques 95 (Proc. of Sixth Eurographics Workshop on Rendering),

Eurographics Symposium on Rendering 2012), Vol.4, No.31 (2012).
48) G. Ward, and P. Heckbert, “Irradiance gradients”, Proc. of Eurographics
Workshop Rendering (1992).
49) E. Tabellion, and A. Lamorlette, “An approximate global illumination

pp. 206-218 (1995).
32) A. Neumann, L. Neumann, Ph. Bekaert, Y. D. Willems, and W.
“Importance-driven47) D. Scherzer, C. H. Nguyen, T. Ritschel, and H.-P. Seidel,“Pre-convolved radiance caching”, Computer Graphics Forum (Proc. ofEurographics Workshop on Rendering, pp.195-206 (1994).

Purgathofer,

pp.259-268 (2009).

stochasticrayradiosity”,

In

Eurographics Rendering Workshop 1996, pp.111-122 (1996).
33) S. Pattanaik and S. Mudur, “Adjoint Equations and Random Walks forIllumination Computation,” ACM Trans. Graphics, Vol.14, No.1,
pp.77-102 (1995).

system for computer generated films”, ACM Trans. Graphics, Vol.23,
No.3, pp.469-476 (2004).
50) J. Schwarzhaupt, H. W. Jensen, and W. Jarosz, “Practical Hessian-based
error control for irradiance caching”, ACM Trans. Graphics (Proc. of
ACM SIGGRAPH Asia 2012), Vol.32, No.6, Article 193 (2012).
51) H. W. Jensen, “Global illumination using photon maps”, Proc. of the34) M. Sbert, “Error and complexity of random walk monte carlo radiosity”,

Eurographics Workshop on Rendering Techniques 96, pp.21-30 (1996).

IEEE Transactions on Visualization and Computer Graphics, Vol.3,

52) V. Havran, R. Herzog, and H.-P. Seidel, “Fast final gathering via reverse

No.1, pp. 23-38 (1997).
35) M. Sbert, “Optimal source selection in shooting random walk montecarlo radiosity”, Computer Graphics Forum (Proc. of Eurographics 97),
pp.301-308 (1997).
36) T. Whitted, “An improved illumination model for shading display”,
Communication of the ACM, Vol.32, No.6, pp.343-349 (1980).
37) R. L. Cook, T. Porter, and L. Carpenter, “Distributed ray tracing”,
Computer Graphics (Proc. of SIGGRAPH 84), Vol.18, No.3,
pp.137-145 (1984).

photon mapping”, Computer Graphics Forum (Proc. of Eurographics
2005), Vol.24, No.3, pp.323-333 (2005).
53) V. Havran, J. Bittner, R. Herzog, and H.-P. Seidel, “Ray maps for global
illumintion”, Eurographics Symposium on Rendering, pp.43-54 (2005).
54) R. Herzog, V. Havran, S. Kinuwaki, K. Myszkowski, and H.-P. Seidel,“Global illumination using photon ray splatting”, Computer Graphics
Forum, Vol.26, No.3, pp.503-513 (2007).
55) L.-H. Chen., T.-C. Tsai, Y.-S. Chen, “Grouped photon mapping”, The
Visual Computer, Vol.26, No.3, pp.217-226 (2010).

38) W. Jarosz, V. Schonefeld, L. Kobbelt, and H. W. Jensen, “Theory,56) B. Spenser and M. W. Jones, “Into the blue: Better caustics throughanalysis and application of 2D global illumination”, ACM Transactions

photon relaxation”, Computer Graphics Forum (Proc. of Eurographics

on Graphics, Vol.31, No.5, Article 125 (2012).

2009), Vol.28, No.2, pp.319-328 (2009).

39) J. Arvo and D. Kirk, “Particle transport and image synthesis”,

57) B. Spenser and M. W. Jones, “Photon parameterization for robust

Computer Graphics (Proc. of SIGGRAPH 90), Vol.24, No.4, pp.63-66

relaxation constraints”, Computer Graphics Forum (Proc. of(1990).

Eurographics 2013), Vol.32, No.2, pp.83-92 (2013).

40) E. Veach and L. Guibas, “Optimally combining sampling techniques for58) B. Spenser and M. W. Jones, “Hierachical photon mapping”, IEEE

Monte Carlo rendering”, Rendering Techniques 94, pp.147-162 (1994).

Transactions on Visualization and Computer Graphics, Vol.15, No.1,

pp.49-61 (2009).
59) T. Hachisuka and H. W. Jensen, “Robust adaptive photon tracing using
photon path visibility”, ACM Transactions on Graphics, Vol.30, No.5
Article 114 (2011).

75) C. P. Robert and G. Casella, Monte Carlo statistical methods,
Springer-Verlag (2004).
76) E. Veach and L. J. Guibas, “Metropolis light transport”, Proc. of
SIGGRAPH 97, pp.65-76 (1997).

60) H. W. Jensen and N. J. Christensen, “Photon maps in bidirectional77) C. Kelemen, L. Szirmay-Kalos, G. Antal, and F. Csonka, “A simple andmonte carlo ray tracing of complex objects”, Computers and Graphics,

robust mutation strategy for the metropolis light transport algorithm”,

Vol.19, No.2, pp.215-224 (1995).

Computer Graphics Forum, 21 (2002).

61) F. Suykens and Y. D. Willems, “Density control for photon maps”, Proc.

78) W. Jakob and S. Marschner, “Manifold exploration: a Markov chainof the Eurographics Workshop on Rendering Techniques, pp.23-34

Monte Carlo technique for rendering scenes with difficult specular(2000).

transport”, ACM Transactions on Graphics (Proc. of SIGGRAPH 2012),

62) I. Peter and G. Pietrek, “Importance driven construction of photonmaps”, Rendering Techniques (Proc. of Eurographics Rendering
Workshop 98), pp.269-280 (1998).
63) R. Schregle, “Bias compensation for photon maps”, Computer Graphics
Forum, Vol.22, No.4, pp.729-742 (2003).
64) L. Schjøth, O. F. Olsen, and J. Sporring, “Diffusion based photon
mapping”, Computer Graphics Forum, Vol.27, No.8, pp.2114-2127(2008).
65) L. Schjøth, “Anisotropic density estimation in global illumination”,
PhD Thesis, University of Copenhagen (2009).
66) A. S. Kaplanyan, C. Dachsbacher, “Adaptive progressive photon
mapping”, ACM Transaction on Graphics, Vol.32, No.2 Article 16(2013).
67) H. Hey and W. Purgathofer, “Advanced radiance estimation of photon
map global illumination”, Computer Graphics Forum, Vol.21, No.3,
pp.541-545 (2002).
68) T. Hachisuka, S. Ogaki, and H. W. Jensen, “Progressive photon
mapping”, ACM Transactions on Graphics (Proc. of SIGGRAPH Asia
2008), Vol.27, No.5 Article 130 (2008).
69) T. Hachisuka and H. W. Jensen, “Stochastic progressive photon
mapping”, ACM Transactions on Graphics (Proc. of SIGGRAPH Asia2009), Vol.28, No.5 Article 141 (2009).
70) C. Knaus and M. Zwicker, “Progressive photon mapping: A
probabilistic approach”, ACM Transactions on Graphics, Vol.20, No.3
Article 25 (2011).
71) T. Hachisuka, W. Jarosz, and H. W. Jensen, “A progressive error
estimation framework for photon density estimation”, ACM
Transactions on Graphics (Proc. of ACM SIGGRAPH Asia 2010),
Vol.29, No.6 Article 144 (2010).
72) I. Georgiev, I. Krivanek, T. Davidovic, and P. Slusallek, “Light transport
simulation with vertex connection and merging”, ACM Transactions on
Graphics (Proc. of SIGGRAPH Asia 2012), Vol.31, No.6 Article 192(2012).
73) T. Hachisuka, J. Pantaleoni, and H. W. Jensen, “A path space extension
for robust light transport simulation”, ACM Transactions on Graphics(Proc. of SIGGRAPH Asia 2012), Vol.31, No.6 Article 191 (2012).
74) S. Meyn and R. Tweedie, Markov chains and stochastic stability,
Springer-Verlag (1993).

Vol.31, No.4 Article 58 (2012).
79) J. Lehtinen, T. Karras, S. Laine, M. Aittala, F. Durand, and T. Aila,“Gradient-domain metropolis light transport”, ACM Transactions on
Graphics (Proc. of SIGGRAPH 2013), Vol.32, No.4 Article 95 (2013).
80) D. Cline, J. Talbot, and P. Egbert, “Energy redistribution path tracing”,
Proc. of SIGGRAPH 2005, pp.1186-1195 (2005).
81) Y.-C. Lai, S. Fan, S. Chenney, and C. Dyer, “Photorealistic image
rendering with population monte carlo energy redistribution”, Proc. ofEurographics Symposium on Rendering 2007, pp.287-296 (2007).
82) W. Jakob, “Light transport on path-space manifolds”, PhD Thesis,
Cornell University (2013).
83) A.

S.

Kaplanyan,

J.

Hanika,andC.

Dachsbacher,“The

Natural-Constraint Representation of the Path Space for Efficient Light
Transport Simulation”, ACM Transactions on Graphics (Proc. of
SIGGRAPH 2014), Vol.33, No.4 (2014).
84) T. Hachisuka, A. S. Kaplanyan, and C. Dachsbacher, “Multiplexed
Metropolis Light Transport”, ACM Transactions on Graphics (Proc. of
SIGGRAPH 2014) to appear (2014).
85) S. Kitaoka, Y. Kitamura, and F. Kishino, “Replice exchange lighttransport”, Computer Graphics Forum, Vol.28, No.8, pp.2330-2342(2009).

Tomoyuki NISHITA(Fellow)
He is a Professor in Department of Economic
Information of Hiroshima Shudo University and a
Director of UEI Research since 2013. He was aProfessor in the Department of Complexity
Science and Engineering (also in the Department
of Information Science) at University of Tokyo,
Japan since 1998, and now is a Professor Emeritusof The University of Tokyo. He taught at
Fukuyama University from 1979 to 1998. He wasan associate researcher in the Engineering
Computer Graphics Laboratory at Brigham Yong
University from 1988 to 1989. His researchinterests center in computer graphics including
lighting models, hidden-surface removal, andantialiasing. He received his B.E., M.E. and Ph.D.
in Engineering in 1971, 1973, and 1985,
respectively, from Hiroshima University. He
received Research Award on Computer Graphics
from Information Processing Society of Japan in
1987, and also received Steven A. Coons Award
from SIGGRAPH in 2005. He was a president of
The Institute of Image Electronics Engineers ofJapan in 2009.

Hisanari OTSU
He received the B.E. degree in Science from
Department of Information Science, Faculty of
Science, The University Tokyo in 2013. He iscurrently working toward the M.S. degree at The
University of Tokyo. His research interests include
physically based rendering and light transport
simulation.