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

Volume 35 (2016), Number 7

Multiple Scattering Approximation in Heterogeneous Mediaby Narrow Beam Distributions
Mikio Shinya† , Yoshinori Dobashi‡ , Michio Shiraishi§ , Motonobu Kawashima¶ and Tomoyuki Nishita∥

Abstract
Fast realistic rendering of objects in scattering media is still a challenging topic in computer graphics. In presence of participating media, a light beam is repeatedly scattered by media particles, changing direction and getting spread out. Explicitly evaluating this beam distribution would enable efficient simulation of multiple scattering events without involving costly
stochastic methods. Narrow beam theory provides explicit equations that approximate light propagation in a narrow incident
beam. Based on this theory, we propose a closed-form distribution function for scattered beams. We successfully apply it to theimage synthesis of scenes in which scattering occurs, and show that our proposed estimation method is more accurate thanthose based on the Wentzel-Kramers-Brillouin (WKB) theory.
Categories and Subject Descriptors (according to ACM CCS): I.3.3 [Computer Graphics]: Picture/Image Generation—

1. Introduction
Realistic simulation of light scattering can produce visually appealing images. Although rendering objects in scattering media
can be performed by computationally expensive Monte-Carlo stylemethods, such as path tracing [KVH84, YIC∗ 10] and photon mapping [JC98], fast realistic rendering of such objects is still a challenging topic in computer graphics.
Explicitly evaluating scattered beam distributions would efficiently simulate multiple scattering events, while avoiding recourse to stochastic procedures. This approach has been very
successful in subsurface scattering simulation. Methods calculate BSSRDFs/BSSTDFs (Bi-Directional Scattering Surface Reflectance/Transmittance Distribution Functions) and apply linear
filtering using distribution functions as filter kernels for the efficient simulation of sub-surface scattering. Dipole [JMLH01], multipole [DJ05], and improved multipole [dI11, FHK14] methods
have proven excellent in rendering layered scattering material suchas human skin. These methods are basically low-pass filters thatblur object textures. However, they experience difficulties when
dealing with heterogeneous media.
Scattering simulations in participating media are more challenging because 3D distributions are required and heterogeneous† Toho University/UEI Research‡ Hokkaido University/UEI Research§ Toho University¶ MomomoWorks Inc.
∥ UEI Research/Hiroshima Shudo Universityc 2016 The Author(s)⃝c 2016 The Eurographics Association and John
Computer Graphics Forum ⃝Wiley & Sons Ltd. Published by John Wiley & Sons Ltd.

structures are more visually important. It is known that WentzelKramers-Brillouin (WKB) theory [LRT82] can estimate beam distributions. Premoze et al. [PAS03] applied the WKB theory to
global illumination in participating media. Their proposed pathintegral method integrates scattered illumination using beam distributions as a filter kernel. Elek et al. [ERS13] applied the WKB
theory to simulate object blur due to scattering, and presented anefficient filtering scheme for fast rendering.
Given a deviated light path in media, which follows a curve in
3D space, the WKB theory estimates its probability. To obtain filter
kernels, the probabilities should be integrated over all possible lightpath curves starting at, and reaching specified points. This involves
very complex computations that require the introduction of major
simplifications on path curve shapes for practical computation. Thiscan cause large errors that limit their applications.
In this paper, we present explicit formulae for beam distributions
based on narrow beam theory [Ish78]. A naive application of thetheory to the filtering process involves 9D integrals for each pixel
calculation. We reduce the computation by evaluating moments andfitting a heuristic function to the moments. From analysis of a single scattering event, we found a distribution function based on theerror function (an integral of the Gaussian function) to be suitable
for this application. Reasonable agreement was experimentally observed between the results from stochastic photon counting and our
proposed estimation. We applied the distribution function to object
blur due to scattering, as in Elek et al.’s approach [ERS13], andconfirmed significant improvements.
The contributions of this paper are summarized as:• we introduce narrow beam theory, a powerful tool for estimating
scattering events, to the field of computer graphics;

M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions

The light transport equation for intensity I(r̂, ŝ) can be expressed
by

‫ܫ‬௢

0
scattering medium(∇ · ŝ)I(r̂, ŝ) = −(ρn σt )I(r̂, ŝ) +(ρn σs )
∫

ࢠො
‫ݖ‬

࢘

‫ݎ‬Ƹ

∫

p(ŝ, ŝ′ )I(r̂, ŝ′ )d ŝ′ ,

(1)

4π

p(ŝ, ŝ′ )d ŝ′ = 1.

(2)

4π

ߠ

࢙

‫ݏ‬Ƹ

Figure 1: Coordinate system and experimental arrangement.

