Rendering discrete participating media using geometrical optics approximation

We consider the scattering of light in participating media composed of sparsely and randomly distributed discrete particles. The particle size is expected to range from the scale of the wavelength to several orders of magnitude greater, resulting in an appearance with distinct graininess as opposed to the smooth appearance of continuous media. One fundamental issue in the physically-based synthesis of such appearance is to determine the necessary optical properties in every local region. Since these properties vary spatially, we resort to geometrical optics approximation (GOA), a highly efficient alternative to rigorous Lorenz—Mie theory, to quantitatively represent the scattering of a single particle. This enables us to quickly compute bulk optical properties for any particle size distribution. We then use a practical Monte Carlo rendering solution to solve energy transfer in the discrete participating media. Our proposed framework is the first to simulate a wide range of discrete participating media with different levels of graininess, converging to the continuous media case as the particle concentration increases.


Introduction
Rendering participating media is a long-standing problem in computer graphics, with much effort devoted to solving it plausibly and efficiently [1,2]. From the physics point of view, radiative transfer is rather complicated and must be rigorously derived from Maxwell's electromagnetic theory. To make simulation tractable, various compromises have been made in physically-based rendering over the past decades. Two main assumptions are independent scattering [3] and local continuity (or statistical homogeneity [4]). The first means that the particles forming the medium are far apart from each other and are mutually unaffected. Recently, this assumption has been relaxed in computer graphics by incorporating spatial correlations between scatterers into radiative transfer frameworks [5][6][7][8][9], leading to non-exponential attenuation of light.
The second assumption of local continuity implies that the medium is homogeneous and compact in each differential volume even if macroscopic heterogeneity exists. Under this circumstance, the light is not sensitive to the discrete spatial distribution of the scatterers, but only to their local average properties. Consequently, light scattering phenomena take place at every point of the medium, resulting in locally smooth rendering. However, many participating media are composed of separate particles distributed randomly within a given volume, and scattering happens only at particle positions [4,10]. These facts indicate that the assumption of continuous media only holds when the particle size is much smaller than the resolution of the sensor (e.g., the human eye) and the quantity is sufficiently large [11]. Otherwise, individual grains can be observed when zooming in. To overcome this constraint, the graininess of the medium should be taken into consideration.
Several recent studies in computer graphics [12][13][14] have noticed the graininess in rendering granular materials. They typically rely on explicit geometries and precomputed transport functions to capture the appearance of discernible grains. As these approaches are designed for a certain number of very large particles using geometric optics, the scattering behaviors of particles forming the media are rather limited. For instance, diffraction, which dominates light scattering from small particles, is often ignored. In this paper, we attempt to put forward a more general approach to model and render discrete participating media with a large number of scatterers whose particle size distribution (PSD) ranges widely. These media are omnipresent in natural and artificial environments, such as airborne dust, blowing snow, powder suspensions, and air bubbles in liquid.
Unlike conventional continuous media, these discrete media have spatially-varying optical properties (e.g., the extinction coefficient) that cannot be determined in advance but must be evaluated on-the-fly. These properties are closely related to the scattering behavior of each individual particle and fluctuations in PSDs. To derive necessary optical properties for any local region of a given discrete participating medium in a physically-based manner, one can resort to Lorenz-Mie theory [15,16] which can provide substantial realism in rendering participating media [17]. However, numerical evaluation of the Lorenz-Mie coefficients is known to be difficult and time-consuming as the particle size increases [18]. To ameliorate this issue, we introduce geometrical optics approximation (GOA) [18,19] and use it to simplify the computation of light scattering when the particle is sufficiently large. We show how to perform a smooth transition between Lorenz-Mie theory and GOA in computing optical properties, enabling both high accuracy and low computational cost.
Due to the variations in local PSDs, optical properties exhibit multi-scale effects with respect to the scene configuration. We derive a novel multi-scale volumetric rendering equation (VRE) and propose a practical Monte Carlo rendering solution to solve it. Our solution only relies on the position and the radius of each particle distributed randomly according to given PSDs, avoiding cumbersome geometric modeling and lengthy precomputation. Experimental results verify that the proposed solution is able to capture the distinct grainy appearance of discrete participating media and can guarantee temporal coherence in animation. We also show that it converges to continuous media in the limit of particle concentration.
In summary, the main contributions of this paper are: • a general and physically-based framework for modeling and rendering discrete participating media, considering diffraction, polarization, a wide range of PSDs, etc., • the use of GOA for efficient and accurate evaluation of the multi-scale bulk optical properties in any local region of a medium, and • a new Monte Carlo rendering solution that captures both low-frequency haziness and highfrequency graininess in discrete participating media.

