Rendering Discrete Participating Media with Geometrical Optics Approximation
Jie Guo, Bingyang Hu, Yanjun Chen, Yuanqi Li, Yanwen Guo, Ling-Qi Yan
IIEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 1
Rendering Discrete Participating Media withGeometrical Optics Approximation
Jie Guo, Bingyang Hu, Yanjun Chen, Yuanqi Li, Yanwen Guo and Ling-Qi Yan
Abstract —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 the scale several orders of magnitude greater than thewavelength, and the appearance shows distinct graininess as opposed to the smooth appearance of continuous media. Onefundamental issue in physically-based synthesizing this appearance is to determine necessary optical properties in every local region.Since these optical properties vary spatially, we resort to geometrical optics approximation (GOA), a highly efficient alternative torigorous Lorenz-Mie theory, to quantitatively represent the scattering of a single particle. This enables us to quickly compute bulkoptical properties according to any particle size distribution. Then, we propose a practical Monte Carlo rendering solution to solve thetransfer of energy in discrete participating media. Results show that for the first time our proposed framework can simulate a widerange of discrete participating media with different levels of graininess and converges to continuous media as the particle concentrationincreases.
Index Terms —Light scattering, Geometrical optics approximation, Discrete participating media, Volume rendering. (cid:70)
NTRODUCTION
Rendering participating media is a long-standing prob-lem in computer graphics, with much effort devoted tosolving this problem plausibly and efficiently [1], [2]. Fromthe physics point of view, radiative transfer is rather com-plicated and should be rigorously derived from Maxwell’selectromagnetic theory. To make the simulation tractable,some compromises are made in physically-based renderingover the past decades. Two main assumptions are indepen-dent scattering [3] and local continuity (or statistical homo-geneity [4]). The first assumption of independent scatteringmeans that the particles forming the medium are far apartfrom each other and are mutually unaffected. Recently,this assumption has been relaxed in computer graphics viaincorporating spatial correlations between scatterers intoradiative transfer frameworks [5], [6], [7], [8], [9], leadingto non-exponential attenuation of light.The second assumption of local continuity implies thatthe medium is homogeneous and compact in each dif-ferential volume even if macroscopic heterogeneity exists.Under this circumstance, the light is not sensitive to thediscrete spatial distribution of the scatterers, but only totheir local average properties. Consequently, light scatteringphenomena take place at any point of the medium, resultingin locally smooth renderings. However, many participatingmedia are composed of separate particles distributed ran-domly within a given volume. The scattering will happenonly at the particle positions [4], [10]. These facts indicatethat the assumption of continuous media only holds when • Jie Guo, Bingyang Hu, Yanjun Chen, Yuanqi Li and Yanwen Guo are withthe State Key Lab for Novel Software Technology, Nanjing University,Nanjing, Jiangsu 210023, P.R. China.E-mail: [email protected], [email protected], [email protected],[email protected] and [email protected] • Ling-Qi Yan is with Department of Computer Science, UC Santa Barbara,United States.E-mail: [email protected] the particle size is much smaller as compared with theresolution of the sensor (e.g., human eyes) and the quan-tity is sufficiently large [11]. Otherwise, individual grainscan be observed when zoomed in. To break through thisconstraint, the graininess of the medium should be takeninto consideration.Several recent studies in computer graphics [12], [13],[14] have noticed the graininess in rendering granularmaterials. They typically rely on explicit geometries andprecomputed transport functions to capture the appearanceof discernible grains. As these approaches are designed fora certain amount of very large particles under geometricoptics, the scattering behaviors of particles forming themedia are rather limited. For instance, diffraction, whichdominates light scattering from small particles, is oftenignored. In this paper, we attempt to put forward a moregeneral approach to model and render discrete participatingmedia with a large amount of scatterers whose particlesize distributions (PSDs) range widely. These media areomnipresent in natural and artificial environments, such asflying dusts, blowing snows, powder suspension, and airbubbles in liquid.Unlike conventional continuous media, these discretemedia have spatially-varying optical properties (e.g., theextinction coefficient) that cannot be determined in advancebut should be evaluated on-the-fly. They are closely relatedto the scattering behavior of each individual particle andthe fluctuation of PSDs. To derive necessary optical proper-ties for any local region of a given discrete participatingmedium in a physically-based manner, one can resort toLorenz-Mie theory [15], [16] which offers substantial realismin rendering participating media [17]. However, numericalevaluation of the Lorenz-Mie coefficients is known to bedifficult and time-consuming when the particle size be-comes large [18]. To ameliorate this issue, we introducegeometrical optics approximation (GOA) [18], [19] and use a r X i v : . [ c s . G R ] F e b EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 2 it to simplify the computation of light scattering when theparticle is sufficiently large. We show how to perform asmooth transition between Lorenz-Mie theory and GOAin computing the optical properties, enabling both highaccuracy and low computational cost.Due to the variations in local PSDs, the optical propertiesexhibit multi-scale effects with respect to the scene config-uration. We derive a novel multi-scale volumetric renderingequation (VRE) and propose a practical Monte Carlo render-ing solution to solve it. Our solution only relies on the po-sition and the radius of each particle distributed randomlyaccording to some PSDs, avoiding cumbersome geometricmodeling and lengthy precomputation. Experimental re-sults verify that the proposed solution is able to capture thedistinct grainy appearance of discrete participating mediaand guarantee temporal coherence in animation. We alsoshow that it converges to continuous media in the limit ofparticle concentration.In summary, the main contributions of this paper are: • a general and physically-based framework for mod-eling and rendering discrete participating media,considering diffraction, polarization, a wide range ofPSDs, etc, • the use of GOA for efficient and accurate evaluationof the multi-scale bulk optical properties at any localregion of a medium, and • a new Monte Carlo rendering solution that capturesboth low-frequency haziness and high-frequencygraininess in discrete participating media. ELATED W ORK
Rendering participating media is a challenging but impor-tant 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 ourcurrent framework, we choose Monte Carlo path integrationby virtue of its elegant simplicity, generality, and accuracy.This technique operates by stochastically constructing alarge number of light paths between sensors and emittersto simulate the light transport in the scene. To facilitate thequery of particles along paths, we augment each ray witha cylinder, in a way similar to the photon beam [25], [26],[31]. Multiple importance sampling [32], [33] algorithms arebeneficial for reducing the large variance caused by MonteCarlo sampling.
Since the original VRE is only a rough approximation to thereal radiative transport in participating media, the range ofappearance that can be faithfully simulated is limited. Re-cent trend in computer graphics tries to capture more detailsin the volume by relaxing the assumptions in the originalVRE, enriching the range of achievable appearances. Forexample, to account for angular anisotropy, the VRE isextended with local directional dependency based on themicroflake model [34], [35], [36]. It is also possible to extendthe 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 statisticallyhomogeneous within each differential volume [4], ignoringany sub-pixel details.To handle participating media with complex 3D struc-tures, volumetric representations of explicit geometries havebeen widely used. By capturing the geometric and opti-cal properties of a fabric down to the fiber level, micro-appearance models, described using high-resolution vol-umes, offer state-of-the-art renderings for fabrics and tex-tiles [37], [38], [39], [40], [41]. Unfortunately, these methodsare highly data-intensive and plagued by heavy computa-tion. To improve the performance while maintaining goodaccuracy, some downsampling strategies [42], [43] are devel-oped. Current rendering solutions for granular materials arealso based on explicit geometries and pre-captured opticalproperties of each individual grain [12], [13], [14]. Theygenerally employ shell tracing to make large jumps insidemedia. Even so, the computational cost is still high. The goalof this work is to develop a general framework for handlingparticipating media with graininess which also allows rapidcomputation and convenient usage.
Our work is also closely related to the simulation of glints onsurfaces. Yan et al. [44], [45] suggested using explicit high-resolution normal maps to model sub-pixel surface detailsand successfully simulated spatially-varying glints with apatch-based normal distribution function. Subsequent workadopted a wave optics model to achieve more accurateresults with noticeable color effects [46]. There are alsoother methods focusing specifically on capturing spatially-varying highlights from scratched surfaces, under eithergeometric optics [47], [48], [49] or wave optics [50]. To easethe burden of computation and storage, Kuznetsov et al.[51] proposed to learn high-frequency angular patterns fromexisting examples, using a generative adversarial network(GAN). Jakob et al. [52] addressed the problem of glit-tery surface simulation using a purely procedural approachwhich requires far less storage and supports on-the-fly pointqueries. This approach has been extended to incorporateiridescence [53] and allow fast global illumination [54].
Lorenz-Mie theory [15], [16] develops a rigorous solution tothe problem of light scattering by spherical particles. It wasintroduced to the graphics community by Rushmeier [55] toaccurately simulate the physics of light transport in partic-ipating media. Later, Callet [56] used this theory to modelpigmented materials consisting of pigmented particles in atransparent solvent. Atmospheric phenomena, such as halosand rainbows, are especially favored by this theory [57],[58], [59], [60]. Frisvad et al. [17] generalized the originalLorenz-Mie theory and used it to compute the appearanceof materials with different mixed particle concentrations.Though accurate, this theory is computationally expensive.We show that GOA [18], [19], [61], [62], [63], [64], [65] ismuch more efficient than Lorenz-Mie theory in computingoptical properties of individual particles in various media,especially when the particle size is large. Moreover, the
EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 3 computation of GOA can be made in non-ideal situationssuch as absorbing particles [65] and non-spherical particles[61], [66], [67], [68]. Compared to Lorenz-Mie theory, GOAis less explored in computer graphics. We choose GOA inour framework, taking advantage of its high performance.
IGHT S CATTERING BY
A S
INGLE P ARTICLE
We first study light scattering by a single particle. Wesuppose that the particle is approximately spherical andhas a set of physical properties including its radius r andthe refractive index η p . Currently, we assume that particlesforming the medium have the same composition and onlytheir 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 particleas η = η p /η m . The size of a spherical particle may alsobe expressed in terms of the dimensionless size parameter α = kr = 2 πη m r/λ , where k is the wave number defined by k = 2 πη m /λ and λ is the wavelength of light in the medium.To describe the scattering, we need two scattering am-plitude functions: S ( θ, ϕ, r ) and S ( θ, ϕ, r ) , where θ is thescattering angle and ϕ is the azimuth angle. The subscripts and denote perpendicular and parallel polarizations,respectively. For spherical particles, S and S are invariantwith respect to ϕ , but they change depending on the radius r . For unpolarized light, these two functions define thephase function of a single particle as [4] f p ( θ, r ) = | S ( θ, r ) | + | S ( θ, r ) | | k | C s ( r ) (1)which is properly normalized by the scattering cross section C s : C s ( r ) = (cid:90) π (cid:90) π | S ( θ, r ) | + | S ( θ, r ) | | k | sin θ d θ d φ. (2)Another important property of the particle is the extinctioncross section C t which is evaluated by C t ( r ) = 4 π Re (cid:26) S (0 , r ) | k | (cid:27) (3)with S (0 , r ) = S (0 , r ) = S (0 , r ) . The notation Re takes thereal 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 ( θ, r ) and S ( θ, r ) are available, we can easily find the scattering,extinction and absorption cross sections as well as the phasefunction of the particle. For light scattering of an elec-tromagnetic wave from a homogeneous spherical particle,exact solutions of the two scattering amplitude functions aregiven by Lorenz-Mie theory [15], [16]. Its accuracy has beenvalidated against real measurements in various literature[69], [70]. Please refer to Appendix A for more details.As a rigorous and general electromagnetic treatment oflight scattering by spherical particles, Lorenz-Mie theorycan precisely handle a wide range of particle sizes. How-ever, as the particle size increases, numerical calculations ofthe Lorenz-Mie coefficients become very tedious and time-consuming, due to the fact that the number of terms to be Fig. 1. Light scattering by a spherical particle in the GOA picture. computed in the series for S ( θ, r ) and S ( θ, r ) is propor-tional to the size parameter α [17], [18]. For this reason,simpler approximate expressions should be developed toreduce the computational complexity. In the case that theparticle size is large with respect to the wavelength of theilluminating light, geometrical optics approximation (GOA)[18], [19], [61], [62], [63], [64], [65] provides a simplified butalso good solution. Within the framework of GOA, light scattering is calculat-ed by a superposition of classical diffraction, geometricalreflection and transmission. The diffraction is independentof the particle’s composition (i.e., the refractive index). Itsamplitude functions for the forward direction are readilydescribed by the Fraunhofer diffraction as [18] S D , ( θ, r ) = S D , ( θ, r ) = α J ( α sin θ ) α sin θ (4)where J is the first-order Bessel function.Leaving out diffraction, a light ray hitting a sphericalparticle at an incident angle θ i is partially reflected andpartially refracted depending on the properties of the in-terface, as sketched in Fig. 1. The refracted ray may undergoa number of internal reflections before leaving the particle.For each emerging ray, we use an integer p to denote thenumber of chords it makes inside the spherical particle.Obviously, the externally reflected ray has p = 0 while theother rays are transmitted with p − internal reflections. Theangle of deflection θ p between the p th emerging ray and thedirection of the incident ray is given by θ p = 2 pθ t − θ i − ( p − π (5)with sin θ i = η sin θ t according to Snell’s law. The scatteringangle θ is further determined by the deflection angle θ p as θ = q ( θ p − πl ) (6)where q ∈ { , − } and l is an integer ensuring that thescattering angle θ is well defined in the range between and π .
1. It is suggested that an appropriate number of terms to sum is (cid:100)| α | +4 . | α | / + 1 (cid:101) [71].2. q = 1 indicates that the incident ray hits the particle on the upperhemisphere and q = − for the lower. EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 4
Clearly, reflected and transmitted rays depend on theshape and composition of the particle. Their scattering am-plitudes for each polarization are derived as [10], [19] S ( p ) j ( θ, r ) = α(cid:15) j ( θ i ) (cid:115) sin 2 θ i θ | d θ p / d θ i | e i φ j = 1 , . (7)Here, the fraction (cid:15) j ( θ i ) , which is due to the reflectionand/or refraction for an emergent ray of order p , is definedas (cid:15) j ( θ i ) = (cid:26) R j ( θ i ) p = 0(1 − R j ( θ i ) )( − R j ( θ i )) p − p > (8)with R j ( θ i ) being the Fresnel reflection coefficients. Thephase difference φ = φ p + φ f includes φ p due to the lengthof optical path: φ p = 2 α (cos θ i − pη cos θ t ) (9)and φ f due to focal line: φ f = π (cid:16) p − l − s − q (cid:17) (10)with s = sgn(d θ p / d θ i ) = sgn(2 p tan θ t / tan θ i − .Putting together S D ,j ( θ, r ) and S ( p ) j ( θ, r ) , we are able toget the total amplitude functions of GOA as S j ( θ, r ) = (cid:40) (cid:80) ∞ p =0 S ( p ) j ( θ, r ) + S D ,j ( θ, r ) θ ∈ [0 , π ) (cid:80) ∞ p =0 S ( p ) j ( θ, r ) θ ∈ [ π , π ] (11)with j = 1 , . These expressions can be evaluated quiteefficiently.In GOA, analytical expression of C t ( r ) can be derived as(see the derivation in Appendix B) C t ( r ) = 2 πr + 2 πr | k | (cid:88) p ∈P (cid:15) j (0) | p/η − | cos( φ p + φ f ) (12)where P = { , , , · · · } . This expression is fast to evaluateand well captures the ripple structures [72].For absorbing particles, we certainly have C a ( r ) = C t ( r ) − C s ( r ) > . Within GOA, the absorption cross section C a ( r ) is faithfully approximated by [4] C a ( r ) = 16 π r η i λη r (cid:104) η − ( η − (cid:105) (13)in which η r and η i are the real and imaginary parts of η , re-spectively. The scattering amplitude functions S ( θ, r ) and S ( θ, r ) are also slightly different. The details are providedin Appendix C. p There is also an infinite summation in computing S and S of GOA. However, unlike that in Lorenz-Mie theory, thenumber of terms needed is independent of the particle size,and a small p suffices in most cases. As shown in Fig. 22,when calculating ( | S | + | S | ) / with p = 3 in GOA, we getan almost identical curve with that of p = 100 , irrespectiveof the particle size. This is because higher-order reflections( p > ) carry much less energy and have very little impacton the scattered light intensities. Regarding this, p is safelyset to 3 in what follows. Please see more discussions inAppendix E. p=100p=3 (a) r = 1 µ m p=100p=3 (b) r = 100 µ m Fig. 2. Impact of p over ( | S | + | S | ) / in GOA. Here, η = 1 . and λ = 0 . µ m . -1 Lorenz-MieGOA (p=3)GOA (p=1) Lorenz-MieGOA (p=3)GOA (p=1)
Fig. 3. Variation of the extinction cross section C t as a function of r .Here, we compare our GOA (red: p = 3 , green: p = 1 ) calculationagainst that of Lorenz-Mie theory (blue) with η = 1 . and λ = 0 . µ m .The ripple structures [72] that are well captured by GOA are clearlyshown in the left diagram. Moreover, we can further simplify the extinction crosssection to C t ( r ) = 2 πr + 4 rλη ( η + 1) | η − | sin (cid:18) πrλ (1 − η ) (cid:19) (14)by setting p = 1 , since C t only relies on the value of S (or S ) evaluated at θ = 0 , and the light rays with p > contribute little to the forward scattering. This is evidencedin Fig. 3 where the C t curves of p = 1 (green) and p = 3 (red)are virtually indistinguishable for a very wide range of r . InAppendix E, we show that the Relative Mean Squared Error(RelMSE) is less than − for η = 1 . and λ = 0 . µ m using p = 1 and p = 3 , respectively. To investigate the range of validity of GOA in simulatingthe scattering patterns of spherical particles, we comparethe results with the rigorous Lorenz-Mie results on a widerange of particle radii in Fig. 4. The wavelength of theincident rays is set to . µ m and the relative index ofrefraction is η = 1 . in all the calculations. Fig. 4 revealsthat the scattering amplitude distributions by GOA alignwell with those obtained with Lorenz-Mie theory for largeparticles with r > µ m . The agreement of these twomethods is especially good in almost all directions when theradius is large enough (e.g., r = 100 µ m ). However, when r ≈ µ m , some discrepancies between the two methodsappear. These discrepancies become large as the radius ofparticle decreases further.To show the influence of these discrepancies on the per-ception of the translucent appearance, we render a smoothmedium comprising monodisperse particles of radius r .
3. The RelMSE for C t is calculated by ( C p =3t − C p =1t ) / ( C p =3t ) EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 5
Lorenz-MieGOA log | S | log | S | (a) r = 0 . µ m Lorenz-MieGOA log | S | log | S | (b) r = 1 µ m Lorenz-MieGOA log | S | log | S | (c) r = 2 µ m Lorenz-MieGOA log | S | log | S | (d) r = 10 µ m Lorenz-MieGOA log | S | log | S | (e) r = 100 µ m Fig. 4. Visual comparisons of log | S | (upper half) and log | S | (lower half) by Lorenz-Mie calculations (blue curves) with those by GOA (red curves)for η = 1 . and λ = 0 . µ m . The particle radius is set to r = 0 . , , , and µ m , respectively. GO A L o r e n z - M i e (a) r = 1 µ m (b) r = 2 µ m Fig. 5. Rendering a smooth medium with phase functions derived fromLorenz-Mie theory (top row) and GOA (bottom row), respectively. Here, η = 1 . . The red arrows highlight the differences. -5 R e l M S E -4 -2 Fig. 6. Comparisons between GOA and Lorenz-Mie theory in RelMSEand runtime complexity. Left: Variation of RelMSE as a function of r .Right: Variation of runtime ratio t Mie /t GOA as a function of r . We determine the phase function using either Lorenz-Mietheory or GOA according to r . The extinction coefficient isset to a constant . for a fair comparison. The synthesizedimages are presented in Fig. 5 with r setting to µ m and µ m , respectively. Clearly, the differences of phasefunctions between Lorenz-Mie and GOA in the case of r = 1 µ m result in inconsistence of appearance. However,this inconsistence almost disappears completely when r goes up to µ m , although there are some mismatches on thebackward peaks of S ( θ ) (see Fig. 4(c), same for S ( θ ) ). Fig.4 and Fig. 5 together have verified the accuracy of choosingGOA to compute S ( θ ) and S ( θ ) when r ≥ µ m . Please
4. These phase functions are precomputed and stored in tables. see more comparisons and discussions in Appendix E.To further show the similarity between Lorenz-Mie theo-ry and GOA in computing S ( θ ) and S ( θ ) , we report theirRelative Mean Squared Error (RelMSE) in Fig. 6 left. TheRelMSE is computed on ( | S ( θ ) | + | S ( θ ) | ) / and differentrelative refractive indexes are tested. The results confirmthat subtle errors exist when r is large: despite some fluctu-ations, these calculations of GOA exhibit errors of less than . as compared with exact Lorenz-Mie calculations when r ≥ µ m .In the calculations of ( | S ( θ ) | + | S ( θ ) | ) / , Lorenz-Mietheory generally consumes much more time than GOA asevidenced in Fig. 6 right. As r increases, the runtime ratio ofLorenz-Mie theory and GOA t Mie /t GOA grows linearly withrespect to r , and can easily reach two orders of magnitudedifference in performance. This is explained by the factthat the number of terms in Lorenz-Mie theory is linearlyproportional to the size parameter α , as we mentioned pre-viously. In comparison, the runtime for GOA is independentof the radius r .Considering the trade-off between accuracy and timecomplexity, we choose GOA when r ≥ µ m and switchto Lorenz-Mie theory otherwise. This makes the runtime ofcomputing S ( θ ) (or S ( θ ) ) almost constant with respect to r while retaining the accuracy as much as possible. ULK O PTICAL P ROPERTIES WITH G RAININESS
Now, we consider light scattering by a cloud of sphericalparticles of the same composition but of different sizes.The particles are assumed to be in each other’s far-fieldregimes and their sizes are likely to range from wavelength-scale to the scale much larger than the wavelength. Inthis section, we first discuss the particle size distributionthat may vary spatially and then study the bulk opticalproperties of the discrete participating medium consideringgraininess. Thanks to the high efficiency of GOA, we areable to evaluate the bulk optical properties on-the-fly.
We use the particle size distribution (PSD) N ( r ) to describethe population of particles in a discrete participating medi-um. Thus N ( r )d r is the total concentration (particle numberper unit volume) of particles with sizes in the domain [ r, r + d r ] . The total particle number concentration withinsome limited interval [ r min , r max ] of sizes is obtained by EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 6
Fig. 7. Illustration of spatially-varying PSD. First row: A realization of polydisperse particles and its global PSD. Second row: ActualPSDs for particles in three small patches. N = (cid:82) r max r min N ( r )d r . To use N ( r ) as a probability densitydistribution (PDF), we have to normalize N ( r ) via N ( r ) /N .It is generally reported that particle sizes follow close toa log-normal distribution [17], [68]: N ( r ) = N √ πr ln σ g exp (cid:34) − (cid:18) ln r − ln ¯ r g ln σ g (cid:19) (cid:35) (15)in which σ g is the geometric standard deviation and ¯ r g is thegeometric mean radius. Obviously, this statistical tendencystems from the observation of a large number of particles. Asmall number of particles will give rise to a size distributiondeviating from the log-normal distribution. To demonstratethis in 2D, we generate particles in a box with uniform-ly distributed positions and log-normal distributed radii.The visualization of these particles and its PSD are shownin the first row of Fig. 7. We then extract three small patchesfrom the box and estimate their actual PSDs by binning.As expected, the PSDs plotted in the second row of Fig.7 vary spatially and contain quite different features. Inwhat follows, we use the notation N ( r, x ) to emphasizethat the PSD varies spatially. Nevertheless, the ensembleaverage over these spatially-varying PSDs converges to thelog-normal distribution shown in the first row of Fig. 7.Rendering with this global PSD yields a smooth appearancesimilar to that from a traditional continuous medium. With the spatially-varying PSD N ( r, x ) , we are able toobtain the bulk optical properties of a local area in whichmany independent particles are immersed. Supposing that V x is a small volume centered around x , the bulk extinctioncoefficient σ t of this volume is evaluated by σ t ( V x ) = 1 µ ( V x ) (cid:90) r max ( V x ) r min ( V x ) C t ( r ) (cid:90) x ∈V x N ( r, x )d x d r (16)in which µ ( V x ) is the measurement of V x . r min ( V x ) and r max ( V x ) return the minimum and maximum particle
5. In rendering literature, the symbol σ refers to the cross sectionsometimes, while using µ for the coefficient. r min = 0 . , r max = 1 r min = 1 , r max = 10 -5 r min = 1 , r max = 1000 Fig. 8. Top row: Visualizations of the phase functions estimated byour method. Bottom row: The corresponding renderings of a smoothhomogeneous medium with these phase functions, while keeping otherproperties fixed. The red arrows highlight the differences. radii inside V x , respectively. For brevity, we simplify both r min ( V x ) and r max ( V x ) by dropping V x henceforth. For thescattering coefficient and the absorption coefficient, they canbe defined in a similar way by replacing C t ( r ) with C s ( r ) and C a ( r ) , respectively . Generally, these properties exhibitmulti-scale effects with respect to the size of V x .The ensemble phase function is derived as f ( V x , θ ) = (cid:82) r max r min C s ( r ) f p ( θ, r ) (cid:82) x ∈V x N ( r, x )d x d rµ ( V x ) σ s ( V x )= (cid:82) r max r min ( | S ( θ, r ) | + | S ( θ, r ) | ) (cid:82) x ∈V x N ( r, x )d x d r | k | µ ( V x ) σ s ( V x ) . (17)in which σ s ( V x ) serves as the normalization factor for f ( V x , θ ) . Fig. 8 visualizes the phase functions generated bythe above formula. Here, we use uniform sampled radiibetween r min and r max . By fixing other properties, weshow the influence of these phase functions on the finalappearance of a smooth homogeneous medium in Fig. 8.Clearly, this formula is only valid for σ s ( V x ) (cid:54) = 0 . When σ s ( V x ) = 0 , i.e., the volume V x is free of particles, f ( V x , θ ) degenerates into a delta function: σ s ( V x ) f ( V x , θ ) = δ ( θ ) .Similarly, we can derive the transmittance along a lightbeam of length s as T A ( s ) = exp (cid:40) − µ ( A ) (cid:90) r max ( C ) r min ( C ) C t ( r ) (cid:90) x ∈C N ( r, x )d x d r (cid:41) (18)in which C = A × s represents a small cylinder aroundthe light beam and A is its cross section. The derivation isprovided in Appendix D.For particles with a monodisperse distribution, the trans-mittance is simplified to T A ( s ) = exp (cid:26) − C t µ ( A ) (cid:90) x ∈C N ( x )d x (cid:27) (19)
6. These properties can be viewed as the properties at position x when V x is infinitely small, i.e., σ t ( x ) = σ t ( V x ) , σ s ( x ) = σ s ( V x ) and σ a ( x ) = σ a ( V x ) .7. This is similar to the log-normal distribution with a very large σ g . EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 7 in which C t is a constant and N ( x ) is a function of x only.The integral in the above formula simply counts the numberof particles located in the query region A × s . ENDERING S OLUTION
With these bulk optical properties, we are able to derivea multi-scale volumetric rendering equation (VRE) describ-ing radiative transfer in discrete random media. Then, wedevelop a Monte Carlo sampling based solution to solvethe VRE. This solution only requires the position and thesize of each particle, avoiding the explicit tessellation of itsshape. In a preprocessing stage, we generate and store N tot particles with random positions and log-normal distributedradii for a discrete participating medium. During rendering,the stored particles are queried to determine the opticalproperties for each traced ray. A uniform grid is developedfor acceleration. Conventionally, the VRE describing macroscopic light scat-tering in participating media is written as L ( x , ω ) = T ( x , z ) L ( z , ω )+ (cid:90) z T ( x , y ) (cid:90) S σ s ( y ) f ( y , ω (cid:48) , ω ) L i ( y , ω (cid:48) )d ω (cid:48) d y (20)in which L ( x , ω ) is the radiance arriving at x ∈ R along adirection ω ∈ S , f represents the scattering phase functioncharacterizing the probability of radiation incident from ω (cid:48) being scattered into direction ω , z is the distance throughthe medium to the nearest boundary at z = x − z ω and y = x − y ω is a point at distance y ∈ (0 , z ) . The conventionaltransmittance between x and y is computed as T ( x , y ) =exp {− (cid:82) y σ t ( x − s ω )d s } .By substituting the multi-scale properties into the aboveequation, we arrive at a multi-scale version of the VRE: L ( x , ω ) = T A ( x , z ) L ( z , ω )+ (cid:90) z T A ( x , y ) (cid:90) S L i ( y , ω (cid:48) ) Q ( θ )d ω (cid:48) d y (21)with Q ( θ ) = (cid:26) σ s ( V y ) f ( V y , θ ) N ( V y ) (cid:54) = 0 δ ( θ ) N ( V y ) = 0 . (22)Since this multi-scale VRE is a general extension to theconventional one, it naturally supports multiple scattering. In our multi-scale VRE, every optical property depends ona PSD while the PSD is defined on a differential volume.To evaluate the transmittance T A ( x , y ) between any twopositions x and y in the medium, we need a differentialvolume around the ray x → y . This volume is used to queryparticles which contribute to the transmittance T A ( x , y ) .In our implementation, we design it to be a thin cylindercentered around x → y , as illustrated in Fig. 9. We namesuch a cylinder as a query cylinder. In this sense, we vieweach ray as a “fat ray” which gathers small particles along itstrajectory. This is quite different from the implementation of x yx y Query cylinderParticleActive voxel p r p r c Ray
Fig. 9. Illustration of a query cylinder and the active voxels in our raytraversal algorithm. Only the particles contained in the active voxels willbe queried. rendering continuous media in which the optical propertiesare determined globally, without explicitly querying parti-cles in a local area.In theory, the cross section A should be infinitely small.However, a too small one may have numerical issues andcause large variance. Conversely, bias will be introducedin T A ( x , y ) when A is very large. In practice, we select A as follows and keep it unchanged as the ray traversesthe medium. Supposing that z near and z med respectivelydenote the depth of the near plane and the smallest depthof the medium in the view frustum, the size of A is selectedaccording to S A = kS pix z med z near . (23)Here, S pix is the pixel’s size and k can be viewed as thepercentage of the pixel’s footprint at the distance z med .Typically, satisfactory results are obtained when k is in therange [0 . , . The influence of k on the visual effect isdiscussed in the next section. Gathering particles within a query cylinder (central ray: x → y , radius: r c and cross section: A ) requires conductingsphere-cylinder intersection test for every particle in themedium. Given a particle with the position p and radius r p ,it is supposed to be inside the query cylinder if the distancefrom p to the central ray x → y is smaller than r p + r c , asillustrated in Fig. 9 left.Testing all particles of the medium is notoriously time-consuming. To boost the performance, we accelerate theprocess of ray traversing the medium using the 3D digitaldifferential analyzer (3D-DDA) [73], [74]. Specifically, weconstruct a uniform grid for the medium and adopt aray traversal algorithm, similar to that in [73], to find theactive voxels intersected by the query cylinder. Fig. 9 rightillustrates all the active voxels corresponding to the orangequery cylinder. Only those particles inside the active voxelswill be tested against the query cylinder. We determine theactive voxels simply by the ray x → y . This introduces neg-ligible bias, because the radius of the cross section is morethan two orders of magnitude smaller than the side lengthof the voxel. This is significantly different to the thick beamsused in beam radiance estimation [25]. After collecting allthe particles inside the query cylinder, we accumulate theircontributions to the transmittance T A ( x , y ) according to Eq.(18).To construct a uniform grid, its resolution should be care-fully determined. We have observed by experiments that EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 8 high performance is achieved when roughly one particleresides in each voxel after space subdivision. Q ( θ ) To solve the multi-scale VRE, we also have to compute Q ( θ ) at any sample position y . Q ( θ ) describes the angulardistribution of scattering at y . To quickly compute it, asmall query region around y should be defined. We set thisquery region to a small sphere with radius r c . If this queryregion contains N q > particles, we evaluate Q ( θ ) with thefollowing formula: Q ( θ ) = 12 | k | µ ( V y ) N q (cid:88) i =1 ( | S ( θ, r i ) | + | S ( θ, r i ) | ) . (24) Similar to the simplifications used in rendering surfaceglints [52], [53], importance sampling is performed accord-ing to the global optical properties of the medium, assumingit to be continuous. Specifically, we use the global extinctioncoefficient for free-flight sampling and use the tabulatedglobal phase function for angular sampling. These globaloptical properties only need to be determined once in thepreprocessing stage, assuming the entire bounding box ofthe medium to be V x in Eq. (16) and Eq. (17). ESULTS
We have implemented the rendering solution on top of theMitsuba renderer [75], with spectral rendering enabled. Weuse 8 spectral samples in the range of the visible spectrumat equally-spaced locations [46]. After rendering, we convertthe spectral values to the sRGB color space. All synthesizedimages are created on a PC with an Intel 16-core i7-6900KCPU and 16G RAM.To compute the bulk optical properties of a discreteparticipating medium, we need to specify the complex re-fractive index ( η ) for the particles involved and the globalPSD ( ¯ r g and σ g ). The physical unit for the particle radiusis µ m. We also provide an upper bound to the particlesize ( r max = 2000 ) to avoid unreasonably large particleswhich are unusual and are no longer suitable to be treatedas participating media. As mentioned previously, we useGOA in the calculation of the scattering amplitude functionswhen r ≥ and switch to Lorenz-Mie theory otherwise.Except the S TAIRCASE scene, the refractive index is setaccording to the data of ice selected from [17], and the hostmedium is set to be air with η m = 1 . We first compare our method with the traditional pathtracing. Previous methods simulating the grainy appearanceof discrete participating media mostly rely on explicit pathtracing (EPT), with potential approximations to simplifythe computation of high-order scattering [12], [13], [14].However, EPT and other approximations are restricted togeometric optics. This means that only surface reflectionand refraction are properly handled. In principle, the meshof every particle should be explicitly generated and costlyray-object intersections are required. O u r s EP T (a) ¯ r g = 100 and σ g = 6 (b) ¯ r g = 600 and σ g = 3 Fig. 10. Comparisons between EPT and our method in simulating thegrainy appearance of discrete participating media with different PSDs( N tot = 5 · ). (a) r = 10 µ m (b) r = 100 µ m (c) r = 1000 µ m Fig. 11. Angular percentages of energy contributed by the Fraunhoferdiffraction for different sized particles.
Compared with EPT, our method offers at least twobenefits. First, particle scattering is considered which in-cludes the Fraunhofer diffraction and phase differences, etc.This significantly expands the range of particles that canbe handled. Notably, these optical phenomena are quiteimportant in correctly simulating light scattering, especiallyfor very small particles. To verify this, we produce andrender · spherical particles with different PSDs inFig. 10. These particles have random positions and log-normal sampled radii. For EPT, small transparent ballsare instantiated in the scene. Since only surface reflectionand refraction are computed for EPT (the top row), theenergy inherently belonging to the Fraunhofer diffractionis not correctly captured by EPT, resulting in overly sparseand specular volumetric glints in this L AMP scene. Theimportance of the Fraunhofer diffraction is visualized in Fig.11. This figure plots the percentages of energy contributedby the Fraunhofer diffraction at different scattering anglesand for different sized particles. Obviously, the Fraunhoferdiffraction cannot be ignored especially for small particles.Moreover, EPT easily misses many small particles that arehard to be gathered along an ordinary light path, leading toa slow convergence rate. Also, even substantially increasingthe sampling rate or using “fat ray” tracing similar to ours,the appearance is still quite sparse. On the contrary, ourmethod (in the bottom row) preserves these energy from theFraunhofer diffraction and produces smoother appearancethat is closer to real scenarios, thanks to the Airy’s pattern caused by the Fraunhofer diffraction and the efficientrendering solution tailored for sparse media.When only very large particles exist in the medium,
8. The Fraunhofer diffraction will distort light passing through theparticle, making it visually large. The smaller the particle is, the largerthe light’s distribution is.
EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 9
EPT Ours Ours (-FD) FD only
Fig. 12. Comparisons between EPT and our method with or without theFraunhofer diffraction (-FD) in rendering large particles of the samesize (1000 µ m). The impact of the Fraunhofer diffraction is highlighted inthe difference image (FD only). HG OursEPTHG -2 Fig. 13. Comparisons between our directional scattering model andthe best-fit Henyey and Greenstein (HG) model [76]. Due to obviousdiscrepancies between our estimated phase function (red curve in theright diagram) and the HG phase function (blue curve in the rightdiagram), the renderings are slightly differen (see the closeups). our method and EPT tend to produce the similar grainyappearance, as shown in Fig. 12. Here, we render a discretemedium with particles of the same size (1000 µ m).Since the radius is sufficiently large, pure geometric opticsbecomes applicable and serves as a valid approximationto the particle scattering. This is further verified in thethird column of Fig. 11. We see that the contribution of theFrauhofer diffraction concentrates in a very narrow anglein this case, making it hard to be observed. Consequently,if we remove the Fraunhofer diffraction in our model andleave only reflection and refraction (the third colume of Fig.12), we will achieve the grainy appearance similar to that ofour full model. However, there are still subtle differences(highlighted in the difference image) contributed by theFraunhofer diffraction.The second strength of our rendering solution lies in theruntime performance. For rendering these sparse media inthe L AMP scene and the C
ORNELL B OX scene in Figs. 10and 12, our method achieves a . × speedup over EPT inthe rendering time. In our framework, we derive the phase function fromLorenz-Mie theory and GOA. In computer graphics, it ismore common to adopt an empirical model, e.g., the Henyeyand Greenstein (HG) model [76], because of its simplicityand well-defined behavior. However, the HG model is notvery accurate, as pointed out by various literature [69], [77],[78]. In Fig. 13 we compare our estimated phase function with the best-fit HG phase function in rendering the samescene as in Fig. 12. Since the HG phase function cannotfaithfully encode the scattering pattern from these relativelylarge particles ( r = 1000 µ m), the rendering with the best-fit HG phase function is slightly different from ours, ascompared in the insets. Recall that our rendering is closeto that generated using EPT in this specific scene with largeparticles. We also compare the grainy appearance of discrete mediasimulated by our method with the smooth appearance ofcontinuous media. In Fig. 14, we assume that the particlesin the D
RAGON scene possess the same radius r and arerandomly distributed in a cube of the volume V . Underthis configuration, it is easy to derive the extinction coeffi-cient and the scattering coefficient of the continuous mediaas C t ( r ) N tot /V and C s ( r ) N tot /V , respectively. The phasefunction can be computed in a similar way and stored ina table. With these global properties, traditional volumetricpath tracing is applicable to render these continuous media.Generally, each discrete medium and its paired continuousmedium have the same overall brightness. For the discretemedia, when the number of particles N tot is small, we canclearly observe individual particles lit by the lamps. As N tot increases, the appearance tends to become hazy, and therendering result gets closer to the corresponding continuousmedium. When N tot is sufficiently large (e.g., N tot = 10 ),the discrete medium and its continuous counterpart willachieve quite similar appearance. In Fig. 15, a much densermedium ( N tot = 10 ) is rendered by our method, whichagain produces smooth appearance matching that from acontinuous medium. These pair-wise comparisons demon-strate that our rendering solution converges to the tradition-al volumetric rendering of continuous media in the limit ofparticle concentration. We determine the query cylinder’s radius r c according toEq. (24) in which the parameter k plays an important role.We suggest to choose its value in the range [0 . , whichyields reasonable grainy appearance as show in the secondrow of Fig. 16. Values in this range allow us to faithfullycapture almost all grains in the scenes with little bias.Generally, a too large k will produce uncomfortable aliasingas shown in the last two rows of Fig. 16. In these two cases,although the overall brightness is similar to that of k = 0 . ,the volumetric glints are overly blurred. On the other hand,a too small k will miss many particles during query andis therefore inefficient, especially when the medium is verysparse, e.g., N tot = 10 . However, k with a value smallerthan 0.5 is also acceptable sometimes. For instance, the firstimage in the bottom row of Eq. 16 is rendered with k = 0 . which achieves the similar effect with that of k = 0 . forthis relatively dense medium ( N tot = 10 ). We employ a uniform grid to accelerate the ray traversalprocess, considering that particles are uniformly distributed
EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 10 C o n ti nu o u s m e d i u m D i sc r e t e m e d i u m (a) r = 1000 µ m, N tot = 10 (b) r = 400 µ m, N tot = 10 (c) r = 800 µ m, N tot = 10 (d) r = 200 µ m, N tot = 10 (e) r = 400 µ m, N tot = 10 Fig. 14. Visual comparisons against continuous media with varying optical properties. As the particle number N tot increases, the appearance ofthe discrete medium exhibits obvious low-frequency haziness and converges to that of a continuous medium. (a) Discrete medium N tot = 10 (b) Continuous medium Fig. 15. Comparison between our rendering solution and its continuouscounterpart on a very dense medium ( r = 400 µ m and N tot = 10 ). N t o t = N t o t = (a) k = 0 . (b) k = 0 . (c) k = 2 . (d) k = 4 . Fig. 16. Impact of the query cylinder’s cross section on the simulatedgrainy appearance. Here, k is used in Eq. (24) which determines thepercentage of the pixel’s footprint at the nearest boundary of the medi-um. in the scene. As shown in the right diagram of Fig. 17,the grid resolution ( Res ) will influence the performance.Although no analytical analysis can be referred to selectthe best resolution, we empirically observe that a gridresolution yielding roughly one particle per voxel achievesthe optimum solution for most scenes. The B
ENCH scenein the left panel of Fig. 17 contains · particles andachieves the best performance at a resolution of × × .As the resolution increases, the runtime grows steadily dueto the additional cost introduced by the grid. However, aresolution much lower than × × will have a verypoor performance since too many particles reside in each
20 40 60 80 100 120
Resolution T i m e [ m i nu t e ] Fig. 17. Evolution of the rendering time (at 2048 spp) over the gridresolution. This B
ENCH scene contains N tot = 2 · particles ( ¯ r g =800 and σ g = 4 ) and achieves the lowest rendering time at a resolutionof × × which is roughly one particle per voxel. voxel. For other scenes, a similar conclusion can be drawn. The proposed rendering solution can be easily generalizedto support absorbing particles. The absorption cross sectionof each particle is computed by Eq. (13) which varies linearlywith the imaginary part of the complex refractive index, i.e., η i . The extinction cross section is slightly modified accordingto the formulas in Appendix C. The results of changing η i are shown in Fig. 18. This S TAIRCASE scene aims to simulateflying dusts in a dirty room lit by a local area light throughthe window. As expected, an increase of η i will increasinglydim the intensity of scattering. Fig. 19 analyzes the impact of the global PSD on the ap-pearance of discrete participating media. Here, we generate N tot = 4 · particles with different global PSDs inthe H OUSE scene. This scene is designed to simulate theappearance of blowing snows in the sky. The first row showsthe impact of the geometric mean radius ¯ r g . The gener-al trend is that the scattering effects become increasinglyprominent as ¯ r g grows, since the scattering coefficient ispositively correlated with the particle radius. Concerningthe geometric standard deviation σ g , it is responsible for thelevel of graininess as shown in the second row of Fig. 19.A small σ g tends to generate smoother appearance than a EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 11 η i = 0 η i = 2 · − η i = 4 · − η i = 6 · − η i = 8 · − Fig. 18. Impact of the imaginary part of the refractive index ( η i ) on the S TAIRCASE scene. The real part is set to 1.5. Other scene parameters are N tot = 3 · , ¯ r g = 800 and σ g = 4 . ¯ r g = 50 , σ g = 2¯ r g = 100 , σ g = 1¯ r g, = 100 , σ g, = 1¯ r g, = 1000 , σ g, = 1 . r g = 200 , σ g = 2¯ r g = 100 , σ g = 3¯ r g, = 200 , σ g, = 1¯ r g, = 1000 , σ g, = 1 . Fig. 19. Impact of the global PSD. First row: Changing ¯ r g while fixing σ g to 2 and N tot to · . Second row: Changing σ g while fixing ¯ r g to100 and N tot to · . Third row: Changing ¯ r g, in the first mode of thebimodal log-normal distribution while fixing others ( N tot = 2 · ). large one. This is to be expected since a large σ g means localPSDs changing widely, leading to stronger graininess.In fact, the PSD is not limited to the log-normal distri-bution. Other distributions also work. For instance, in thethird row of Fig. 19, we show a cloud of blowing snowsfollowing a bimodal log-normal distribution. We generate . · small particle with ¯ r g, = 100 (or ¯ r g, = 200 )and σ g, = 1 , and also generate large particles with ¯ r g, = 1000 , σ g, = 1 . , leading to · particles in total.Following this complex distribution, we can observe bothhaziness from massive small particles and graininess from ahandful of large particles. TABLE 1Rendering time performance (in minutes) and memory consumption forsix typical scenes shown in this paper. Absorption is not included in thecomputation.
Scene N tot Spp Memory Rendering timeEPT OursL
AMP · ORNELL B OX RAGON ENCH · TAIRCASE · OUSE · T i m e [ m i nu t e ] OursEPT
Fig. 20. Rendering time comparison between our method and EPT onthe C
ORNELL B OX scene shown in Fig. 12. The runtime performance of some test scenes can be foundin Table 1. As we mentioned previously, our renderingsolution achieves a roughly . × speed improvement overEPT for the L AMP scene and the C
ORNELL B OX scene. Asthe particle number N tot increases, the improvement will bemore evident. As shown in Fig. 20, × speed improvementis achieved when N tot increases to in the C ORNELL B OX scene. Since only the position and the radius of each particleare required, the memory consumption of our renderingsolution is affordable even when N tot is very large. The stor-age scales with the number of particles and the resolution ofthe grid. IMITATIONS AND F UTURE W ORK
Although our framework has successfully simulated thegrainy appearance of discrete participating media, it hasseveral limitations deserving further research.
Non-spherical particles.
Our current framework focuses ondiscrete participating media composed of spherical parti-
EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 12 cles. However, non-spherical particles are also very com-mon, e.g., in the context of rendering rainbow [68] or largesnow flakes. Extending our framework to non-spherical par-ticles requires to derive new expressions for the scatteringamplitude functions S and S . For some special particleshapes, determining S and S is rather straightforward andanalytical expressions exist [61], [66], [67]. For more generalshapes, precomputation would be required in practice. Spatial correlations.
As we assume the particles to besparsely distributed, spatial correlations between particlesare not considered and we directly extend the conventionalVRE to support multi-scale graininess. However, for denselypacked particles the conventional VRE becomes question-able due to strongly correlated scattering effects [12], [13],[14]. Recently, new radiative transfer frameworks dedicatedfor correlated media are available in computer graphics [5],[6], [9]. It would be an interesting future work to investigatea more general framework supporting both effects.
Mixture of different particles.
Another possible directionof future work is to explore an efficient strategy to handleparticle mixtures with different compositions. For instance,dusts in real world may be made up of soil particles,textile fibers, human skin cells, etc. A physically-correctparticipating medium should consider such heterogeneousgranular mixtures. For continuous media, this is relativelysimple since the concentrations of different grains are fixed[17]. However, for discrete participating media, the con-centrations are dynamic and change spatially in a similarway as PSDs [14]. Therefore, determining the bulk opticalproperties should take spatially-varying concentrations intoconsideration.
ONCLUSION
We have developed a general and physically-based frame-work for modeling and rendering discrete participatingmedia composed of massive assemblies of independentparticles. Notable characteristics of these media include awide range of PSDs and the appearance of being graininess.To faithfully simulate their appearances, we have derived anovel multi-scale VRE in which a combination of Lorenz-Mie theory and GOA is used to enable high-efficient eval-uation of the important optical properties. A Monte Carlorendering solution is developed to solve the multi-scaleVRE with high accuracy and low computational cost. Wehave extensively evaluated our framework and comparedagainst conventional methods, demonstrating that the pro-posed framework allows us to reproduce a variety of grainyappearances stemming from different discrete participatingmedia and guarantee temporal coherence in animation.Therefore, for the first time in computer graphics, we havegreatly extended the participating media rendering frame-work to handle a much larger range of particle size statistics.We believe the proposed framework is a further step incomputer graphics to manage the details of participatingmedia and expect a more insightful exploration of thisphenomenon in the future research. R EFERENCES [1] E. Cerezo, F. P´erez, X. Pueyo, F. J. Seron, and F. X. Sillion, “Asurvey on participating media rendering techniques,”
The VisualComputer , vol. 21, no. 5, pp. 303–328, Jun 2005. [2] J. Nov´ak, I. Georgiev, J. Hanika, and W. Jarosz, “Monte carlo meth-ods for volumetric light transport simulation,”
Computer GraphicsForum (Proceedings of Eurographics - State of the Art Reports) , vol. 37,no. 2, may 2018.[3] S. Chandrasekhar,
Radiative transfer . Dover, 1960.[4] C. F. Bohren and D. R. Huffman,
Absorption and Scattering of Lightby Small Particles . Wiley, 1983.[5] A. Jarabo, C. Aliaga, and D. Gutierrez, “A radiative transferframework for spatially-correlated materials,”
ACM Trans. Graph. ,vol. 37, no. 4, pp. 83:1–83:13, Jul. 2018.[6] B. Bitterli, S. Ravichandran, T. M ¨uller, M. Wrenninge, J. Nov´ak,S. Marschner, and W. Jarosz, “A radiative transfer framework fornon-exponential media,”
ACM Trans. Graph. , vol. 37, no. 6, pp.225:1–225:17, Dec. 2018.[7] E. d’Eon, “A reciprocal formulation of nonexponential radiativetransfer. 1: Sketch and motivation,”
Journal of Computational andTheoretical Transport , vol. 47, no. 1-3, pp. 84–115, 2018.[8] ——, “A reciprocal formulation of nonexponential radiative trans-fer. 2: Monte carlo estimation and diffusion approximation,” 2018.[9] J. Guo, Y. Chen, B. Hu, L.-Q. Yan, Y. Guo, and Y. Liu, “Fractionalgaussian fields for modeling and rendering of spatially-correlatedmedia,”
ACM Trans. Graph. , vol. 38, no. 4, pp. 45:1–45:13, Jul. 2019.[10] H. C. van de Hulst,
Light scattering by small particles . Dover, 1981.[11] J. Arvo, “Transfer equations in global illumination,” in
GlobalIllumination, SIGGRAPH 93 Course Notes , 1993.[12] J. T. Moon, B. Walter, and S. R. Marschner, “Rendering discreterandom media using precomputed scattering solutions,” in
Pro-ceedings of the 18th Eurographics Conference on Rendering Techniques ,ser. EGSR’07, 2007, pp. 231–242.[13] J. Meng, M. Papas, R. Habel, C. Dachsbacher, S. Marschner,M. Gross, and W. Jarosz, “Multi-scale modeling and rendering ofgranular materials,”
ACM Trans. Graph. , vol. 34, no. 4, pp. 49:1–49:13, Jul. 2015.[14] T. M ¨uller, M. Papas, M. Gross, W. Jarosz, and J. Nov´ak, “Efficientrendering of heterogeneous polydisperse granular media,”
ACMTrans. Graph. , vol. 35, no. 6, pp. 168:1–168:14, Nov. 2016.[15] L. Lorenz, “Lysbevægelser i og uden for en af plane lysbølgerbelyst kugle,”
Det kongelig danske Videnskabernes Selskabs Skrifter ,pp. 2–62, 1890.[16] G. Mie, “Beitr¨age zur optik tr ¨uber medien, speziell kolloidalermetall¨osungen,”
Annalen der Physik , vol. 330, no. 3, pp. 377–445,1908.[17] J. R. Frisvad, N. J. Christensen, and H. W. Jensen, “Computingthe scattering properties of participating media using lorenz-mietheory,” in
ACM SIGGRAPH 2007 Papers , ser. SIGGRAPH ’07, 2007.[18] W. J. Glantschnig and S.-H. Chen, “Light scattering from waterdroplets in the geometrical optics approximation,”
Appl. Opt. ,vol. 20, no. 14, pp. 2499–2509, Jul 1981.[19] A. Ungut, G. Grehan, and G. Gouesbet, “Comparisons betweengeometrical optics and lorenz-mie theory,”
Appl. Opt. , vol. 20,no. 17, pp. 2911–2918, Sep 1981.[20] E. P. Lafortune and Y. D. Willems, “Rendering participating mediawith bidirectional path tracing,” in
EGWR , Vienna, Jun. 1996, pp.91–100.[21] E. Veach, “Robust monte carlo methods for light transport simula-tion,” Ph.D. dissertation, Stanford, CA, USA, 1997.[22] M. Pauly, T. Kollig, and A. Keller, “Metropolis light transport forparticipating media,” in
EGWR , Vienna, 2000, pp. 11–22.[23] L. Szirmay-Kalos, M. Magdics, and M. Sbert, “Multiple scatteringin inhomogeneous participating media using rao-blackwellizationand control variates,”
Computer Graphics Forum , vol. 37, no. 2, pp.63–74, 2018.[24] H. W. Jensen and P. H. Christensen, “Efficient simulation of lighttransport in scenes with participating media using photon maps,”in
SIGGRAPH , Jul. 1998, pp. 311–320.[25] W. Jarosz, M. Zwicker, and H. W. Jensen, “The beam radiance esti-mate for volumetric photon mapping,”
Computer Graphics Forum ,vol. 27, no. 2, pp. 557–566, Apr. 2008.[26] W. Jarosz, D. Nowrouzezahrai, I. Sadeghi, and H. W. Jensen, “Acomprehensive theory of volumetric radiance estimation usingphoton points and beams,”
ACM Trans. Graph. , vol. 30, no. 1, pp.5:1–5:19, Feb. 2011.[27] W. Jarosz, D. Nowrouzezahrai, R. Thomas, P.-P. Sloan, andM. Zwicker, “Progressive photon beams,”
ACM Transactions onGraphics (Proceedings of SIGGRAPH Asia) , vol. 30, no. 6, Dec. 2011.
EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 13 [28] T. Hachisuka, W. Jarosz, I. Georgiev, A. Kaplanyan, andD. Nowrouzezahrai, “State of the art in photon density estima-tion,” in
ACM SIGGRAPH Asia Courses , nov 2013.[29] B. Bitterli and W. Jarosz, “Beyond points and beams: Higher-dimensional photon samples for volumetric light transport,”
SIG-GRAPH , vol. 36, no. 4, pp. 112:1–112:12, Jul. 2017.[30] X. Deng, S. Jiao, B. Bitterli, and W. Jarosz, “Photon surfaces for ro-bust, unbiased volumetric density estimation,”
ACM Transactionson Graphics (Proceedings of SIGGRAPH) , vol. 38, no. 4, Jul. 2019.[31] J. Kˇriv´anek, I. Georgiev, T. Hachisuka, P. V´evoda, M. ˇSik,D. Nowrouzezahrai, and W. Jarosz, “Unifying points, beams,and paths in volumetric light transport simulation,”
ACM Trans.Graph. , vol. 33, no. 4, pp. 103:1–103:13, Jul. 2014.[32] E. Veach and L. J. Guibas, “Optimally combining sampling tech-niques for monte carlo rendering,” in
Proceedings of SIGGRAPH’95 , ser. SIGGRAPH ’95. Association for Computing Machinery,1995, pp. 419–428.[33] B. Miller, I. Georgiev, and W. Jarosz, “A null-scattering pathintegral formulation of light transport,”
ACM Trans. Graph. , vol. 38,no. 4, Jul. 2019.[34] W. Jakob, A. Arbree, J. T. Moon, K. Bala, and S. Marschner,“A radiative transfer framework for rendering materials withanisotropic structure,”
ACM Trans. Graph. , vol. 29, no. 4, pp. 53:1–53:13, Jul. 2010.[35] E. Heitz, J. Dupuy, C. Crassin, and C. Dachsbacher, “The sggxmicroflake distribution,”
ACM Trans. Graph. , vol. 34, no. 4, pp.48:1–48:11, Jul. 2015.[36] J. Dupuy, E. Heitz, and E. d’Eon, “Additional Progress Towardsthe Unification of Microfacet and Microflake Theories,” in
Eu-rographics Symposium on Rendering - Experimental Ideas & Imple-mentations , E. Eisemann and E. Fiume, Eds. The EurographicsAssociation, 2016.[37] S. Zhao, W. Jakob, S. Marschner, and K. Bala, “Building volumetricappearance models of fabric using micro ct imaging,”
ACM Trans.Graph. , vol. 30, no. 4, pp. 44:1–44:10, Jul. 2011.[38] ——, “Structure-aware synthesis for predictive woven fabric ap-pearance,”
ACM Trans. Graph. , vol. 31, no. 4, pp. 75:1–75:10, Jul.2012.[39] P. Khungurn, D. Schroeder, S. Zhao, K. Bala, and S. Marschner,“Matching real fabrics with micro-appearance models,”
ACMTrans. Graph. , vol. 35, no. 1, 2016.[40] S. Zhao, F. Luan, and K. Bala, “Fitting procedural yarn models forrealistic cloth rendering,”
ACM Trans. Graph. , vol. 35, no. 4, Jul.2016.[41] C. Aliaga, C. Castillo, D. Gutierrez, M. A. Otaduy, J. Lopez-Moreno, and A. Jarabo, “An appearance model for textile fibers,”
Computer Graphics Forum , vol. 36, no. 4, pp. 35–45, 2017.[42] S. Zhao, L. Wu, F. Durand, and R. Ramamoorthi, “Downsamplingscattering parameters for rendering anisotropic media,”
ACMTrans. Graph. , vol. 35, no. 6, pp. 166:1–166:11, Nov. 2016.[43] G. Loubet and F. Neyret, “A new microflake model with mi-croscopic self-shadowing for accurate volume downsampling,”
Computer Graphics Forum , vol. 37, no. 2, pp. 111–121, May 2018.[44] L.-Q. Yan, M. Haˇsan, W. Jakob, J. Lawrence, S. Marschner, andR. Ramamoorthi, “Rendering glints on high-resolution normal-mapped specular surfaces,”
ACM Trans. Graph. , vol. 33, no. 4, pp.116:1–116:9, Jul. 2014.[45] L.-Q. Yan, M. Haˇsan, S. Marschner, and R. Ramamoorthi,“Position-normal distributions for efficient rendering of specularmicrostructure,”
ACM Trans. Graph. , vol. 35, no. 4, pp. 56:1–56:9,Jul. 2016.[46] L.-Q. Yan, M. Haˇsan, B. Walter, S. Marschner, and R. Ramamoorthi,“Rendering specular microgeometry with wave optics,”
ACMTrans. Graph. , vol. 37, no. 4, 2018.[47] S. M´erillou, J. Dischler, and D. Ghazanfarpour, “Surface scratches:measuring, modeling and rendering,”
The Visual Computer , vol. 17,no. 1, pp. 30–45, Feb 2001.[48] C. Bosch, X. Pueyo, S. M´erillou, and D. Ghazanfarpour, “Aphysically-based model for rendering realistic scratches,”
Comput-er Graphics Forum , vol. 23, no. 3, pp. 361–370, 2004.[49] B. Raymond, G. Guennebaud, and P. Barla, “Multi-scale renderingof scratched materials using a structured sv-brdf model,”
ACMTrans. Graph. , vol. 35, no. 4, Jul. 2016.[50] S. Werner, Z. Velinov, W. Jakob, and M. B. Hullin, “Scratch iri-descence: Wave-optical rendering of diffractive surface structure,”
ACM Trans. Graph. , vol. 36, no. 6, Nov. 2017. [51] A. Kuznetsov, M. Haˇsan, Z. Xu, L.-Q. Yan, B. Walter, N. K.Kalantari, S. Marschner, and R. Ramamoorthi, “Learning gener-ative models for rendering specular microgeometry,”
ACM Trans.Graph. , vol. 38, no. 6, 2019.[52] W. Jakob, M. Haˇsan, L.-Q. Yan, J. Lawrence, R. Ramamoorthi, andS. Marschner, “Discrete stochastic microfacet models,”
ACM Trans.Graph. , vol. 33, no. 4, pp. 115:1–115:10, Jul. 2014.[53] J. Guo, Y. Chen, Y. Guo, and J. Pan, “A Physically-based Appear-ance Model for Special Effect Pigments,”
Computer Graphics Forum ,2018.[54] B. Wang, L. Wang, and N. Holzschuch, “Fast Global Illuminationwith Discrete Stochastic Microfacets Using a Filterable Model,”
Computer Graphics Forum , 2018.[55] H. Rushmeier, “Input for participating media,” in
ACM SIG-GRAPH 1995 Courses , ser. SIGGRAPH ’95, 1995.[56] P. Callet, “Pertinent data for modelling pigmented materials inrealistic rendering,”
Computer Graphics Forum , vol. 15, no. 2, pp.119–127, 1996.[57] D. Jack`el and B. Walter, “Modeling and rendering of the atmo-sphere using mie-scattering,”
Computer Graphics Forum , vol. 16,no. 4, pp. 201–210, 1997.[58] T. Nishita and Y. Dobashi, “Modeling and rendering of variousnatural phenomena consisting of particles,” in
Computer GraphicsInternational 2001 (CGI’01) , 2001, pp. 149–156.[59] K. Riley, D. S. Ebert, M. Kraus, J. Tessendorf, and C. Hansen,“Efficient rendering of atmospheric phenomena,” in
Proceedingsof the Fifteenth Eurographics Conference on Rendering Techniques , ser.EGSR’04, 2004, pp. 375–386.[60] P. Laven, “Simulation of rainbows, coronas, and glories by use ofmie theory,”
Appl. Opt. , vol. 42, no. 3, pp. 436–444, Jan 2003.[61] E. A. Hovenac, “Calculation of far-field scattering from nonspher-ical particles using a geometrical optics approach,”
Appl. Opt. ,vol. 30, no. 33, pp. 4739–4746, Nov 1991.[62] X. Zhou, S. Li, and K. Stamnes, “Geometrical-optics code forcomputing the optical properties of large dielectric spheres,”
Appl.Opt. , vol. 42, no. 21, pp. 4295–4306, Jul 2003.[63] L. Wu, H. Yang, X. Li, B. Yang, and G. Li, “Scattering by large bub-bles: Comparisons between geometrical-optics theory and debyeseries,”
Journal of Quantitative Spectroscopy and Radiative Transfer ,vol. 108, no. 1, pp. 54 – 64, 2007.[64] H. Yu, J. Shen, and Y. Wei, “Geometrical optics approximation oflight scattering by large air bubbles,”
Particuology , vol. 6, no. 5, pp.340 – 346, 2008.[65] ——, “Geometrical optics approximation for light scattering byabsorbing spherical particles,”
Journal of Quantitative Spectroscopyand Radiative Transfer , vol. 110, no. 13, pp. 1178 – 1189, 2009.[66] H. He, W. Li, X. Zhang, M. Xia, and K. Yang, “Light scatteringby a spheroidal bubble with geometrical optics approximation,”
Journal of Quantitative Spectroscopy and Radiative Transfer , vol. 113,no. 12, pp. 1467 – 1475, 2012.[67] Y. F. Lu, Y. P. Han, J. J. Wang, and Z. W. Cui, “Geometricaloptics approximation for forward light scattering by a large chiralsphere,”
Journal of Quantitative Spectroscopy and Radiative Transfer ,vol. 228, pp. 90 – 96, 2019.[68] I. Sadeghi, A. Munoz, P. Laven, W. Jarosz, F. Seron, D. Gutierrez,and H. W. Jensen, “Physically-based simulation of rainbows,”
ACM Transactions on Graphics , vol. 31, no. 1, pp. 3:1–3:12, 2012.[69] I. Gkioulekas, S. Zhao, K. Bala, T. Zickler, and A. Levin, “Inversevolume rendering with material dictionaries,”
ACM Trans. Graph. ,vol. 32, no. 6, 2013.[70] A. Dal Corso, J. R. Frisvad, T. K. Kjeldsen, and J. A. Bærentzen,“Interactive Appearance Prediction for Cloudy Beverages,” in
Workshop on Material Appearance Modeling , R. Klein and H. Rush-meier, Eds. The Eurographics Association, 2016.[71] V. E. Cachorro and L. L. Salcedo, “New improvements for miescattering calculations,”
Journal of Electromagnetic Waves and Appli-cations , pp. 913–926, 1991.[72] R. A. Dobbins and T. I. Eklund, “Ripple structure of the extinctioncoefficient,”
Applied Optics , vol. 16, pp. 281–282, Feb. 1977.[73] J. Amanatides and A. Woo, “A fast voxel traversal algorithm forray tracing,” in
In Eurographics ’87 , 1987, pp. 3–10.[74] I. Wald, T. Ize, A. Kensler, A. Knoll, and S. G. Parker, “Ray trac-ing animated scenes using coherent grid traversal,”
ACM Trans.Graph.
EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 14 [76] L. G. Henyey and J. L. Greenstein, “Diffuse radiation in theGalaxy.”
The Astrophysical Journal , vol. 93, pp. 70–83, Jan 1941.[77] D. Toublanc, “Henyey–greenstein and mie phase functions inmonte carlo radiative transfer computations,”
Appl. Opt. , vol. 35,no. 18, pp. 3270–3274, Jun 1996.[78] T. Hawkins, P. Einarsson, and P. Debevec, “Acquisition of time-varying participating media,”
ACM Trans. Graph. , vol. 24, no. 3,pp. 812–815, Jul. 2005.[79] C. F. Bohren and D. P. Gilra, “Extinction by a spherical particlein an absorbing medium,”
Journal of Colloid and Interface Science ,vol. 72, no. 2, pp. 215 – 221, 1979.[80] J. Randrianalisoa, D. Baillis, and L. Pilon, “Modeling radiationcharacteristics of semitransparent media containing bubbles orparticles,”
J. Opt. Soc. Am. A , vol. 23, no. 7, pp. 1645–1656, Jul2006.[81] J. Yin and L. Pilon, “Efficiency factors and radiation characteristicsof spherical scatterers in an absorbing medium,”
J. Opt. Soc. Am.A , vol. 23, no. 11, pp. 2784–2796, Nov 2006. A PPENDIX
AA B
RIEF I NTRODUCTION OF L ORENZ -M IE T HEORY
In this section, we briefly describe Lorenz-Mie theory [15],[16] which has already been employed in computer graphics[17], [55], [56], [57], [59]. For light scattering of an elec-tromagnetic wave from a homogeneous spherical particle,exact solutions of the two scattering amplitude functions S and S are given by: S ( θ, r ) = ∞ (cid:88) n =1 n + 1 n ( n + 1) [ a n ( r ) π n (cos θ ) + b n ( r ) τ n (cos θ )] (26) S ( θ, r ) = ∞ (cid:88) n =1 n + 1 n ( n + 1) [ b n ( r ) π n (cos θ ) + a n ( r ) τ n (cos θ )] (27)which express the scattered fields in terms of an infiniteseries of spherical multipole partial waves. Here, a n ( r ) and b n ( r ) are the Lorenz-Mie coefficients of particle size r ; π n (cos θ ) and τ n (cos θ ) are derived from the Legendre func-tions. Please refer to [17] for the details and the expressionsof a n ( r ) , b n ( r ) , π n (cos θ ) , and τ n (cos θ ) .Inserting the expression of S (0 , r ) = S (0 , r ) = S (0 , r ) into Eq. (3), we can obtain a well-defined form of theextinction cross section as [79] C t ( r ) = λ π ∞ (cid:88) n =1 (2 n + 1)Re (cid:26) a n ( r ) + b n ( r ) η (cid:27) . (28)For the scattering cross section, no simple closed-form for-mula is available. It is generally approximated by [80], [81] C s ( r ) = λ e − πr Im { η m } /λ πγ | η m | ∞ (cid:88) n =1 (2 n + 1) (cid:0) | a n ( r ) | + | b n ( r ) | (cid:1) (29)with γ = 2(1 + ( β − e β ) /β and β = 4 πr Im { η m } /λ . Thenotation Re and Im take the real and imaginary part of acomplex number, respectively. The phase function is givenby [10] f p ( θ, r ) = | S ( θ, r ) | + | S ( θ, r ) | π (cid:80) ∞ n =1 (2 n + 1)( | a n ( r ) | + | b n ( r ) | ) . (30) A PPENDIX BD ERIVATION OF C t ( r ) IN E Q . (12) Substituting Eq. (11) into Eq. (3), we have C t ( r ) = 4 π | k | Re S D ,j (0 , r ) + ∞ (cid:88) p =0 S ( p ) j (0 , r ) = 4 π | k | Re α ∞ (cid:88) p =0 α(cid:15) j (0) (cid:115) p/η − e i φ = 4 π | k | Re α (cid:88) p ∈P α(cid:15) j (0) (cid:115) p/η − cos( φ ) = 2 πr + 2 πr | k | (cid:88) p ∈P (cid:15) j (0) | p/η − | cos φ (31)in which P = { , , , · · · } . A PPENDIX
CGOA
FOR P ARTICLES WITH A BSORPTION
For particles with absorption, the refractive index is a com-plex number, which could be written as η p = η r + η i i .Defining the effective refractive index [65]: η (cid:48) = (cid:26)
12 ( η − η + η sin θ i )+ 12 [4 η η + ( η − η − η sin θ i ) ] (cid:27) (32)we have η m sin θ i = η (cid:48) sin θ (cid:48) t (33)where θ (cid:48) t is the effective refractive angle. When particles areabsorbing, the refractive angle θ t should be replaced by θ (cid:48) t .The overall phase shift is changed to φ = (cid:26) φ p + φ f + φ r ,j p = 0 φ p + φ f + φ t ,j p > j = 1 , . (34)The analytical expressions of phase shifts due to reflection φ r ,j and refraction φ t ,j are provided in [65].Moreover, the amplitude functions in Eq. (7) should bemultiplied with the attenuation factor ξ p [65]: ξ p = e − χpα cos θ (cid:48) t /η m (35)considering amplitude attenuation in the absorbing particle.Here, χ is the effective absorption coefficient defined as χ = (cid:26)
12 ( − η + η + η sin θ i )+ 12 [4 η η + ( η − η − η sin θ i ) ] (cid:27) . (36) EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 15 A PPENDIX DD ERIVATION OF THE T RANSMITTANCE IN E Q . (18) Considering a light ray x → y passing through a discreteparticipating medium, the transmittance between x and y = x − s ω is calculated by T ( x , y ) = exp (cid:26) − (cid:90) s σ t ( V x − s (cid:48) ω )d s (cid:48) (cid:27) = exp (cid:40) − (cid:90) s (cid:82) r max r min C t ( r ) (cid:82) x ∈V x − s (cid:48) ω N ( r, x )d x d rµ ( V x − s (cid:48) ω ) d s (cid:48) (cid:41) = exp (cid:40) − (cid:90) s (cid:82) r max r min C t ( r ) (cid:82) x ∈A× d s (cid:48) N ( r, x )d x d rµ ( A ) × d s (cid:48) d s (cid:48) (cid:41) = exp (cid:40) − (cid:82) r max r min C t ( r ) (cid:82) x ∈A× s N ( r, x )d x d rµ ( A ) (cid:41) . (37)Here we set V x − s (cid:48) ω to A × d s (cid:48) . A PPENDIX EM ORE D ISCUSSIONS ON p In GOA, the parameter p is the number of chords that eachray makes inside the particle. The ray is internal reflected p − times before leaving the particle. Since higher-order rays( p > ) have negligible light intensities as compared withother lower-order rays ( p ≤ ), they can be removed in thecomputation of the scattering amplitude functions S and S . To validate this, we plot ( | S | + | S | ) / with increasingvalues of p in Fig. 22 for η = 1 . and in Fig. 23 for p =1 . . As see, when p is small (i.e., p = 1 ), the simulated ( | S | + | S | ) / curves have remarkable differences comparedwith the ground truths generated with a very high order(i.e., p = 100 ). However, the ( | S | + | S | ) / curves with p = 3 and p = 4 are almost identical, and closely matchthe ground truths. This implies that p = 3 is sufficient incomputing S and S with GOA.For evaluating the extinction cross section C t , we canfurther reduce p to 1. This simplification will lower thecomputational cost while introducing negligible error, asverified in Fig. 21. Here, we adopt the Relative MeanSquared Error (RelMSE) between p = 3 and p = 1 : RelMSE { C t } = (cid:16) C p =3t − C p =1t (cid:17) (cid:16) C p =3t (cid:17) (38)to measure the error of C t . As seen, the RelMSE of C t is verylow, especially for r > µ m . A PPENDIX FM ORE C OMPARISONS BETWEEN
GOA
AND L ORENZ -M IE T HEORY
This section provides more visual comparisons betweenGOA and Lorenz-Mie theory. In Fig. 24, we visualize thecurves of log S ( θ ) generated by GOA (red curves) andLorenz-Mie theory (blue curves), respectively. Here, we testtwo different relative refractive indexes: η = 1 . and η = 1 . . The particle radius ranges from . µ m to µ m . -1 R e l M S E -5 (a) η = 1 . -1 R e l M S E -3 (b) η = 1 . Fig. 21. Variation of C t ’s RelMSE as a function of the radius r . Again, close agreements are found when r > µ m withsome differences existing mainly on the backward peaks.When r = 0 . µ m , large errors occur in any direction,indicating that GOA does not work properly in this case.Similar conclusions are reached when comparing GOA andLorenz-Mie theory for the generation of log S ( θ ) curves inFig. 25.Although there are some mismatches between GOA andLorenz-Mie theory in the case of r = 2 µ m , the influenceon the appearance of rendered media is subtle. To see this,we render a smooth cubic medium in Fig. 26 and Fig. 27with different scene configurations. The medium is assumedto comprise monodisperse particles. The extinction coeffi-cient and the phase function are respectively determinedby Lorenz-Mie theory and GOA in a preprocessing stage,according to the particle radius r and the particle number N . However, for r = 0 . µ m we use the same extinctioncoefficient derived from Lorenz-Mie theory in both casessince GOA yields a negative value. This guarantees thefairness of comparison. Nonetheless, quite different appear-ances are achieved by Lorenz-Mie theory and GOA when r = 0 . µ m due to the large discrepancy in S ( θ ) and S ( θ ) . The difference of translucent appearance becomesless noticeable when r goes up to µ m and shrinks furtheras r increases. EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 16 p=100p=1p=2p=3p=4 (a) r = 1 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (b) r = 1 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (c) r = 1 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (d) r = 10 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (e) r = 10 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (f) r = 10 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (g) r = 100 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (h) r = 100 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (i) r = 100 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (j) r = 1000 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (k) r = 1000 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (l) r = 1000 µ m , λ = 0 . µ m Fig. 22. Visual comparisons of ( | S | + | S | ) / with increasing values of p in GOA. Here, we show different combinations of particle size r andwavelength λ , while the relative refractive index of the particle η is set to 1.33. EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 17 p=100p=1p=2p=3p=4 (a) r = 1 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (b) r = 1 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (c) r = 1 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (d) r = 10 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (e) r = 10 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (f) r = 10 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (g) r = 100 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (h) r = 100 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (i) r = 100 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (j) r = 1000 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (k) r = 1000 µ m , λ = 0 . µ m p=100p=1p=2p=3p=4 (l) r = 1000 µ m , λ = 0 . µ m Fig. 23. Visual comparisons of ( | S | + | S | ) / with increasing values of p in GOA. Here, we show different combinations of particle size r andwavelength λ , while the relative refractive index of the particle η is set to 1.40. EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 18
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA (a) r = 0 . µ m Lorenz-MieGOA (b) r = 1 µ m Lorenz-MieGOA (c) r = 2 µ m Lorenz-MieGOA (d) r = 10 µ m Lorenz-MieGOA (e) r = 100 µ m Fig. 24. Visual comparisons of log | S | by Lorenz-Mie calculations (blue curves) with those by GOA (red curves). First row: η = 1 . and λ =0 . µ m . Second row: η = 1 . and λ = 0 . µ m . The particle radius is set to r = 0 . , , , and µ m , respectively. Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA
Lorenz-MieGOA (a) r = 0 . µ m Lorenz-MieGOA (b) r = 1 µ m Lorenz-MieGOA (c) r = 2 µ m Lorenz-MieGOA (d) r = 10 µ m Lorenz-MieGOA (e) r = 100 µ m Fig. 25. Visual comparisons of log | S | by Lorenz-Mie calculations (blue curves) with those by GOA (red curves). First row: η = 1 . and λ =0 . µ m . Second row: η = 1 . and λ = 0 . µ m . The particle radius is set to r = 0 . , , , and µ m , respectively. EEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, IN SUBMISSION. 19 (a) r = 0 . µ m , N = 10 (b) r = 1 µ m , N = 10 (c) r = 2 µ m , N = 2 · (d) r = 10 µ m , N = 10 (e) r = 100 µ m , N = 10 Fig. 26. Rendering a smooth cubic medium with optical quantities derived from Lorenz-Mie theory (top row) and GOA (bottom row), respectively.Here, the relative refractive index η is set to . . (a) r = 0 . µ m , N = 10 (b) r = 1 µ m , N = 10 (c) r = 2 µ m , N = 2 · (d) r = 10 µ m , N = 10 (e) r = 100 µ m , N = 10 Fig. 27. Rendering a smooth cubic medium with optical quantities derived from Lorenz-Mie theory (top row) and GOA (bottom row), respectively.Here, the relative refractive index η is set to .56