Light scattering by objects is generally characterized by a bidirectional scattering distribution function (BSDF) [ Gla95 ].

    (1)

which is the ratio of the scattered radiance dL₀ in the outgoing direction (θ₀, φ₀) to the irradiance L_icos θ_idΩ_i in the direction (θ_i, φ_i) (Figure 1) at wavelength λ. When referring to reflection or transmission, a BSDF becomes a bidirectional reflectance distribution function (BRDF) or a bidirectional transmittance distribution function (BTDF). This paper studies the case of transmission.

Figure 1. Light scattering at a surface (transmission case).

In computer graphics applications, materials may be classified into three major types: opaque, transparent and translucent. An opaque object only involves reflection, a transparent object involves both reflection and transmission, and a translucent object has volumetric scattering in addition to reflection and transmission at the object surface. Thus, a transmission model is needed for not only transparent but also translucent objects. Such objects include glass wares, plastics, ices, biological tissues, marbles, waxes, and so on.

There has been extensive research on modelling BRDFs in computer graphics, but studies on BTDFs are limited. Different from at an opaque surface, a scattering process at a surface of some transparent or translucent material is generally a combination of reflection and transmission events, and the number of the events may be one (single scattering), two or more (multiple scattering). Solving the case of single scattering is a basis for solving the case of multiple scattering.

This paper presents a realistic transmission model of rough surfaces. We follow the same assumptions of rough surfaces as presented in the surface reflection model [ Sun07 ]. The key point of these assumptions is that the surface is sufficiently smooth locally (tangent plane approximation) and statistical techniques can be applied to calculate light transmission through a local illumination area (statistical approach efficiency). We have obtained an analytical expression for single scattering. This model has been compared to our Monte Carlo simulations as well as to the previous simulations, and good agreements have been achieved.

Existing BRDF models commonly consist of the diffuse and specular terms. The diffuse component is typically Lambertian, but the specular term differs in various models. A simple approach describes the specular component with an empirical function, such as the models of Phong [ Pho75 ], Ward [ War92 ], and Lafortune et al. [ LFTG97 ].

Deriving accurate models needs physically based approaches. One approach uses the Kirchhoff theory with the tangent plane approximation of the surface [ BS63, HTSG91 ]. Another approach is based on the microfacet assumption of Torrance and Sparrow [ TS67 ]. In this approach, the specular term is expressed as a product of the Fresnel coefficient, masking and shadowing factor, and surface orientation probability [ Bli77, CT82 ]. Ashikhmin et al. [ APS00 ] developed an analytic model to remove the limitation of V-shaped grooves needed for the traditional microfacet model. Recently, Sun [ Sun07 ] reviewed all the previous illumination models, refined the assumptions of rough surface, and derived an illumination model of rough surface by using statistical techniques.

To our best knowledge, two transmission models exist in computer graphics. The first was proposed by He [ He93 ] based on the Kirchhoff theory, and the second by Stam [ Sta01 ] as an extension from Cook-Torrance′s reflection model [ CT82 ]. In practice, the rendering of light transmission is rather simple, typically based on a formula that extends Phong′s reflection model to the case of transmission.

Beyond computer graphics, some research has been conducted to numerically simulate transmission. One example is the work of Nieto-Vesperinas et al. [ NVSGSD90 ] where light transmission at rough surfaces was computed using a Monte Carlo method.

Since transmission at a surface is a part of the problem of object translucency, we briefly review some research on translucency models. Hanrahan and Krueger [ HK93 ] developed a pioneering model of subsurface scattering using the linear transport theory. Jensen and Christensen [ JC98 ] studied light transport in participating media using Monte Carlo bi-directional ray tracing and volumetric photon mapping. Dorsey et al. [ DEJ99 ] simulated subsurface scattering of weathered stones using Monte Carlo ray tracing. Pharr and Hanrahan [ PH00 ] developed a Monte Carlo approach to solve generic scattering equations. Stam [ Sta01 ] used the radiative transfer equation to model subsurface scattering of human skins. Koenderink and van Doorn [ KvD01 ] studied subsurface scattering with a diffusion approximation of light transport theory. Jensen et al. [ JMLH01 ] proposed an analytic model of BSSRDF, and later Jensen et al. [ JB02 ] developed a two-pass technique to efficiently render translucent objects. Recently, Wang et al. [ WTL05 ] presented a technique based on pre-computed light transport to render translucent objects, and Mertens et al. [ MKB05 ] proposed an efficient algorithm to render the local effect of subsurface scattering. These studies focused on the subsurface or volumetric scattering, and light transmission at the surface was not considered.