where ŝ and r̂ represent the direction and position in 3D, ρn is thedensity of the media, σt and σs are the extinction and scattering
coefficients, and p is the phase function normalized to 1. The integral domain 4π indicates the unit sphere, and the domain is infinite
unless explicitly specified.
We set the intensity I to be the sum of the direct component Iriand indirect component Id
I = Iri + Id .

• we propose a heuristic distribution function based on narrow
beam theory;• we improve a screen-space scattering simulation by using ourdistribution function.
The proposed distribution function is also suitable for other image synthesis applications, such as shafts of light, Premoze et al.’s
global illumination [PAS03], and BSSTDF calculations in heterogeneous media.
2. Narrow Beam Theory
This section gives a brief introduction to narrow beam theory [Ish78]. The major symbols used in the equations are listed inTable 1, and the associated coordinate system is shown in Figure 1.

Let us assume that• the intensity I is concentrated near the ẑ direction (i.e., ŝ · ẑ =cos θ ≃ 1);• the variation in the phase function is smooth, allowing Id to beapproximated by a second-order Taylor expansion in direction ŝ
in the spherical integral in Equation (1).
We also assume the following situations, to obtain explicit solutions:• there is only one type of scattering particles in the media and σt ,σs , and p(s, s′ ) are spacially invariant;• the scattering media is structured into layers, and the particle
density may change along the z axis, described as ρ(z).
We setAs described in Section 1 of the supplemental document [PAS03,
APRN04, PARN04], the WKB theory estimates the probability fora given light path, but the path-space formulation is not a simple
and convenient tool to obtain distribution functions in 3D space.
The narrow beam theory, on the other hand, analytically approximates solutions in real 3D space and leads to distribution functions
in more explicit forms.
Table 1: Symbols used in narrow beam theory.
r̂, r, z
ŝ, s
σt
σs
σa
ρn
ζ
τ
κ
q
p(s)
P(κ)
I0 (r, s)
F0 (κ, q)position (r̂ = r + zẑ)
direction (ŝ = s + ẑ)
extinction coefficientscattering coefficientabsorption coefficient (σa = σt − σs )particle density
integral of ρn along path
optical depth (τ = σt ζ)
Fourier variable corresponding to s
Fourier variable corresponding to rphase function
Fourier transform of p
incident light distribution
Fourier transform of I0ŝ ≃ s + ẑ,r̂ = r + zẑ,

where ẑ is the unit vector in the z direction. Using this approximation, the spherical integral with ŝ in Equation (1) can be approximated by a planar integral with a 2D vector s. We also assume anaxially symmetric phase function around the forward direction:p(s, s′ ) = p(|s − s′ |).

(3)

Using these simplifications, the integro-differential equation in
Equation (1) can be replaced by a second-order partial differential
equation. Using Fourier transform, we can obtain the solution as
Iri (z, r, s) = 1/(2π)4

∫ ∞

−∞

∫ ∞

dκ

−∞

exp(−iκ · r) · exp(−is · q) ·F0 (κ, q + κz) · exp(−σt ζ(z))dq,
4

Id (z, r, s) = 1/(2π) exp(−ρn σa z)exp(−is · q)dq
Fd (z, κ, q) = σt F0 (κ, q)

∫ z
0
′

∫

∫

Fd (z, κ, q + κz) ·exp(−iκ · r)dκ,ζ(z) =0′′ρn (z )dz ,(5)ρn (z′ )P(q − z′ κ)exp(−σs ζ(z )) exp(−q2 Â(z, z′ ))dz′ ,

∫ z

(4)

(6)
(7)

c 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions

τ(z) = σt ζ(z),Â(z, z′ ) =∫ zz′A(z′′ )dz′′ ,(8)3.1. Incident Ray and Gaussian Phase Function(9)

Let us specify the properties of participating media and the lighting conditions, and proceed with the calculations. Since we aim toevaluate the distribution around a viewing ray, it is natural to assume that the incident light be a ray, which we represent by Diracdelta functions

A(z) = ρn (z)(σs /4) < θ2 >,

(10)

2

(11)

<θ > =∫

24π

θ p(s)ds/∫

p(s)ds,4π

where < . > represents an average, i indicates the imaginary unit,
F0 is the Fourier transform of the incident light distribution I0 (r, s)∫

F0 (κ, q) =exp(iq · s)ds∫

I0 (r, s) exp(iκ · r)dr,and P is the Fourier transform of the phase function p. τ(z) indicates optical depth at z, and κ and q are the Fourier variables corresponding to r and s, respectively. Integrals without a specified
domain are planar integrals over the whole plane. For more details,
readers are invited to refer to the book of Ishimaru [Ish78].
Note that Equation (5) involves a 5D integral. If we want to apply
filtering to the image of an object R(r, s) with the distribution Id
given by∫ ∫Id (r − r ′ , s′ )R(r ′ , s′ )dr′ ds′ ,′I0 (r, s) = δ(r)δ(s).
The phase function is assumed to be Gaussian† and is set top(s) = 1/(4πB) exp(−|s|2 /4B),∫