Participating media rendering
Rendering participating media is a challenging but important problem, which requires efficiently solving the VRE [1,3] by means of Monte Carlo path integration [2,[20][21][22][23], photon density estimation [24][25][26][27][28][29][30], or a combination of both [31]. In our current framework, we choose Monte Carlo path integration because of its elegant simplicity, generality, and accuracy. This technique operates by stochastically constructing a large number of light paths between receivers and emitters to simulate the light transport in the scene. To facilitate the query of particles along paths, we augment each ray with a cylinder, in a way similar to the photon beam approach [25,26,31]. Multiple importance sampling [32,33] algorithms are beneficial for reducing the large variance caused by Monte Carlo sampling.

Detailed volumetric modeling
Since the original VRE is only a rough approximation to the real radiative transport in participating media, the range of appearance that can be faithfully simulated is limited. The recent trend in computer graphics is to try to capture more details in the volume by relaxing the assumptions in the original VRE, enriching the range of achievable appearances. For example, to account for angular anisotropy, the VRE can be extended with local directional dependency based on a microflake model [34][35][36]. It is also possible to extend the VRE to simulate the effects of spatial correlations [5,6,9], yielding non-exponential attenuation of light. Notably, these methods still assume the media to be statistically homogeneous within each differential volume [4], ignoring any sub-pixel details.
To handle participating media with complex 3D structures, volumetric representations of explicit geometries have been widely used. By capturing the geometric and optical properties of a fabric down to the fiber level, micro-appearance models, described using high-resolution volumes, offer stateof-the-art rendering for fabrics and textiles [37][38][39][40][41]. Unfortunately, these methods are highly dataintensive and plagued by heavy computation. To improve the performance while maintaining good accuracy, some downsampling strategies [42,43] have been developed. Current rendering solutions for granular materials are also based on explicit geometries and pre-captured optical properties of each individual grain [12][13][14]. They generally employ shell tracing to make large jumps inside media. Even so, the computational cost is still high. The goal of this work is to develop a general framework for handling participating media with graininess which also allows rapid computation and convenient usage.

Glittery surface simulation
Our work is also closely related to the simulation of glints on surfaces. Yan et al. [44,45] suggested using explicit high-resolution normal maps to model sub-pixel surface details and successfully simulated spatially-varying glints with a patch-based normal distribution function. Subsequent work adopted a wave optics model to achieve more accurate results with noticeable color effects [46]. Other methods focus specifically on capturing spatially-varying highlights from scratched surfaces, under either geometric optics [47][48][49] or wave optics [50]. To ease the burden of computation and storage, Kuznetsov et al. [51] proposed to learn high-frequency angular patterns from existing examples, using a generative adversarial network (GAN). Jakob et al. [52] addressed the problem of glittery surface simulation using a purely procedural approach which requires far less storage and supports on-the-fly point queries. This approach has been extended to incorporate iridescence [53] and allow fast global illumination [54].

Lorenz-Mie theory and GOA
Lorenz-Mie theory [15,16] develops a rigorous solution to the problem of light scattering by spherical particles. It was introduced to the graphics community by Rushmeier [55] to accurately simulate the physics of light transport in participating media. Later, Callet [56] used this theory to model pigmented materials consisting of pigmented particles in a transparent solvent. Atmospheric phenomena, such as halos and rainbows, are especially well handled by this theory [57][58][59][60]. Frisvad et al. [17] generalized the original Lorenz-Mie theory and used it to compute the appearance of materials with different mixed particle concentrations. Though accurate, this theory is computationally expensive. We show that GOA [18,19,[61][62][63][64][65] is much more efficient than Lorenz-Mie theory in computing optical properties of individual particles in various media, especially when the particle size is large. Moreover, the computation of GOA can be made in nonideal situations including absorbing particles [65] and aspherical particles [61,[66][67][68]. Compared to Lorenz-Mie theory, GOA is less explored in computer graphics. We choose GOA in our framework, taking advantage of its high performance.

Preliminaries
We first consider light scattering by a single particle. Meanings of some abbreviations used in the paper are listed in Table 1. We suppose that the particle is approximately spherical and has a set of physical properties including its radius r and a refractive index η p . Currently, we assume that particles forming the medium have the same composition and only their sizes vary. In this case, the refractive index η p is fixed. Supposing that the host medium has the refractive index η m , we can define the relative refractive index of the particle as η = η p /η m . The size of a spherical particle may also be expressed in terms of the dimensionless size parameter α = kr = 2πη m r/λ, where k is the wave number defined by k = 2πη m /λ To describe scattering, we need two scattering amplitude functions: S 1 (θ, ϕ, r) and S 2 (θ, ϕ, r), where θ is the scattering angle and ϕ is the azimuth angle. The subscripts 1 and 2 denote perpendicular and parallel polarization, respectively. For spherical particles, S 1 and S 2 are invariant with respect to ϕ, but they change depending on the radius r. For unpolarized light, these two functions define the phase function of a single particle as [4]: which is properly normalized by the scattering cross section C s : (2) Another important property of the particle is the extinction cross section C t which is given by with S(0, r) = S 1 (0, r) = S 2 (0, r). The notation Re takes the real part of a complex number. For particles with absorption, the absorption cross section is given by C a (r) = C t (r) − C s (r). As seen, once the scattering amplitude functions S 1 (θ, r) and S 2 (θ, r) are available, we can easily find the scattering, extinction and absorption cross sections as well as the phase function of the particle. For scattering of light, an electromagnetic wave, from a homogeneous spherical particle, exact solutions of the two scattering amplitude functions are given by Lorenz-Mie theory [15,16]. Its accuracy has been validated against real measurements in various literature [69,70].
As a rigorous and general electromagnetic treatment of light scattering by spherical particles, Lorenz-Mie theory can precisely handle a wide range of particle sizes. However, as the particle size increases, numerical calculations of the Lorenz-Mie coefficients become very tedious and time-consuming, due to the fact that the number of terms to be computed in the series for S 1 (θ, r) and S 2 (θ, r) is proportional to the size parameter α [17,18]; it is suggested that an appropriate number of terms to sum is |α| + 4.3|α| 1/3 + 1 [71]. For this reason, simpler approximate expressions are useful to reduce the computational complexity. In the case that the particle size is large with respect to the wavelength of the illuminating light, geometrical optics approximation (GOA) [18,19,[61][62][63][64][65] provides a simplified but good solution.

Geometrical optics approximation
Within the framework of GOA, light scattering is calculated by a superposition of classical diffraction, geometrical reflection, and transmission. The diffraction is independent of the particle's composition (i.e., its refractive index). Its amplitude functions for the forward direction are readily described by Fraunhofer diffraction theory as [18]: where J 1 is the first-order Bessel function.
Leaving out diffraction, a light ray hitting a spherical particle at an incident angle θ i is partially reflected and partially refracted depending on the properties of the interface, as sketched in Fig. 1. The refracted ray may undergo a number of internal reflections before leaving the particle. For each emerging ray, we use an integer p to denote the number of chords it makes inside the spherical particle. Obviously, the externally reflected ray has p = 0 while the other rays are transmitted with p − 1 internal reflections. The angle of deflection θ p between the pth emerging ray and the direction of the incident ray is given by with sin θ i = η sin θ t according to Snell's law. The scattering angle θ is further determined by the deflection angle θ p as where q = 1 indicates that the incident ray hits the particle on the upper hemisphere and q = −1 on the lower, and l is an integer ensuring that the scattering angle θ is restricted to the range 0-π. Clearly, reflected and transmitted rays depend on the shape and composition of the particle. Their scattering amplitudes for each polarization can be derived to be [10,19]: Here, the fraction j (θ i ), which is due to the reflection and/or refraction for an emergent ray of order p, is defined as and φ f due to focal line: Putting together S D,j (θ, r) and S (p) j (θ, r), we get the total amplitude functions of GOA as for j = 1, 2. These expressions can be evaluated quite efficiently.
In GOA, an analytical expression for C t (r) can be derived: where P = {1, 3, 5, · · · }. This expression is fast to evaluate and captures ripple structures well [72]. For absorbing particles, we certainly have C a (r) = C t (r) − C s (r) > 0. Within GOA, the absorption cross section C a (r) is faithfully approximated by [4]: in which η r and η i are the real and imaginary parts of η, respectively. The scattering amplitude functions S 1 (θ, r) and S 2 (θ, r) are also slightly different. To our knowledge, this is the first time that GOA has been introduced into computer graphics. The benefit of GOA is its comparable accuracy to rigorous Lorenz-Mie theory and its high efficiency in describing and evaluating key optical properties (e.g., the extinction and absorption coefficients) that are necessary for rendering.

Number of terms
An infinite summation is needed in principle in GOA to compute S 1 and S 2 . However, unlike the one needed in Lorenz-Mie theory, the number of terms needed is independent of the particle size, and a small p suffices in most cases. As shown in Fig. 2, when calculating (|S 1 | + |S 2 |)/2 with p = 3 in GOA, we get an almost identical curve with that for p = 100, irrespective of the particle size. This is because higher-order reflections (p > 3) carry much less energy and have very little impact on the scattered light intensities. Thus, p may be safely set to 3 in what follows.
Moreover, we can further simplify the extinction cross section to by setting p = 1, since C t only relies on the value of S 1 (or S 2 ) evaluated at θ = 0, and light rays with p > 1 contribute little to forward scattering: see Fig. 3, where the C t curves for p = 1 (green) and p = 3 (red) are virtually indistinguishable for a very wide range of r.  Variation of the extinction cross section Ct as a function of r. We compare our GOA calculation (red: p = 3, green: p = 1) to that of Lorenz-Mie theory (blue) with η = 1.33 and λ = 0.6 μm. The ripple structures [72] are well captured by GOA as clearly shown in the left diagram.

GOA and Lorenz-Mie theory
To investigate the range of validity of GOA in simulating the scattering patterns of spherical particles, we compare its results with the rigorous Lorenz-Mie results for a wide range of particle radii in Fig. 4. The wavelength of the incident rays is set to 0.6 μm and the relative index of refraction is η = 1.33 in all calculations. Figure 4 reveals that the scattering amplitude distributions by GOA align well with those obtained with Lorenz-Mie theory for large particles with r > 1 μm. The agreement of these two methods is especially good in almost all directions when the radius is large enough (e.g., r = 100 μm). However, when r ≈ 1 μm, some discrepancies between the two methods appear. These discrepancies become large as the particle radius decreases further.
To show the influence of these discrepancies on the perception of translucent appearance, we rendered a smooth medium comprising monodisperse particles of radius r. We determined the phase function using either Lorenz-Mie theory or GOA according to r; these phase functions were precomputed and stored in tables. The extinction coefficient was set to a constant 0.6 for a fair comparison. The synthesized images are presented in Fig. 5 with r set to 1 and 2 μm. Clearly, the differences of phase functions between Lorenz-Mie and GOA in the case of r = 1 μm result in inconsistent appearance. However, this inconsistency almost disappears completely as r increases to 2 μm, although there are some mismatches on the backward peaks of S 1 (θ) (see Fig. 4(c); S 2 (θ) is similar). Figures 4 and 5 together validate choosing GOA to compute S 1 (θ) and S 2 (θ) when r 2 μm.
To further show the similarity between Lorenz-Mie theory and GOA in computing S 1 (θ) and S 2 (θ), we report their relative mean squared error (RelMSE) in Fig. 6(a). RelMSE is computed on (|S 1 (θ)| + |S 2 (θ)|)/2 for different relative refractive indexes. The results confirm that subtle errors exist when r is large: despite some fluctuations, these calculations of GOA exhibit errors of less than 0.01 as compared with exact Lorenz-Mie calculations when r 2 μm.
When calculating (|S 1 (θ)| + |S 2 (θ)|)/2, Lorenz-Mie theory generally consumes much more time than GOA as evidenced in Fig. 6(b). As r increases, the runtime ratio of Lorenz-Mie theory to GOA t Mie /t GOA grows   linearly with respect to r, and can easily reach two orders of magnitude difference in performance. This is explained by the fact that the number of terms in Lorenz-Mie theory is linearly proportional to the size parameter α, as mentioned previously. In comparison, the runtime for GOA is independent of the radius r.
Considering the trade-off between accuracy and time complexity, we choose GOA when r 2 μm and switch to Lorenz-Mie theory otherwise. This makes the runtime of computing S 1 (θ) (or S 2 (θ)) almost constant with respect to r while retaining an accurate method when required.

Bulk optical properties with graininess
Next, we consider light scattering by a cloud of spherical particles of the same composition but different sizes. The particles are assumed to be in each other's far-field regimes and their sizes are likely to range from wavelength-scale to much larger. In this section, we first discuss the particle size distribution that may vary spatially and then study the bulk optical properties of the discrete participating medium considering graininess. Thanks to the high efficiency of GOA, we are able to evaluate bulk optical properties on-the-fly.

Particle size distribution
We use the particle size distribution (PSD) N (r) to describe the population of particles in a discrete participating medium. N (r) dr is the total concentration (number of particles per unit volume) of particles with sizes in the range [r, r+dr]. The total particle concentration within some limited interval [r min , r max ] of sizes is obtained by N 0 = r max r min N (r) dr. To use N (r) as a probability density distribution (PDF), we have to normalize N (r) via N (r)/N 0 .
It is generally reported that particle sizes closely follow a log-normal distribution [17,68]: in which σ g is the geometric standard deviation andr g is the geometric mean radius. Obviously, this statistical tendency stems from observation of a large number of particles. A small number of particles will give rise to a size distribution deviating from the log-normal distribution. To demonstrate this in 2D, we generated 2000 particles in a box with uniformly distributed positions and log-normal distributed radii. A visualization of these particles and their PSD are shown in Fig. 7(above). We then extracted three small patches from the box and estimate their actual PSDs by binning. As expected, the PSDs plotted in Fig. 7(below) vary spatially and contain quite different features. In what follows, we use the notation N (r, x) to emphasize that the PSD varies spatially. Nevertheless, the ensemble average over these spatially-varying PSDs converges to the log-normal distribution shown in Fig. 7(above). Rendering with this global PSD yields a smooth appearance similar to that from a traditional continuous medium.

Bulk optical properties
With the spatially-varying PSD N (r, x), we are able to obtain the bulk optical properties of a local area in which many independent particles are immersed. Supposing that V x is a small volume centered around x, the bulk extinction coefficient x σ t of this volume is evaluated by (16) x In rendering literature, the symbol σ refers to the cross section sometimes, while using μ for the coefficient.
in which μ(V x ) is the measurement of V x , and r min (V x ) and r max (V x ) are the minimum and maximum particle radii inside V x , respectively. For brevity, we simplify both r min (V x ) and r max (V x ) by dropping V x henceforth. The scattering coefficient and the absorption coefficient can be defined in a similar way by replacing C t (r) with C s (r) and C a (r), respectively y . Generally, these properties exhibit multi-scale effects with respect to the size of V x .
The ensemble phase function is derived as in which σ s (V x ) serves as the normalization factor for f (V x , θ). Figure 8 visualizes the phase functions generated by the above formula. Here, we use uniform sampled radii between r min and r max (which is similar to the log-normal distribution with a very large σ g ). By fixing other properties, we show the influence of these phase functions on the final appearance of a smooth homogeneous medium in Fig. 8. Clearly, this formula is only valid for σ s (V x ) = 0. When σ s (V x ) = 0, i.e., the volume V x is free of particles, f (V x , θ) degenerates to a delta function: Similarly, we can derive the transmittance along a light beam of length s as T A (s) (18) in which C = A×s represents a small cylinder around the light beam with cross section A.
For particles with a monodisperse distribution, the transmittance is simplified to (19) in which C t is a constant and N (x) is a function of x only. The integral above formula simply counts the number of particles in the query region A × s.

Rendering solution
With these bulk optical properties, we are able to derive a multi-scale volumetric rendering equation (VRE) describing radiative transfer in discrete random media. We also develop a Monte Carlo sampling based method to solve the VRE; it only requires the position and size of each particle, avoiding explicit tessellation of its shape. In a preprocessing stage, we generate and store N tot particles with random positions and log-normal distributed radii for a discrete participating medium. During rendering, the stored particles are queried to determine the optical properties for each traced ray. A uniform grid is used for acceleration.

Multi-scale VRE
Conventionally, the VRE describing macroscopic light scattering in participating media is written as (20) in which L(x, ω) is the radiance arriving at x ∈ R 3 along a direction ω ∈ S 2 , f represents the scattering phase function characterizing the probability of radiation incident from ω being scattered into direction ω, z is the distance through the medium to the nearest boundary at z = x − zω, and y = x − yω is a point at distance y ∈ (0, z). The conventional transmittance between x and y is computed as T (x, y) = exp(− y 0 σ t (x − sω) ds). By substituting the multi-scale properties into the above equation, we arrive at a multi-scale version of the VRE: Since this multi-scale VRE is a general extension to the conventional one, it naturally supports multiple scattering.

Query cylinder
In our multi-scale VRE, every optical property depends on a PSD while the PSD is defined on a differential volume. To evaluate the transmittance T A (x, y) between any two positions x and y in the medium, we need a differential volume around the ray x → y. This volume is used to query particles which contribute to the transmittance T A (x, y). In our implementation, this is a thin cylinder centered around x → y, as illustrated in Fig. 9. We call such a cylinder a query cylinder. In this sense, we view each ray as a "fat ray" which gathers small particles along its trajectory. This is quite different from the implementation of rendering continuous media in which the optical properties are determined globally, without explicitly querying particles in a local area. In theory, the cross section A should be infinitely small. However, making it too small one may cause numerical issues and large variance. Conversely, bias will be introduced in T A (x, y) if A is too large. In practice, we select A as follows and keep it unchanged as the ray traverses the medium. Supposing that z near and z med respectively denote the depth of the near plane and the smallest depth of the medium in the view frustum, the size of A is selected according to Here, S pix is the pixel's size and k can be viewed as the proportion of the pixel's footprint at distance z med . Typically, satisfactory results are obtained when k is in the range [0. 5,1]. The influence of k on the visual effect is discussed in the next section.

Gathering particles
Gathering particles within a query cylinder (central ray: x → y, radius: r c and cross section: A) requires conducting sphere-cylinder intersection tests for every particle in the medium. A particle with position p and radius r p is inside the query cylinder if the distance from p to the central ray x → y is smaller than r p + r c , as illustrated in Fig. 9(left). Testing all particles of the medium is notoriously time-consuming. To boost the performance, we accelerate the process of ray traversal in the medium using a 3D digital differential analyzer (3D-DDA) [73,74]. Specifically, we construct a uniform grid for the medium and adopt a ray traversal algorithm, similar to that in Ref. [73], to find the active voxels intersected by the query cylinder. Figure 9(right) illustrates all active voxels corresponding to the orange query cylinder. Only those particles inside the active voxels are tested against the query cylinder. We determine the active voxels simply by the ray x → y. This introduces negligible bias, because the radius of the cross section is over two orders of magnitude smaller than the side length of the voxel. This is significantly different to the thick beams used in beam radiance estimation [25]. After collecting all particles inside the query cylinder, we accumulate their contributions to the transmittance T A (x, y) according to Eq. (18).
To construct a uniform grid, its resolution should be carefully determined. We have observed by experiments that high performance is achieved when roughly one particle resides in each voxel after spatial subdivision.

Computing Q(θ)
To solve the multi-scale VRE, we also have to compute Q(θ) at each sample position y. Q(θ) describes the angular distribution of scattering at y. To quickly compute it, a small query region around y should be defined. We set this query region to a small sphere with radius r c . If this query region contains N q > 0 particles, we evaluate Q(θ) using the following formula:

Importance sampling
Following the simplifications used in rendering surface glints [52,53], importance sampling is performed according to the global optical properties of the medium, assuming it to be continuous. Specifically, we use the global extinction coefficient for free-flight sampling and use the tabulated global phase function for angular sampling. These global optical properties only need to be determined once during preprocessing, assuming the entire bounding box of the medium to be V x in Eqs. (16) and (17).

Setting
We have implemented our rendering solution on top of the Mitsuba renderer [75], with spectral rendering enabled. We use 8 spectral samples in the range of the visible spectrum at equally-spaced locations [46]. After rendering, we convert the spectral values to the sRGB color space. All synthesized images were created on a PC with an Intel 16-core i7-6900K CPU and 16 GB RAM.
To compute the bulk optical properties of a discrete participating medium, we need to specify the complex refractive index (η) for the particles involved and the global PSD (r g and σ g ). The physical unit for the particle radius is μm. We also provide an upper bound to the particle size (r max = 2000) to avoid unreasonably large particles which are unusual and unsuitable for treatment as participating media. As mentioned previously, we use GOA in the calculation of the scattering amplitude functions when r 2 and switch to Lorenz-Mie theory otherwise. Except for the Staircase scene, the refractive index is set using data for ice from Ref. [17], and the host medium is set to air with η m = 1.

Comparison to explicit path tracing
We first compare our method to traditional path tracing. Previous methods simulating the grainy appearance of discrete participating media mostly rely on explicit path tracing (EPT), possibly with approximations to simplify the computation of highorder scattering [12][13][14]. However, EPT and other approximations are restricted to geometric optics. This means that only surface reflection and refraction are properly handled. In principle, a mesh should be explicitly generated for every particle and costly ray-object intersections are required.
Compared to EPT, our method offers at least two benefits. First, particle scattering is considered which includes Fraunhofer diffraction and phase differences, etc. This significantly expands the range of particles that can be handled. Notably, these optical phenomena are quite important in Fig. 10 Comparisons between EPT and our method in simulating the grainy appearance of discrete participating media for different PSDs (Ntot = 5 × 10 5 ).
correctly simulating light scattering, especially for very small particles. To verify this, we produced and rendered 5 × 10 5 spherical particles with different PSDs in Fig. 10. These particles have random positions and log-normal sampled radii. For EPT, small transparent balls are instantiated in the scene. Since only surface reflection and refraction are computed for EPT (above), the energy inherently belonging to Fraunhofer diffraction is not correctly captured by EPT, resulting in overly sparse and specular volumetric glints in this Lamp scene. The importance of Fraunhofer diffraction is shown in Fig. 11, which plots the percentage of energy contributed by Fraunhofer diffraction at different scattering angles and for different sized particles. Obviously, Fraunhofer diffraction cannot be ignored, especially for small particles. Moreover, EPT easily misses many small particles that are hard gather along an ordinary light path, leading to a slow convergence rate. Also, even when substantially increasing the sampling rate or using "fat ray" tracing as in our method, the appearance is still quite sparse. On the contrary, our method (below) preserves the energy from Fraunhofer diffraction and produces a smoother appearance that is closer to reality, thanks to the Airy pattern x caused by Fraunhofer diffraction and the efficient rendering solution tailored for sparse media.
When only very large particles exist in the medium, our method and EPT tend to produce a similar grainy appearance, as shown in Fig. 12. Here, we render a discrete medium with 10 5 particles of the same size (1000 μm). Since the radius is sufficiently large, pure geometric optics becomes applicable and serves as a valid approximation to particle scattering. This is further verified in Fig. 11(c). We see that the contribution of Fraunhofer diffraction is concentrated in a very narrow angle in this case, making it hard to observe. Consequently, if we remove Fraunhofer diffraction from our model and leave only reflection and refraction (see the third column of Fig. 12), we achieve a grainy appearance similar to that of our full model. However, there are still subtle differences (highlighted in the difference image) contributed by Fraunhofer diffraction.
The second strength of our rendering solution lies in its runtime performance. For rendering these sparse media in the Lamp scene and the Cornell Box scene in Figs. 10 and 12, our method achieves a 1.2× speedup over EPT in rendering time.

Comparison to the Henyey and Greenstein model
In our framework, we derive the phase function from Lorenz-Mie theory and GOA. In computer graphics, it is more common to adopt an empirical model, e.g., the Henyey and Greenstein (HG) model [76], because of its simplicity and well-defined behavior. However, the HG model is not very accurate, as pointed out by various works [69,77,78]. In Fig. 13 we compare our estimated phase function with the best-fit HG phase function when rendering the same scene as in Fig. 12. Since the HG phase function cannot faithfully encode the scattering pattern from these relatively large particles (r = 1000 μm), the rendering from the best-fit HG phase function is slightly different from ours (see insets). Recall that our rendering is close to that generated using EPT for this specific scene with large particles.

Comparison to continuous media
We also compared the grainy appearance of discrete media simulated by our method with the smooth appearance of continuous media. In Fig. 14, we Fig. 13 Comparisons between our directional scattering model and the best-fit Henyey and Greenstein (HG) model [76]. The scattering model generated by Lorenz-Mie theory (green curve, right) is also provided for reference. Due to obvious discrepancies between our estimated phase function (red curve, right) and the HG phase function (blue curve, right), the renderings differ slightly (see the close-ups).

Fig. 14
Visual comparisons with continuous media for varying optical properties. As the particle number Ntot increases, the appearance of the discrete medium exhibits obvious low-frequency haziness and converges to that of a continuous medium.
assume that the particles in the Dragon scene possess the same radius r and are randomly distributed in a cube of volume V . In this configuration, it is easy to derive the extinction coefficient and the scattering coefficient for the continuous media as C t (r)N tot /V and C s (r)N tot /V , respectively. The phase function can be computed in a similar way and stored in a table. With these global properties, traditional volumetric path tracing is applicable to render these continuous media. Generally, each discrete medium and its paired continuous medium have the same overall brightness. For the discrete media, when the number of particles N tot is small, we can clearly observe individual particles lit by the lamps. As N tot increases, the appearance becomes hazy, and the rendering result becomes closer to that of the corresponding continuous medium. When N tot is sufficiently large (e.g., N tot = 10 7 ), the discrete medium and its continuous counterpart provide quite similar appearance. In Fig. 15, a much denser medium (N tot = 10 8 ) is rendered by our method, which again produces a smooth appearance matching that from a continuous medium. These pair-wise comparisons demonstrate that our rendering solution converges to traditional volumetric rendering of continuous media in the limit of particle concentration.

Query cylinder cross section
We determine the query cylinder's radius r c according to Eq. (24) in which the parameter k plays an important role. We suggest choosing its value in the range [0.5, 1] which yields a reasonable grainy appearance as shown Fig. 16(b). Values in this range allow us to faithfully capture almost all grains in scenes with

Fig. 15
Comparison between our rendering solution and its continuous counterpart for a very dense medium (r = 400 μm and Ntot = 10 8 ).
little bias. Generally, setting k too large produces undesirable aliasing as shown in Figs. 16(c) and 16(d). In these two cases, although the overall brightness is similar to that for k = 0.5, the volumetric glints are overly blurred. On the other hand, a too small k will miss many particles during querying and is therefore inefficient, especially when the medium is very sparse, e.g., N tot = 10 4 . However, k with a value smaller than 0.5 may sometimes be acceptable. For instance, the first image in the bottom row of Fig. 16  is rendered with k = 0.2 which achieves a similar effect to that for k = 0.5 for this relatively dense medium (N tot = 10 5 ).

Choice of grid resolution
We employ a uniform grid to accelerate the ray traversal process, considering that particles are uniformly distributed in the scene. As shown Fig. 17(right), the grid resolution (Res) influences the performance. Although no analytical analysis can be provided to select the best resolution, we empirically observe that a grid resolution yielding roughly one particle per voxel achieves the optimum solution for most scenes. The Bench scene in Fig. 17(left) contains 2 × 10 5 particles and best performance is achieved at a resolution of 60 × 60 × 60. As the resolution increases, the runtime grows steadily due to the additional cost introduced by the grid. However, a resolution much lower than 60 × 60 × 60 has very poor performance since too many particles reside in each voxel. For other scenes, a similar conclusion can be drawn.

Impact of the refractive index
The proposed rendering solution can be easily generalized to support absorbing particles. The absorption cross section of each particle is computed by Eq. (13) which varies linearly with the imaginary part of the complex refractive index, η i . The extinction cross section must be slightly modified according to Ref. [65]. Results of doing so are shown in Fig. 18. This Staircase scene aims to simulate floating dust in a dirty room lit by a local area light through the window. As expected, increasing η i increasingly reduces the intensity of scattering. Figure 19 analyzes the impact of the global PSD on the appearance of discrete participating media. Here, we generated N tot = 4 × 10 6 particles with different global PSDs in the House scene. This scene is designed to simulate the appearance of blowing snow in the sky. The first row shows the impact of the geometric mean radiusr g . The general trend is that the scattering effects become increasingly prominent asr g grows, since the scattering coefficient is positively correlated with particle radius. The geometric standard deviation σ g is responsible for the level of graininess as shown in the second row of Fig. 19. A small σ g tends to generate smoother appearance than a large one. This is to be expected since a large σ g means local PSDs changing widely, leading to stronger graininess.

Impact of global PSD
In fact, the PSD is not limited to a log-normal distribution. Other distributions also work. For instance, in the third row of Fig. 19, we show a cloud of blowing snow following a bimodal log-normal distribution. We generated 1.99 × 10 7 small particles withr g,1 = 100 (orr g,1 = 200) and σ g,1 = 1, and also generated 10 5 large particles withr g,2 = 1000, σ g,2 = 1.2, leading to 2 × 10 7 particles in total. Using this complex distribution, we can observe both haziness from massive small particles and graininess from a handful of large particles.

Performance
The runtime performance for some test scenes can be found in Table 2. As mentioned previously, our rendering solution achieves roughly 1.2× the speed of EPT for the Lamp scene and the Cornell Box scene. As the number of particles N tot increases, the improvement becomes more evident. As shown in Fig. 20, twice the speed is achieved when N tot is increased to 10 6 for the Cornell Box scene. Since only the position and the radius of each particle are  required, the memory consumption of our rendering solution is affordable even when N tot is very large. The storage scales with the number of particles and the resolution of the grid.

Limitations and future work
Although our framework successfully simulates the grainy appearance of discrete participating media, it has several limitations deserving further research.

Aspherical particles
Our current framework focuses on discrete participating media composed of spherical particles. However, aspherical particles are also very common, e.g., in the context of rendering rainbows [68] or large snow flakes. Extending our framework to aspherical particles requires the derivation of new expressions for the scattering amplitude functions S 1 and S 2 . For some special particle shapes, determining S 1 and S 2 is rather straightforward and analytical expressions exist [61,66,67]. For more general shapes, precomputation would be required in practice.

Spatial correlations
As we assume the particles to be sparsely distributed, spatial correlations between particles are not considered and we directly extend the conventional VRE to support multi-scale graininess. However, for densely packed particles the conventional VRE becomes questionable due to strongly correlated scattering effects [12][13][14]. Recently, new radiative transfer frameworks dedicated for correlated media have become available in computer graphics [5,6,9]. It would be interesting future work to investigate a more general framework supporting both effects.

Mixed particles
Another possible direction for future work is to explore an efficient strategy to handle particle mixtures with different compositions. For instance, dust in the real world may be made up of soil particles, textile fibers, human skin cells, etc. A physically-correct participating medium should consider such heterogeneous granular mixtures. For continuous media, this is relatively simple since the concentrations of different grains are fixed [17]. However, for discrete participating media, the concentrations are dynamic and change spatially as for PSDs [14]. Determining the bulk optical properties should take spatially-varying concentrations into consideration.

Conclusions
We have developed a general physically-based framework for modeling and rendering discrete participating media composed of massive assemblies of independent particles. Notable characteristics of these media include a wide range of PSDs and the appearance of graininess. To faithfully simulate their appearances, we have derived a novel multiscale VRE in which a combination of Lorenz-Mie theory and GOA is used to enable high-efficiency evaluation of the important optical properties. A Monte Carlo rendering solution is used to solve the multi-scale VRE with high accuracy and low computational cost. We have extensively evaluated our framework and compared it to conventional methods, demonstrating that the proposed framework allows us to reproduce a variety of grainy appearances stemming from different discrete participating media. We can guarantee temporal coherence in animation. We have greatly extended the participating media rendering framework to handle a much larger range of particle sizes. We believe that the proposed framework is a further step in computer graphics to manage the details of participating media and expect further insightful exploration of this phenomenon in future research.
Ling-Qi Yan is an assistant professor of computer science at UC Santa Barbara, co-director of the MIRAGE Lab, and affiliated faculty in the Four Eyes Lab. Before that, he received his doctoral degree from the Department of Electrical Engineering and Computer Sciences at UC Berkeley and his bachelor degree in computer science from Tsinghua University. His research interests include physically-based rendering, real-time ray tracing, and realistic appearance modeling.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Other papers from this open access journal are available free of charge from http://www.springer.com/journal/41095. To submit a manuscript, please go to https://www. editorialmanager.com/cvmj.