Light transmission processes at a rough surface can be classified into single and multiple scattering. In single scattering (ray 1 in Figure 2), a light ray is scattered one time (this is in fact a refraction at the local area). In multiple scattering (ray 2 or 3 in Figure 2), there are multiple times of reflection and transmission. The total BTDF may be expressed as

ρ_total = ρ_single + ρ_multiple     (2)

where ρ_single and ρ_multiple are the contributions from single and multiple scattering, respectively.

Figure 2. Light transmission processes of single scattering (ray 1) and multiple scattering (rays 2 and 3).

The assumptions and conditions of our considered surface are similar to the surface reflection model [ Sun07 ] where the focus was on reflection, but now the focus is on transmission. For convenience, the assumptions are listed below:

For most natural surfaces, the surface height distribution can be described well with a Gaussian function [ Tho99 ]. The probability density function of surface height is

    (3)

where ζ is the surface height, and σ is the standard deviation or RMS.

To describe surface roughness, we need to consider the surface height correlation. A two-point correlation function is generally defined as

C(r) = < h(r₀)h(r₀ + r) >/σ¹  (4)

which involves the average of the product of heights at points r₀ and r₀ + r on the z = 0 plane. Since the surface is isotropic, we can write C(r) = C(r). We adopt the following correlation function [ Ogi91 ]

C(r) = exp(-r²/τ²)    (5)

where τ is the correlation length. Now we define the surface smoothness as

s = τ/σ    (6)

The smaller is s, the rougher the surface; vice versa.

Given surface profile ζ = h(x,y) the orientation of a micro-area δA is described by the partial derivatives (ζ'_x, ζ'_y)

    (7)

The probability for the orientation of a micro-area δA in dζ'_xdζ'_y is given as (in [ Sun07 ])

    (8)

where dΩ_n is differential solid angle of n(θ_n,φ_n)

dΩ_n = sinθ_ndθ_ndφ_n .    (9)

Given a micro-area δA (Figure 3), the incident radiance L_i(θ_i, φ_i) and the transmitted radiance L₀(θ₀, φ₀) are related as,

L₀(θ₀, φ₀)cosβdΩ₀ = F_t(α,λ)L_i(θ_i,φ_i)δ(n,e_i,e₀)cosαdΩ_i     (10)

Figure 3. Ray transmission at a micro-area.

where β is the transmission angle for the incident angle α, dΩ₀ is the solid angle in the transmission direction, dΩ_i is the solid angle in the incident direction, and F_t(α, λ) is the Fresnel coefficient of transmission averaged over polarizations at the wavelength λ. δ(n,e_i,e₀) is a Dirac delta function. That is, when n,e_i , ande₀ are coplanar and sinα = n sinβ (n is the relative index of refraction), δ(n,e_i,e₀) = 1; otherwise δ(n,e_i,e₀) = 0. The radiant flux δΦ(θ₀,φ₀) through a micro-area δA is given as

δΦ(θ₀,φ₀) = L₀(θ₀,φ₀)cosβδAdΩ₀ = F_t(α, λ)(θ_i,φ_i)δ(n,e_i,e₀)cosαδAdΩ_i     (11)

Since a local illumination area ΔA over which the BTDF is defined (Figure 1) contains many micro-areas δA, the total radiant flux overΔA contains contributions from all possible micro-areas,

ΔΦ(θ₀,φ_,0) = ∑_ΔAV(δA)δΦ(θ₀,φ₀) = L_i(θ_i,φ_i)cosαdΩ_i∑_ΔAV(δA)F_t(α, λ)δ(n,e_i,e₀)δA    (12)

where V(δA) is a visibility function describing the probability of a micro-area δA that is visible in both directions e_i(θ_i,φ_i) and e₀(θ₀,φ₀). The radiance to the transmission direction e₀(θ₀,φ₀) is given as

    (13)

Because θ₀ is measured from the positive z-axis (see Figure 1) and its value is within [π / 2, π], we take its absolute value in the above equation. Substituting Eq. 13 into Eq. 1, we obtain

    (14)

Here the summation over the local illumination area ΔA has been decomposed into the summation over all micro-areas at a fixed height ζ and the summation over all micro-areas at different heights. Since the visibility function V(δA) at a fixed height remains the same for given incident and outgoing directions, it has been put outside the inner summation for a fixed height. Considering that the projected area of δA on the z = 0 plane is