p(s)ds = 1.
Then, both Fourier transforms have analytical forms, such that
F0 (κ, q) = 1,

∫ z(π/Σ2G ) exp(−ρn σs z′ )dz′ (π/C2 (z′ )){}exp −(r − zs + Bz′ s/Σ2G )2 /(4C2 (z′ )){}= KG (z) exp −|s|2 /(4Σ2G ) · π2

′

0

0

The theory assumes that Id is distributed close to the ẑ direction.
If the variation in the reflection from the object R(r ′ , s′ ) is small,
we can approximate the filtering integral by′′dz′ exp(−ρn σs z′ )/(Σ2G (z′ ) ·C2 (z′ )){}exp −(r + (−z + Bz′ /Σ2G (z′ )s)2 /(4C2 (z′ )) ,
(16)

where
Â(z, z′ ) =′Id (r − r , s )R(r , s )dr ds
)
′ ′
′
≃
Id (r − r , s )ds R(r ′ , 0)dr′ .

Σ2G (z′ )2 ′Id (r − r ′ , s′ )ds′ ,z′A(z′′ )dz′′

KG (z) =(12)

we can reduce the filtering computation.

We will now detail these three reductions in the following sections.
c 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝′2

Bz [1 − B/Σ2G ],(1/2π)4 σs exp(−ρn σa z).

(17)
(18)
(19)
(20)

The direct component Iri can be calculated fromIri (z, r, s) = δ(r − zs)δ(s) exp(−τ(z)).

For further reduction, we took the following approaches:
1. Analytic Fourier transform: We selected the incident light I0 atz = 0 and the phase function p so that we could analyticallycompute their Fourier transforms F0 and P.
2. Moment calculation: We calculated the zeroth- and secondorder moments of Id , which, fortunately, can be calculated very
efficiently.
3. Fitting with a heuristic function: We carefully selected a heuristic function that gives a good approximation to the distribution.
We set the parameters of the function so as to have the same
calculated moment values.

′= B + Â(z, z ),C (z ) =In this case, if we can efficiently calculate the filter kernelΨ(r) =∫ z= (σt /4) < θ2 > (ζ(z) − ζ(z′ )),

∫ (∫

∫

(15)

Id (z, r, s) = KG (z) exp(−|s|2 )/(4Σ2G (z))

As seen in the previous section, drastic reductions of the computation complexity is necessary to apply narrow beam theory to image
synthesis.

′

2

Putting these expressions into Equation (5) yields a line integral

∫ z

′

′

P(q − z κ) = exp(−B(q − z κ) ).

3. Distribution Function

′

(14)

′

the 5D integral should be repeated for all r and s , requiring a
9D integral just to evaluate a single pixel. Consequently, we needfurther simplifications.

∫ ∫

(13)

(21)

Note that τ = (ρn σt )z when ρn is constant. More details can befound in Section 2 of the associated supplemental document.

3.2. Moments
To calculate the distribution, the line integral of Equation (16) has
to be evaluated for all r and s. In order to estimate it more quickly,† A sum of Gaussians would also be possible for a more accurate representation.

M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions

let us calculate the zeroth- to second-order moments of Id

0.5

∫ ∫

Id (z, r, s)dsdr,

M1 (z) = 0,

∫ ∫M2 (z) =

(22)

0.4

ρ2 Id (z, r, s)dsdr.

(24)

Elek

PI

(23)
RMS

M0 (z) =0.3

direct0.2

Here, the first-order moment (center of mass) is null because of its
symmetry.

0.1

When the density ρn is constant, we have closed-form formulae0singleProposed
0.5

M0 (z) = exp(−ρn σa z)(1 − exp(−ρn σs z)),(25)

M2 (z) = 4 exp(−ρn σa z)}[{exp(−ρn σs z) Az2 /a − 2B/a2 +{}]2B/a2 − 2Bz/a + (−A/a + B)z2 + Az3 , (26)a = ρn σs .

(27)

When density is not uniform but is dependent on z, we have aline integral
M0 (z) = exp(−σa ζ(z))(1 − exp(−σs ζ(z))),
∫ z

′

(28)

′

M2 (z) = 4 (ρn (z )σs ) exp(−σs ζ(z ))[ 0]C(z′ )2 + (z − Bz′ /Σ2G (z′ ))2 Σ2G (z′ )) dz′ .

0.6

0.8

0.9

1

albedoFigure 2: RMS errors and albedo (σs /σt ).

We regard this as a fundamental feature of the distribution function, and adopt a heuristic distribution function described by
Ψ(z, r) = α(z)(1/|r|)G(|r|/β(z)).