(δA)^⊥ = δAcosθ_n    (15)

where θ_n is the polar angle of the normal n(θ_n, φ_n) of δA (see Figure 3), the portion of the total projected areas ∑_{fixed ζ}(δA)^⊥)δ(n,e_i,e₀) in ΔA is the probability of a surface point with height in differential interval [ζ,ζ+dζ] and with orientation in intervals [ζ'_x,ζ'_x+dζ'_x] and [ζ'_y,ζ'_y+dζ'_y]. Thus

    (16)

From Assumption 5, the combined probability density function can be decomposed as

p(ζ,ζ'_x,ζ'_y) = p(ζ)p(ζ'_x,ζ'_y).    (17)

Applying Eqs. 16, 17 into Eq. 14, we obtain

    (18)

This equation may be further expressed as

    (19)

where the function χ(α,β) describes dΩ_n/dΩ₀ (see Appendix A), and

    (20)

is the averaged bistatic visibility function. A bistatic visibility function simultaneously involves the incident direction e_i and the outgoing direction e₀ . For light transmission, since e_i points into the original medium and e₀ into the new medium, the correlation between the two directions can be ignored, as stated in Assumption 6. Therefore

V(ζ, e_i,e₀) ≈ V(zeta, θ_i)V(ζ, θ₀),    (21)

where V(ζ, θ) is an individual visibility function that describes the probability of being visible for a ray starting at height ζ and with angle θ (Figure 4), and accordingly,

    (22)

We further approximate Eq. ( 22) as

    (23)

where V(0,θ_i) and V(0,θ₀) are the individual visibility functions for the incident and outgoing directions when the ray starts from ζ = 0. From the previous study [ Sun07 ],

    (24)

where k₀ = 0.7. Thus, we finally obtain

    (25)

where χ(α, β) is given in the Appendix.

Figure 4. A ray starts at height ζ and with polar angle θ.

Figure 5.  BTDF for different n and s. Parameters are θ_i = 30°, n = 1/1.4 for the first row, and n = 1.4 for the second row. From the left to right, the values of s are 6, 3, and 1, respectively.

Figure 5 shows ρ_single for different values of relative index of refraction (IOR) n and smoothness s. The solid straight lines (in green) in the upper hemisphere indicate the incident direction, and the solid straight lines (in blue) in the lower hemisphere indicate the transmission direction. ρ_single has a sharp lobe and shows the off-specular effect. When n < 1 (the first row), as the outgoing direction changes from θ = -90° to θ = 180°, ρ_single increases gradually and reaches a maximum, then decreases rapidly. Also, the direction that ρ_single has the maximum shifts toward θ = 180° with the decrease of s. In contrast, for n > 1 (the second row), when the outgoing direction changes from θ = -90° to θ = 180°, ρ_single increases rapidly and reaches the maximum, then decreases gradually. Moreover, the direction for the maximum ρ_single shifts towardθ = -90° with the decrease of s.

The plots in Figure 5 can be explained as below. First, when the surface is smooth, most micro-areas distribute around θ_n = 0° and they contribute to ρ_single with α → θ_i . Second, Fresnel's transmission coefficient has the maximum for incident angle α = 0°, and decreases with the increase of α. Therefore, those micro-areas with orientation around the incident direction have large Fresnel's transmission coefficients. These two factors compete with each other. And also, for n < 1, the refraction angle β is larger than the incident angle α, and vice versa. These result in the plot shapes in Figure 5. With the decrease of s, the maximum distribution of orientations of micro-areas tends to shift from θ_n = 0° toward θ_n = 90°, which results in a shift of the direction for the maximum ρ_single .

Figure 6.  BTDF for different n and s. Here θ_i = 0°, and other parameters and notations are the same as Figure 5

For the normal incidence, ρ_single for different values of n and s is shown in Figure 6. For n < 1, the sharp lobe becomes wider with the decrease of s, same as Figure 5. However, for the cases n > 1 in both Figures 5 and 6, although the sharp lobes for s = 3 are all wider than those for s = 6, the sharp lobes for s = 1 have different shapes. Consider the rotational geometry, ρ_single for s = 1 in Figure 6 is actually a lobe with an indented peak. We can understand this by the micro-area model. For rough surfaces (s = 1), most micro-areas distribute with orientations θ_n ≫ 0°, and therefore the transmitted light by single scattering tends to travel along the direction with 90° ≪ θ₀ ≪ 180°. This results in the indentation of the lobe in Figure 6. However, the probability of ray blocking is higher for the rays propagating along this direction. This results in the sharper shape for s = 1 in Figure 5.

In our Monte Carlo simulation, given a Gaussian rough surface with its mean equal to zero and standard deviation σ, totally N light rays are shot from the incident direction e_i , each ray carrying a weight W_l(l=1,2,...,N) that represents its radiance flux intensity. Once a shot ray hits the surface profile, it typically splits into a reflected and a transmitted ray. When the total internal reflection occurs, only a reflected ray is generated.

The surface height at which a shot ray intercepts with the profile is determined by the probability density function of surface height and a generated random number (all the generated random numbers in this paper are uniformly distributed between 0 and 1). The normal direction of this intersection point is obtained by the orientation probability density function with two generated random numbers.

We set the incident flux density to 1. Then the weight of the lth shot ray is given as

W_l = V(ζ₁, θ_i)cosα/cosθ_n     (26)

where ζ₁ is the surface height that the shot ray first intercepts with, α is the incident angle in the local area, and cosθ_n is involved because Eq. 16 just describes the probability distribution of δA^⊥ at a fixed height.

When a ray with weight W hits the surface profile, it splits into a reflected and transmitted ray, and the weight of the reflected ray is F_r(α, λ)W and that of the transmitted ray is F_t(α, λ)W. Therefore, after each ray splitting, the generated rays will decrease in intensities. Once the weight of a newly generated ray is lower than the threshold, the tracking process terminates. Otherwise, it will be tracked continuously; whether it is blocked or not depends on its propagation direction, visibility function, and a generated random number.

The radiance to the transmission direction e₀(θ₀,φ₀) is obtained as

    (27)

where ΔΩ₀ is the solid angle along e₀(θ₀,φ₀), and ∑_{l ∈ ΔΩ0} W_l calculates the sum of the weights of those rays transmitted into ΔΩ₀ . Consider that the incident irradiance is cosθ_i (since incident flux density is set to 1), the BTDF can be calculated by

    (28)

In discussion below we may replace ρ with ρ ∣ cosθ₀ ∣ based on two considerations. First, the previous simulations by Nieto-Vesperinas et al. [ NVSGSD90 ] calculated the transmitted light intensity, which is proportional to ρ ∣ cosθ₀ ∣. For convenience of comparing the results, we need to use ρ ∣ cosθ₀ ∣ instead of ρ. Second, Eq. ( 28) contains 1/ ∣ cosθ₀ ∣ and ∑_{l ∈ ΔΩ0} W_l . When θ₀ → π, ∣ cosθ₀ ∣ → 0. However, we cannot take ΔΩ₀ → 0 for the calculation of Eq. (28). Therefore, ρ might diverge at θ₀ → 90°.

Figure 7 compares our analytical model and simulations. Nieto-Vesperinas et al. [ NVSGSD90 ] considered perpendicular and parallel polarizations separately. For comparison, we calculate the average of the two polarizations. In our analytical model and simulation, light intensity can be calculated by ρ ∣ cosθ₀ ∣ ΔA. Since we do not know the value of ΔA used for the simulation of Nieto-Vesperinas et al., we find it by matching our analytical model with their results for θ_i = 0°. In Figure 7, the comparison shows a very good match.

Figure 7.  Comparison between our analytical model and simulations. The curves with x marks are from the analytical model, the dot curves from the simulations of Nieto-Vesperinas et al. [ NVSGSD90 ], and the solid curves from our simulations. Here, n = 1.411, s = 2.522, (a) θ _i = 0°, (b) θ _i = 20°, (c) θ _i = 40°, and (d) θ _i = 60°.

Figure 8.  Comparison between simulation and analytical model. The solid curves are from our simulation, and the curves with x marks from analytical model. Here, θ _i = 30°, n = 1/1.4, (a) s = 6, (b) s = 3, (c) s = 1, and (d) s = 0.5.

Figure 9.  Comparison between simulation and analytical model. The solid curves are from our simulation, and the curves with x marks from analytical model. Here, θ _i = 30°, n = 1.4, (a) s = 6, (b) s = 3, (c) s = 1, and (d) s = 0.5.