(34)

Parameters α and β can be determined from the moments Mi ,calculated using Equations (26) or (29)(29)

As seen in Equations (26) and (29), the moments can be efficiently
estimated in constant time in homogeneous media, and in linear
time in heterogeneous media. More details can be found in Section 3 of the supplemental document.

0.7α(z) = 2M0 /β2 (z),

(35)

2

(36)

β (z) = 3(M2 (z)/M0 (z)).

Using these formulae, we can determine the filtering kernel inconstant time for homogeneous density distributions, and in linear
time for heterogeneous distributions. The details are provided in
Sections 4 and 5 of the supplemental document.

3.3. Heuristic Distribution Function
The distribution function Ψ(|r|) is set according to the moments
M0 and M2 . This allows for efficient estimation of the distribution without integrating Equation (16) for each value of r and s.
Because the distribution functions include single-scattering components, the shape of the single-scattering distributions can help us
to design our heuristic function.
Let us calculate the single scattering component I1 of Id fromI1 (z, r, s) = (ρn σt )
∫
4π

∫ z
0

exp(−σt ρn y)dy

I0 (z − y, r + ys, s′ )p(s, s′ )ds′

= (ρn σt ) exp(−σt ρn z)
The integral over s
Ψ1 (z, r) =∫ z0

δ(r + ys)p(s)dy. (30)∫I1 (z, r, s)ds,(31)

can be calculated from
Ψ1 (z, r) = (ρn σs ) exp(−σt ρn z) · (1/B) ·√√(1/|r|)(1/ π)G(|r|/(z 2B)),(32)

using the finite integral of a Gaussian function (error function) G√ ∫ ∞G(x) = 1/ 2π
exp(u2 /2)du.
(33)
x

4. Experiments
4.1. Point Spread Functions
With the configuration shown in Figure 1, we randomly launchedmany photons at the origin along the z-axis, simulated their scattering, recorded their traces, and counted the number of photons passing through the cells (we refer to these simulations as photon counting). We compared the results calculated from our heuristic function and the moment equations (Equations (34), (35), and (36)),
with the results measured by photon counting. We fixed σt = 1.0
and used a Gaussian phase function with various spread anglesdθ2 = 4B. Figure 2 shows the measured RMS errors with respect
to albedo (σs /σt ) values. The phase function was fixed to dθ = 30
degrees. We also calculated the direct light component using Equation (21), the single-scattering component from Equation (32), thedistribution based on the path-integral (PI) method, and Elek et al.’s
version [ERS13]. The distribution functions used in the PI method
and Elek’s method are provided in Section 1 of the supplement
document (Equations (14) to (16)). As seen in Figure 2, for lower
albedo values, the errors caused by all methods are relatively small
because single-scattering components are dominant. Although errors grew up with higher albedo values in general, the proposed
method keeps a reasonable accuracy.
Figure 3 shows as color-coded the distributions measured byphoton counting and those calculated by our proposed method,c 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions݀ߠ ൌ ͵Ͳι݀ߠ ൌ ͸Ͳι

y1.0emitterz

camera

x
10-3

ȁ࢘ȁ
2.5

Figure 5: Simple scene used for the comparison.

ࢠ

density ߩ௡

5.0

1

Figure 3: Distribution functions measured by photon counting(path tracing), calculated by our (proposed) method, and by
the methods of (PI) Premoze et al. [PAS03] and (Elek) Elek et
al. [ERS13], all displayed with the same pseudo-colors. The spreadangle of the phase function, dθ, is 30 degrees (left) and 60 degrees(right). Note that the axes are permuted, compared to Figure 1.

5

10

0.8

Single

RMS

PI

0.5

Elek

Elek

Proposed

Path tracing

Figure 6: Resulting images of the square light emitter. Appliedcomputational methods are (from left to right) single scattering, thepath integral, Elek et al.’s, the proposed method, and path tracing.
The density is set to (from top to bottom) 1.0, 5.0, and 10.0.

0.7

0.6

Path integral

0.4

direct0.3

multiple-scattering terms. The figure indicates significant improvements by the proposed distribution functions.

single0.2

0.1

4.2. Rectangular Emitter Images

proposed
0
0

10

20

30

40

50

60

70

݀ߠ(degree)

To evaluate accuracy, we rendered images of a very basic scene
with path tracing and with our proposed method, and compared

Figure 4: RMS errors and spread angle of the phase function dθ.
density=10

Figure 4 shows root-mean-square (RMS) errors. The errors fromour proposed method are mostly around 10%. Observe how theaccuracy has largely been improved by taking into account thec 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝

0.8

1
0.8

Elek
PI

0.6

RMS

RMS