Figures 8 and 9 compare our simulation and the analytical model for different values of n and s. For smooth and moderately smooth surfaces (s is 3 or 6), the analytical model agrees well with the simulation. With small s, the difference between the model and simulation increases. This is because our analytical model only considers single scattering. For smooth surfaces, light transmission is dominated by single scattering. Overall, the model has a good match with the simulation. For rough surfaces (s is small), multiple scattering plays an important role and should be considered.

We have implemented this analytical model for image rendering. In each image, we have linearly scaled the color components so that the maximum value of all color components is 255.

Figure 10 shows the rendered transparent object with semi-infinite depth for different n and s. The face of this transparent object is a rough, plate surface. The distance from the eye to the center of the face is 0.8 times width of the face. A far point light source is set inside the object medium. Therefore, the incident light is almost perpendicular to the medium surface. Figure 11 shows the rendered images of the same transparent object with the same parameters and notations as Figure 11 , except the point light source is closer to the face.

Figures 12 and 13 show rendered sculptures for different s. In each image, a far point light source is set in the back. The distance from the eye to the sculpture is about 0.75 times the size of the sculpture. The index of refraction of the sculptures is 1.4. In Figure 13 , the foreheads in the images for s = 4 and s = 5 do not look smooth, and this is due to the fact that the mesh of this sculpture is not dense.

Figure 10.  Rendered transparent objects with semi-infinite depths for different n and s. Here the parameters and notations are the same as Figure 5.

Figure 11.  Rendered transparent objects with semi-infinite depths for different n and s. A near point light source is set inside the object medium (the distance from the point light source to the face is 0.8 times width of the face). Other parameters and notations are the same as Figure 10.

Figure 12.  Rendered images of Atenea for different s. From the left to the right, the values of s are 2, 3, 4, and 5, respectively.

Figure 13.  Rendered images of Angel for different s. The values of s are taken as same as Figure 12.

This paper presents a physically based transmission model of rough surface. We have obtained an analytical expression for single scattering. The analytical expression has been compared to our Monte Carlo simulations as well as to the previous simulations, and good agreements have been achieved. This model has been used for image rendering.

In future work, the model can be applied to render realistic transmission effects. And also, it could be taken into consideration to study object translucency. To further verify this model, we may generate 2D surfaces for given surface height standard deviation and correlation length, and compute the average of transmission through the surfaces. The current analytical model has not considered multiple scattering, and both the analytical model and simulation have not considered polarization effects. We will consider them in our further work.

Here we derive the relationship between the differential solid angles dΩ₀ and dΩ_n . Given an unit sphere (Figure 14), the area ABCD corresponds to dΩ_n and the area A'B'C'D' to dΩ₀ . The points A, B, A', and B' are coplanar, and similarly the points C, D, C' and D'. The planes ABA'B' and DCD'C' intersects at the line POQ, and the angle between them is dγ. Therefore, the length of the curve segment AD is

| AD | = | AE | ⋅ dγ = sinα ⋅ dγ.    (A1)

Since | AB | = dα, so we obtain

dΩ = | AD | ⋅ | AB | = sinα ⋅ dγdα    (A2)

Figure 14. Relationship between dΩ_n and dΩ₀ (n > 1).

The Snell's law gives the following relations:

sinα = n⋅sinβ    (A3)

and

sin α' = sin(α + dα) = n⋅sinβ' = n⋅sin(β + dβ)    (A4)

where n is the relative index of refraction n > 1 and dβ is defined as

dβ = β' - β    (A5)

From Eq. A4 we can obtain

sinα ⋅ cos(dα) + cosα ⋅ sin(dα) = n ⋅ [sinβ ⋅ cos(dβ) + cosβ ⋅ sin(dβ)]    (A6)

We make the following approximations:

cos(dα) ≈ 1, sin(dα) ≈ dα, cos(dβ) ≈ 1, sin(dβ) ≈ dβ.    (A7)

Substituting Eqs. (A3) and (A7) into Eq. (A6), we obtain

    (A8)

From Figure 14(b) we obtain

β' = β - ∠A'OB' + dα.    (A9)

Therefore, we obtain

    (A10)

The length of the curve segment A'B' is

    (A11)

The length of the curve segment A'D' is

| A'D' | = | A'E' | ⋅ dγ = sin(∠A'OQ)dγ = sin(α - β)dγ    (A12)

Therefore, we obtain

    (A13)

Finally, we obtain

    (A14)

Although Eq.(A14) is derived for n > 1, it is easy to prove that this expression also holds for n < 1.