with two different phase functions. The circular integral 2π · r · Ψ
is shown instead of Ψ, because Ψ itself is singular at the origin (r = 0). The distributions become broader with respect to thedepth z and the variance of the phase function dθ. Observe how the
distributions calculated by our proposed method agree much betterwith the simulated reference solutions than the other two methods.

density=5

single1

0.4
proposed

0.4

0.2

0.2

0

0
15

25

35

45

single

Elek

0.6
PI

proposed15
݀ߠ (deg.)

25

35

45

Figure 7: RMS errors for particle density values (ρn ) of 10.0 and
5.0.

M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributionsȟ߫௜ ȟ߫௜ାଵ…..
pixel
viewing ray

‫ݖ‬௠

‫ݖ‬଴

෍ ȟ߫௜

screen‫ݖ‬଴objectobjectmediumscreen‫ݖ‬஻ ൌ ‫ݖ‬ሺ‫)’ݔ‬ߩ௡

Figure 8: Calculation of the distribution functions in (left) homogeneous media and in (right) heterogeneous media.

Ȳ஻ ሺ‫ݖ‬஻ ǡ ‫ݎ‬ሻ‫ݎ‬Ƹ஻ ሺ‫ݖ‬ሻ

‫’ݔ‬
‫ݔ‬

‫ݎ‬Ƹ஺ ሺ‫ݖ‬ሻQB Ȳ஺ ሺ‫ݖ‬஺ ǡ ‫ݎ‬ሻ㼐㻌Ȳ஺ ሺ‫ݖ‬஻ ǡ ‫ݎ‬ሻmediumscreenFigure 9: Kernel selection.
the synthesized images. Figure 5 illustrates the configuration of thescene, where a red square light emitter is located in a transparent
box filled with homogeneous scattering media. For path tracing, we
traced about 10K paths per pixel, which took a couple of hours per
image. For reference purposes, we also rendered images with single
scattering, path-integration filter kernels, and Elek et al.’s version(Equation 12 in the supplemental document).
Figure 6 shows the corresponding results. When the density islow, the direct ray components are dominant, and differences between the generated images are rather subtle. However, with increases in density, differences become more significant: the images
produced by single scattering and path integration are too dark,
while our proposed method captures the blurred light caused bymultiple scattering. Figure 7 shows RMS errors measured in two
conditions. The filtering based on path integration does a better job
than the single-scattering approximation, but still causes large errors at higher densities. Elek et al.’s distribution function has larger
errors than those of the path integral in general, possibly becauseof its additional simplification. With our proposed method, on theother hand, errors lie around 10 to 20%, and we observe reasonableagreements with the path-tracing results.

5.1. Method
5.1.1. Distribution
The distribution Ψ can be calculated from the moments M0 and
M2 . In the case of homogeneous media, the moments are calculated with Equation (24), using only depth values of media andobject surfaces. To do so, we save depth values in a preprocessing
stage and refer to these values when the moments are evaluated, as
illustrated in Figure 8 (left).
In heterogeneous media, the moments are calculated with Equation (29), where the integral is calculated by referring to the accumulated density ζ. Therefore, we perform the integral by setting
sample planes parallel to the screen, and simultaneously sum ζ andM2 in an alpha-blending fashion, as illustrated in Figure 8 (right).

‫ݖ‬஺ ൌ ‫ݖ‬ሺ‫)ݔ‬

‫ݔ‬

‫ݔ‬

(a) gathering filter(b) distributing filter

Figure 10: Filtering implementation with the gathering and distributing schemes.

5.1.2. Filtering
Elek et al. [ERS13] propose an efficient multi-scale filtering
scheme for real-time execution. In this paper, we employ a simpler
and more generic filtering scheme to avoid potential artifacts.
Consider the configuration shown in Figure 9, in which pixel
value I f ilter at pixel x is evaluated fromI f ilter (x) = ∑ w(x, x′ )R(x′ ),5. Application to Screen-space Scattering Simulations
To demonstrate advantages of the proposed distribution function, it
is applied to screen-space scattering simulations, as introduced by
Elek et al. [ERS13].

QA

x′

′

where x and x denote pixel positions, zA and zB are the object depth
values, and w and R represent the filter kernel and the object image,respectively.
Along the central ray r̂A (z) passing through pixel x, the distribution ΨA (z, |r|) varies with respect to depth z. Without loss of generality, let us assume zB < zA . In this case, some of the scatteredviewing rays reach object point B, which can be evaluated from the
distribution ΨA (zB , r). Therefore, we can set weight w asw(x, x′ ) = ΨA (zB , dr ),dr = |r̂A (zB ) − r̂B |.

(37)
(38)

Furthermore, if we assume local similarity for nearby pixels
ΨA (z, |r|) ≃ ΨB (z, |r|),
we can set w asw(x, x′ ) ≃ ΨB (zB , d)(39)

c 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions

pre-process

object image synthesis(SPP method)

calculation of
moments, ߫ǡ Ƚǡ Ⱦgathering filter

distributing filter
calculation of scattered external light(optional)
blending

Figure 11: Pipeline of our rendering system.

for zB < zA . This switching is necessary to avoid “illuminance leaking” as discussed by Elek et al. [ERS13].
In general, filtering can be implemented in two ways: gathering
contributions from other pixels, or distributing their contribution
to other pixels, as shown in Figure 10. We implemented filteringin both ways depending on the z-value comparison. For those pixels with larger depth, we apply the gathering filter; for those withsmaller depth, we use the filtering result from the distributing filter.
5.1.3. Implementation
We implemented a rendering system, which is outlined in Figure 11.
Preprocess: In a preprocessing stage, the shaded image col[] anddepth map zbu f [] of the scene’s objects are calculated using conventional methods. The shaded image will be filtered by the distribution function estimated in the following main process. Although illumination from an external source may be evaluated bysimple attenuation, we adopted the simplified plane-parallel (SPP)
model [SSD∗ 10] to capture multiple-scattering features of the light
distributed in the scene.
Calculation of moments: Next, the moments and parameters of
distribution Ψ, α and β, are calculated at each pixel according toEquations (25), (29), (35), and (36). We implemented this in threepasses (Algorithm 1). First, we prepare a set of sample planes perpendicular to the viewing direction. Each sample plane is drawn
and the density ρn is accumulated via alpha-blending up to the object depth, resulting in ζ0 . Second, we accumulate the second-ordermoments M2 , using ζ0 . Finally we obtain the parameters of distributions, α and β.
Gathering filter: Using the calculated filter kernels, both gathering filter and distributing filter are applied to every pixel. For thosepixels with a larger depth value than that of the pixel in process,we apply the gathering filter (Algorithm 2). We draw points corresponding to each pixel through the vertex shader, and the fragc 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝Algorithm 1 Calculation of moments
Require: zbu f [], ρn [] /* [] indicates an array */Ensure: α[], β[], ζ0 []initialize: M0 [], M2 []ζ0 [], ζ ⇐ 0/* Draw sample planes through vertex shader */for all sample planes z p lane do/* Process in fragment shader */for all pixel p do
if zbu f [p] > z plane then
ζ0 + = ρn ∗ dz
end ifend forend forfor all sample planes z plane dofor all pixel p do
if zbu f [p] > z plane thenz = zbu f [p], z′ = z plane
ζ+ = ρn ∗ σt ∗ dz
Â = σs B(ζ0 [p] − ζ)
Σ2G (z′ ) = B + Â /* Eq. (18) */C2 (z′ ) = Bz′2 [1 − B/Σ2G ] /* Eq. (19) */e0 = (ρn (z′ )σs ) exp(−σs ζ)
e1 = [C(z′ )2 + (z − Bz′ /Σ2G )2 Σ2G )]M2 [p]+ = e0 ∗ e1 /* Eq. (29) */end ifend forend forfor all pixel p do
M0 [p] = σt exp(−σa ζ(z))(1.0 − exp(−σs ζ0 ))β[p] = 3 ∗ (M2 [p]/M0 ) /* Eq. (36) */α[p] = 2 ∗ M0 /β[p] /* Eq. (35) */end forment shader processes the filtering. We limit the size of the filter at pixel p by simply multiplying β[p] by user-specified factor
cuto f f , which was set to 1.0 during our experiments. For all pixels
p′ within the filter, we calculate the filter weight, weight, and sum
up the weighted colors if the depth zbu f [p′ ] is larger than zbu f [p].
The error function G(x) can be evaluated by Hastings’ approximation [Has55], for example.
Distributing filter: For those pixels with a larger depth value thanthat of the pixel in process (p), we apply the distributing filter (Algorithm 3). In the vertex shader, we calculate the filter size for eachpixel and draw a rectangle centered at the location of the pixel inprocess covering the filter. In the fragment shader, we calculate
the filter weight if the depth zbu f [p′ ] is larger than zbu f [p]. The
weighted colors are summed up by alpha-blending.
Calculation of scattered external light: Since this screen-space
filtering cannot capture effects from light sources out of the image, we optionally calculate their contributions using conventional
methods. This may be conducted by a standard ray-marching
method capturing single-scattering components. In our implementation, we adopted the SPP method as in the pre-processing phase.

M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions

Algorithm 2 Gathering filter
Require: zbu f [], col[], α[], β[], cuto f f
Ensure: Igather []initialize: Igather ⇐ 0/* draw pont and process in fragment shader */for all pixel p do
z0 = zbu f [p]r = WorldCoord(p) /* transform to world coordinate */FilterSize = cuto f f ∗ β[p] /* filter size */for all pixel p′ within the filter centered at p do
if zbu f [p′ ] > z0 thenr ′ = WorldCoord(p), dr = |r ′ − r|weight = (α[p]/dr )G(dr /β[p]) /* Eq. (34) */Igather [p]+ = weight ∗ col[p′ ]end ifend forend for(a) path-tracing (2K rays/pixel)

Algorithm 3 Distributing filter
Require: zbu f [], col[], α[], β[], cuto f f
Ensure: Idistr []initialize: Idistr [] ⇐ 0/* draw rectangles centered at each pixels via vertex shader */for all pixel p do
z0 = zbu f [p], I0 = col[p]r = WorldCoord(p) /* transform to world coordinate*/
FilterSize = cuto f f ∗ β[p] /* filter size *//* draw a rectangle and process filtering in fragment shader */for all pixel p′ within the filter rectangle doif zbu f [p′ ] > z0 thenr ′ = WorldCoord(p), dr = |r ′ − r|weight = (α[p]/dr )G(dr /β[p]) /* Eq. (34) */Idistr [p′ ]+ = weight ∗ I0
end ifend forend for(b) proposed

Blending: Finally, we bend the calculated components into the final image (Algorithm 4).
Algorithm 4 Blending
Require: ζ0 [], col[], Idistr [], Igather [], Iext []Ensure: I f inal []/* draw points corresponding to pixels*/
for all pixel p do
T = exp −ζ0 [p]I f inal [p] = T ∗ col[p] + Idistr [p] + Igather [p]Iext [p]end for(c) Elek’s
Figure 12: Comparison for the CityLight images.

5.2. Synthesized Images
We synthesized images of scenes with heterogeneous scattering
media.
c 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions(a) CityLight(a) path-tracing (4K rays/pixel)(b) Westminster(b) proposed

Figure 14: Example from the video sequence.

Figure 12 shows a night scene, rendered by (a) path tracing,(b) the proposed distribution function, and (c) Elek et al.’s [ERS13]
approximation. In the scene, buildings with colored lit windows are
partially immersed in thick heterogeneous fog. Light in the lower
parts of buildings are more attenuated and blurred by the fog, which
was successfully simulated by our proposed method, which reasonably agrees with the result from path tracing. The RMS difference
from the path-traced image is 28.6%. Elek et al.’s kernel is moreconcentrated in its center, and over-estimates the blurring near lightsources; its RMS difference is 52.7%.

(c) Elek’s
Figure 13: Comparison for the Westminster images.

Figure 13 shows another familiar example, where a famous clock
tower stands in heavy heterogeneous fog. The image rendered byour proposed function reasonably agrees with the one by path tracing, demonstrating significant improvements.
We also created fly-through sequences using these scenes. Representative frames from the sequences are shown in Figure 14.
Since this screen-space method cannot capture effects from outof-screen light sources, we set a weak parallel light, calculated its
scattered intensity by the SPP method [SSD∗ 10], and blended it toc 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝M. Shinya, et. al. / Multiple Scattering Approximation in Heterogeneous Media by Narrow Beam Distributions

the filtered images. Varying fog density was created from a fluid
simulation, and stored as voxel data.
As an indication of typical rendering times, complete processing times were 110 msec and 100 msec, for Figure 14-(a) and (b),respectively, on a Windows PC (Intel Core i7-3770K @3.50GHz,
GeForce GTX 690), while path tracing took several days. The image size is 640 × 480. The most time-consuming part of the algorithm is the linear filtering process, especially for large kernel
sizes. To accelerate this process, we utilized multi-resolution images, similar to a MIP map, and applied a filter to a coarser image
when the kernel size is large (larger than 9 pixels in our implementation).
5.3. Limitations
As mentioned in Section 2, the theory assumes narrow and smooth
distribution of scattering light as well as layered, single-type scattering media. Screen-space algorithms, in general, assume that everything necessary is visible on screen and do not count contributions from invisible/out-of-screen light sources. These assumptions
limit applicable scenes and the accuracy of the method, for example, in the following situations:• Mixed material: Scattering media composed with more than two
different materials, such as milk poured into coffee, may havespatially varying mixture ratios, which makes the phase function
spatially variant and hard to be correctly handled by the theory.
• Occlusion/shadowing: Occlusion causes sharp changes in light
distribution, which the theory cannot correctly handle. The proposed method tries to reduce artifacts by switching gathering anddistributing filters based on depth comparisons, but errors still
remain.
• Large variation of density: When the density function is not layered and variation perpendicular to viewing direction is large, theestimated distribution may involve large errors.
• Very small particles: Particles much smaller than the wavelength
tend to have strong backscattering properties as in Rayleigh scattering. The theory ignores backward-scattering components andthe method may produce large errors.
• External light sources: Lighting effects due to external lightsources out of screen are not taken into account by the filtering process. The proposed method separately estimates the external light components and simple blending is applied, which
may cause errors.
Note that these general limitations are related to the filtering approach, which affect also Elek’s method [ERS13].
6. Conclusion
In this paper, we proposed useful beam distribution functions inheterogeneous scattering media. Based on the narrow beam theory, we formulated the moments of the distribution and proposed aheuristic distribution function to fit to the calculated moments. Experiments show that the proposed functions are more accurate than
WKB-based functions previously used. The proposed function was
successfully applied to screen-space simulations of object blurring
due to multiple scattering, and demonstrated significant image improvements over previous functions.

The proposed distribution functions can be also applied to multiple scattering of directional light sources, such as car headlights,
shafts of light through clouds, or city light beams in the night sky.
They are also suitable for other image synthesis applications, suchas Premoze et al.’s global illumination [PAS03], and BSSTDF calculations in heterogeneous media.
Acknowledgements
The authors would like to thank the reviewers for their careful reviews and helpful comments. They also thank Pierre Poulin for his
technical dicussions and linguistic advices, and Tatsuki Kawaguchi
for his scene models used in the experiments. This research waspartially supported by JSPS KAKENHI, Grant-in-Aid for Scientific Research on Innovative Areas, Grant Number JP15H05924.
References[APRN04] A SHIKHMIN M., P REMOZE S., R AMAMOORTHI R., NAYAR
S.: Blurring of light due to multiple scattering by the medium: a path Integral approach. Tech. Rep. CUCS-017-04, Columbia University, 2004.
2[dI11] D ’E ON E., I RVING G.: A quantized-diffusion model for renderingtranslucent materials. ACM Trans. on Graphics 30, 4 (2011), 56. 1[DJ05] D ONNER C., J ENSEN H. W.: Light diffusion in multi-layeredtranslucent materials. ACM Trans. on Graphics 24, 3 (2005), 1032–1039.
1[ERS13] E LEK O., R ITSCHEL T., S EIDEL H.: Real-time screen-spacescattering in homogeneous environments. IEEE Computer Graphics andApplications 33, 3 (2013), 53–65. 1, 4, 5, 6, 7, 9, 10[FHK14] F RISVAD J. R., H ACHISUKA T., K JELDSEN T. K.: Directional
dipole model for subsurface scattering. ACM Trans. on Graphics 34, 1(2014), 5:1–12. 1[Has55] H ASTINGS C.: Approximations for Digital Computers. Princeton University Press, 1955. 7[Ish78] I SHIMARU A.: Wave propergation and scattering in random media. Academic Press, 1978. 1, 2, 3
[JC98] J ENSEN H. W., C HRISTENSEN P.: Efficient simulation of lighttransport in scenes with participating media using photn mapping. In
Proc. SIGGRAPH (1998), ACM, pp. 311–320. 1[JMLH01] J ENSEN H. W., M ARSCHNER S. R., L EVOY M., H ANRA HAN P.: A practical model for subsurface light transport. In Proc. SIGGRAPH (2001), ACM, pp. 511–518. 1[KVH84] K AJIYA J., VON H ERZEN B. P.: Ray tracing volume densities. In Computer Graphics (Proc. SIGGRAPH) (1984), vol. 18, ACM,
pp. 165–174. 1[LRT82] L ANGOUCHE F., ROEKAERTS D., T IRAPEGUI E.: Functional
Integration and Semiclassical Expansions. D. Reidel Publishing Company, 1982. 1[PARN04] P REMOZE S., A SHIKHMIN M., R AMAMOORTHI R., NAYAR
S.: Practical rendering of multiple scattering effects in participating media. In Proc. Eurographics Symp. on Rendering (2004), Eurographics,
pp. 130–140. 2[PAS03] P REMOZE S., A SHIKHMIN M., S HIRLEY P.: Path integration
for light transport in volumetric materials. In Proc. of Eurographics
Symp. on Rendering (2003), Eurographics, pp. 130–140. 1, 2, 5, 10[SSD∗ 10] S HINYA M., S HIRAISHI M., D OBASHI Y., I WASAKI K.,
N ISHITA T.: A simplified plane-parallel scattering model and its application to hair rendering. In Proc. Pacific Graphics (2010), ACM, pp. 85–
92. 7, 9[YIC∗ 10] Y UE Y., I WASAKI K., C HEN B.-Y., D OBASHI Y., N ISHITA
T.: Unbiased, adaptive stochastic sampling for rendering inhomogeneous participating media. ACM Trans. on Graphics 29, 6 (2010), 177:1–
8. 1
c 2016 The Author(s)⃝c 2016 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum ⃝