=Paper=
{{Paper
|id=Vol-2744/paper16
|storemode=property
|title=Methods Calculating the Slab Radiance Factor
|pdfUrl=https://ceur-ws.org/Vol-2744/paper16.pdf
|volume=Vol-2744
|authors=Vladimir Budak,Dmitry Efremenko
}}
==Methods Calculating the Slab Radiance Factor==
https://ceur-ws.org/Vol-2744/paper16.pdf
Methods Calculating the Slab Radiance Factor*
Vladimir Budak1[0000-0003-4750-0160], Dmitry Efremenko2[0000-0002-7449-5072]
1National Research University “Moscow Power Engineering institute”, Moscow, 111250 Russia
budakvp@gmail.com
2Remote Sensing Technology institute (IMF), German Aerospace Center (DLR), Oberpfaffenhofen,
82234 Wessling, Germany
dmitry.efremenko@dlr.de
Abstract. One of the most critical problems of realistic visualization of the real-world objects
is physically adequate modeling of their reflection of light. Reflection of light by objects
occurs both from the surface and the bulk of matter (scattering). Accounting for the light
reflection from the surface of objects was solved almost a century ago based on its
representation as a Fresnel randomly rough surface. Scattering by a bulk of matter is the
subject of radiation transfer theory, which has only recently received its known completion in
the form of discrete transfer theory. Strict analytical methods for solving the radiation
transport equation (RTE) are often not highly effective for calculating the radiance factor. For
a long time, in the absence of effective numerical methods for solving problems and the
unavailability of high-speed computers for practical calculations, approximate methods for
solving RTE were developed. However, their accuracy and applicability limits were poorly
defined. The discrete transfer theory allowed us to evaluate the existing approximate methods
for solving the UPI, their accuracy, and the efficiency of application for calculating the
radiance factor. It is shown that the most effective method is the method of synthetic
iterations. The method is based on the selection of the solution anisotropic part based on a
small-angle approximation of the RTE solution. The solution regular part can be calculated
by any approximation. Then a simple iteration from the complete solution is performed to
refine the angular distribution of the radiance factor.
Keywords: Radiance Factor, Reflectance, Radiative Transfer, Synthetic Iteration,
Discrete Ordinates.
1 Radiance factor
One of the most pressing problems of computer graphics (CG) today is the mathematical
model of reflection by objects of the real world [1, 2]. The concept of light reflection was
introduced to optics by the first physicist on earth, Abu Ali al-Haysam [3], in the XI century.
Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons License
Attribution 4.0 International (CC BY 4.0).
2 V. Budak et al
It is difficult to list all the names in a short article that the theory of reflection owes its
development over 1000 years. Common reflection characteristics were introduced by C.
Fabry [4]. A complete reflection characteristic is the surface radiance factor, which,
according to the International lighting vocabulary [5], is determined as the ratio of the
radiance of the surface element in the given direction to that of the perfect reflecting or
transmitting diffuser identically irradiated and viewed.
Since the lighting conditions are not predefined, the radiance factors for diffuse and
directional lighting should be distinguished. Let us consider the radiance L(r, ˆl ) at the
point defined by the radius-vector r and propagating along the direction ˆl . Then, the
radiance factor (r; ˆl, ˆl ) in the viewing direction l̂ due to the directed illumination of the
studied surface area at the point r in the direction ˆl is defined as follows:
𝐿(𝒓,𝒍̂′ )
𝜌(𝒓; 𝒍̂, 𝒍̂′ ) = 𝜋 ̂ ,𝒍̂)
, (1)
𝐸(𝒓)(𝑵
where L(r, ˆl ) is a visible surface radiance, E (r ) is normal irradiance from incident
radiation, N̂ is the surface normal at the point r.
It is straight forward to see that according to (1) the radiance factor is the same as the
definition of bidirectional reflectance distribution function (BRDF). However, according to
[4], the term radiance factor is applicable for both reflected and transmitted radiation.
Accordingly, if a thin horizontal layer is irradiated from above L (r, ˆl ) and below
L (r, ˆl ), the radiance L(r, ˆl ) of the radiation coming out of the layer can be represented
as:
𝐿(𝒓, 𝒍̂)|𝜇<0 1 𝜌(𝒓; 𝒍̂′ , 𝒍̂)𝛩(𝜇′ > 0) 𝜏(𝒓; 𝒍̂′ , 𝒍̂)𝛩(𝜇′ < 0) 𝐿↓ (𝒓, 𝒍̂′ ) ′ ̂ ′
[ ] = ∮[ ][ ] |𝜇 |𝑑𝒍 , (2)
𝐿(𝒓, 𝒍̂)|𝜇>0 𝜋 𝜏(𝒓; 𝒍̂′ , 𝒍̂)𝛩(𝜇′ > 0) 𝜌(𝒓; 𝒍̂′ , 𝒍̂)𝛩(𝜇′ < 0) 𝐿↑ (𝒓, 𝒍̂′ )
where (r; ˆl , ˆl ) is the radiance factor by transmission, = (N ˆ , ˆl ) , = (N
ˆ , ˆl ) . The
function is equal to 1 if the argument value is true, and 0 if it is false.
Thus, in the general case, the slab is characterized by a matrix of radiance factor. For
radiance factors, the principle of reciprocity is observed: transmission and reflection from
top to bottom and from bottom to top are equal.
In [6], a model of slab reflection was proposed, which determined all further efforts in
this direction. The model [6] is based on the representation of the reflection process as the
sum of two processes: a reflection from the Fresnel randomly rough surface (RRS)of the
interface between media and light scattering by the thickness of substances. The
development of the RRS reflection model has been very productive. In the article [7], the
primary relations of reflection from the RRS were formulated, and subsequent works were
mainly its refinements.
Methods Calculating the Slab Radiance Factor 3
The problem of reflection (scattering) of light by the entire volume of the layer is
reduced to a boundary value problem for the radiation transfer equation (RTE), which
turned out to be a mathematically more complex problem, practical algorithms for solving
which in fact only by the end of the XX century. In [8, 9], a complete, analytically rigorous
model of radiation reflected by a layer was proposed. With all the advantages of this model,
it requires setting environment such characteristics that are not clear to the average user of
the software.
However, before precise methods for solving RTE were developed, approximate
methods were developed that were limited in their scope. They were commonly referred to
as engineering methods since they were often used in the design of various devices. These
methods were often developed separately, based on some narrow assumptions about the
properties of the medium and the light field inside. From the present time, when a strict
analytical solution of the RTE is available, it is of considerable interest to assess the
accuracy and applicability limits of these engineering approximations. It will allow us to
determine their suitability for use in CG algorithms.
2 Discrete radiative transfer theory
We start our analysis with the boundary value problem formulation for the radiance field
L(, ˆl ) in the case of a homogeneous slab of the optical thickness 0 . Here is the vertical
coordinate measured in dimensionless optical thickness and is zero at the layer top, while
ˆl = , is the direction given in spherical coordinates by the cosine of polar angle and
the azimuthal angle . The slab is illuminated by a mono-directional plane source at the
top of the layer with ˆl = , being the direction of incidence, the cosine of the
0 0 0 0
polar incident angle, while 0 = 0 . The boundary value problem for L(, ˆl ) reads as
follows [10]:
𝜕𝐿(𝜏,𝒍̂) 𝛬
𝜇 + 𝐿(𝜏, 𝒍̂) = ∮ 𝐿(𝜏, 𝒍̂′ )𝑥(𝒍̂, 𝒍̂′ )𝑑𝒍̂′ ,
{ 𝜕𝜏 4𝜋
(3)
𝐿(𝜏, 𝒍̂)|𝜏=0,𝜇>0 = 𝛿(𝒍̂ − 𝒍̂0 ), 𝐿(𝜏, 𝒍̂)|𝜏=𝜏 ,𝜇<0 = 0,
0
where x is the single scattering phase function.
Direct numerical solution of Eq. (3) is impossible due to singularities in L(, ˆl ) (i.e. the
delta-function) coming from the upper boundary conditions. To overcome this problem, the
total radiation is expressed as a sum of the direct radiation and the diffuse radiation [11].
Note that the first (singular) part can be evaluated analytically, while the diffuse part can be
found numerically. In practice, this approach is not efficient when the single scattering
phase function is forward peaked (e.g. due to large particles compared to the wavelength of
4 V. Budak et al
the light [12]), for which several hundreds of Legendre expansion coefficients have to be
taken into account. For the sake of numerical efficiency, in addition to the singular part, the
so called small-angle part of the radiance is treated analytically, while the remaining part of
the radiance field is smooth. This framework is based on the following representation of the
radiance field [13]:
𝐿(𝜏, 𝒍̂) = 𝐿𝑎 (𝜏, 𝒍̂) + 𝐿̃(𝜏, 𝒍̂), (4)
where La (, ˆl ) is the anisotropic part which incorporates all singularities in the angular
domain, while L%(, ˆl ) is the regular part of the solution. As shown in [13], substituting Eq.
(4) into Eq. (3) gives a new boundary problem with an additional source term. The
numerical solution based on the discrete ordinate method involves the following steps:
• Fourier cosine expansion of the radiance and phase functions;
• discretization of the radiance field in the angular domain by introducing a set of
Gaussian points i and weights wi , i = 1,..., N , where N is the number of discrete
ordinates per hemisphere, the upper index ‘+’ refers to the downwelling radiance, while
‘-’ corresponds to the upwelling radiance.
• replacing the scattering integral with a finite sum [14] in the discrete ordinate space;
• transformation of the radiative transfer equation into the differential matrix equation
with respect to column vectors consisting of the radiance values at Gaussian nodes.
After performing these steps, the original problem is transformed into the following
equation:
𝑑𝐶⃗ (𝜏)
⃡𝐶⃗(𝜏) + 𝑀
= −𝐵 ⃡⃗−1 Δ⃗(𝜏), (5)
𝑑𝜏
𝑇
where 𝐶⃗(𝜏) ≡ [𝐶⃗− (𝜏), 𝐶⃗+ (𝜏)] , 𝐶⃗± (𝜏) ≡ {𝐿𝑚 (𝜏, 𝜇𝑖± )}𝑖=1,𝑁 , B is the layer matrix, Δ⃗ is the
t
source term, while M = diag(i− , i+ ) . Expressions for B and Δ⃗ are given in [13]. The
general solution of the inhomogeneous equation (5) is given as the sum of the particular
solution of the inhomogeneous equation and general solution of the homogeneous equation:
⃡(0, 𝜏0 )𝐶⃗(𝜏0 ) = ∫𝜏0 𝑃
−𝐶⃗(0) + 𝑃 ⃡⃗−1 Δ⃗(𝜏, 𝜇0 )𝑑𝜏,
⃡(0, 𝜏)𝑀 (6)
0
where
⃡(𝑡, 𝜏) ≡ 𝑒 𝐵⃡(𝜏−𝑡)
𝑃 (7)
is the solution to the homogeneous equation. 𝑃 ⃡(𝑡, 𝜏) is referred to as the propagator and
relates the radiance fields at spatial coordinates t and . That means, if a radiance field at
is known, it can be ‘propagated’ to the point t by applying the propagator (7). However,
our original problem (3) is a two-point boundary value problem, i.e. at each boundary only
a half of boundary conditions (for one hemisphere) is provided. Therefore, solution (6)
Methods Calculating the Slab Radiance Factor 5
based on the propagators cannot be readily applied.
To solve the two-point boundary value problem, the eigenvalue decomposition of the
matrix exponential in Eq. (7) should be performed. Then, applying the so called scaling
transformation [15] and taking into account the boundary conditions, a robust solution of
Eq. (6) can be obtained, which can be written in the form of scatters (see detailed derivations
in [13]):
𝐶⃗ (0) 𝐹⃗ ⃡
𝑅 ⃡− 𝐶⃗+(0)
𝑇
[ − ] = [ −] + [ − ][ ], (8)
𝐶⃗+ (𝜏0 ) 𝐹⃗+ ⃡+
𝑇 ⃡+ 𝐶⃗− (𝜏0 )
𝑅
where 𝐹⃗± are the source vectors, which correspond to the intrinsic radiation of the layer,
while 𝑅⃡ and 𝑇⃡ are the reflection and transmission matrices incorporating the discrete values
of reflection and transmission coefficients, respectively. Note that 𝐹⃗± , 𝑅 ⃡ and 𝑇⃡ are
expressed in terms of eigenvectors of the layer matrix (explicit expressions can be found,
e.g. in [16]). Obviously, the column-vectors 𝐶⃗+ (0), 𝐶⃗− (𝜏0 ) in Eq. (8) correspond to the
radiances, sticking the slab, and are defined by the boundary conditions. The column-
vectors 𝐶⃗− (0), 𝐶⃗+ (𝜏0 ) correspond to the transmitted and reflected radiances. Thus, Eq.(8)
relates the incoming and outgoing radiances and, thus, can be considered as a generalization
of the radiance factor (2)
Equation (8) can be transformed back into the propagator-like form suited for one-point
boundary value problem:
⃡ −𝑅
𝑇 ⃡− 𝑇
⃡−−1 𝑅
⃡+ ⃡− 𝑇
𝑅 ⃡−−1
𝐶⃗(𝜏0 ) = [ + ] 𝐶⃗(0) + 𝐹
⃡. (9)
⃡ −1 ⃡
−𝑇− 𝑅+ 𝑇⃡−−1
In [17] such kind of transformation was referred to as the stellar product. The entries of
scatter matrices satisfy Riccati equation, which correspond to the one-point boundary value
problem. Note that Eqs. (8)-(9) are the strict analytical solution to Eq. (3) in the form of
brightness coefficient matrices.
Two comments are in order. First, following Kolmogorov [18], the physical laws are
established by discrete measurements, and then formulated in term in a continuous setting.
One can argue that nowadays it is more productive and efficient to formulate physical laws
(at least in light engineering) directly in a discrete form. Secondly, the accuracy of the
solution is determined by the number of ordinates used. The use of many ordinates places
high demands on computer power. As the solution is based on the matrix multiplication,
the matrix inversion and the eigen decomposition, the computational time increases faster
as N2. Therefore, engineering applications require high-performance approximate solution
techniques, which would be less demanding on computational power.
6 V. Budak et al
3 Radiance field inside a slab
After the solution of the discrete RTE, taking into account the boundary conditions, the
two-points boundary value problem is transformed into two one-point boundary value
problems (with initial conditions) that allows using each of them to determine the radiance
field inside the slab [19]:
• through the field on the upper boundary
𝑈 ⃡−1 𝐶⃗(0) + 𝑒 −Γ⃡𝜏 ∫𝜏 𝑒 ⃡Γ𝑡 𝑈
⃡−1 𝐶⃗(𝜏) = 𝑒 −Γ⃡𝜏 𝑈 ⃡⃗−1 Δ⃗(𝑡, 𝜇0 )𝑑𝑡,
⃡−1 𝑀 (10)
0
• through the field on the lower boundary
𝑈 ⃡−1 𝐶⃗(𝜏0 ) − 𝑒 −Γ⃡𝜏 ∫𝜏0 𝑒 ⃡Γ𝑡 𝑈
⃡−1 𝐶⃗(𝜏) = 𝑒 ⃡Γ(𝜏0−𝜏) 𝑈 ⃡⃗−1 Δ⃗(𝑡, 𝜇0 )𝑑𝑡.
⃡−1 𝑀 (11)
𝜏
Note that the same function is on the left side of both formulae (10) - (11), and the
expressions on the right sides are different ones. In this case, expression (10) contains only
positive exponents for the upper hemisphere of sighting directions, and expression (11)
contain the same for the lower one. Since left sides are equal, then right sides are also equal.
That allows us to take the expression for the lower hemisphere from (10), and for the upper
hemisphere from (11):
⃡
⃡−1 𝐶⃗ = [𝑒
Γ−(𝜏0 −𝜏)
u11 𝐶⃗− (𝜏0 ) + 𝑒 ⃡Γ−(𝜏0−𝜏)⃡
⃡ u12 𝐶⃗+ (𝜏0 )
𝑈 ⃡+ 𝜏
]+
𝑒 −Γ
u21 C⃗− (0)
⃡
𝜏
− (∫𝜏 0 𝑒 ⃡Γ𝑡 𝑈 ⃡⃗−1 Δ⃗(𝑡, 𝜇0 )𝑑𝑡)
⃡−1 𝑀
⃡𝜏
−Γ −
𝑒 [ 𝜏
]. (12)
(∫0 𝑒 ⃡Γ𝑡 𝑈 ⃡⃗−1 Δ⃗(𝑡, 𝜇0 )𝑑𝑡)
⃡−1 𝑀
+
The obtained expression contains only negative exponents. Therefore, the radiance field
calculation does not present any problems for the arbitrary slab optical depth. Substituting
the expression for the source function in (6) and calculating the corresponding integrals in
(12), we get the resultant expression for the radiance field inside the slab:
⃡
Γ (𝜏 −𝜏)
⃡−1 𝐶⃗(𝜏) = [𝑒 − 0
𝑈 0 ] 𝐺⃗ (𝜏 ) + 𝐹⃗ (𝜏),
0 (13)
0 𝑒 −Γ⃡+ 𝜏
where
𝜏 𝜏
−𝑑𝑘𝜇 −𝜇
2𝑘+1 𝑑𝑘 𝑒 0 𝑒 0
𝐹⃗ (𝜏) ≡ ∑𝐾 ⃡−1 ⃡ ⃡⃗−1 ⃗𝑚
𝑘=𝑚 4𝜋 𝑞𝑘 𝑑𝑖𝑎𝑔 ( 𝛾 𝜇 −𝑑 − 𝛾 𝜇 −1) 𝑈 (1 − 𝜇0 𝑀 )𝑄𝑘 , (14)
𝑖 0 𝑘 𝑖 0
⃡ 𝐶⃗ (𝜏 ) + 𝑢
𝑢 ⃡12 𝐶⃗+ (𝜏0 ) (𝐹⃗ (𝜏0 ))−
𝐺⃗ (𝜏0 ) ≡ [ 11 − 0 ]−[ ]. (15)
⃡21 𝐶⃗− (0)
𝑢 (𝐹⃗(0))+
Methods Calculating the Slab Radiance Factor 7
The radiance angular distributions at some points of the turbid medium are presented in
Fig. 1 as polar diagrams.
Fig. 1. Normalized radiance angular distribution inside the slab: aerosol, t0 = 10, L = 1.0, q0 = 60º.
Red line – t=0.01t0, green – t=0.25 t0, blue – t=0.50t0, magenta – t=0.75t0, black – t=0.90t0.
4 Analysis of approximate methods
Most of the approximate solution techniques were derived in the XX century, when
there was a lack of computational power even to perform summation across thousands of
Legendre expansion coefficients. Even the methods outlined below are not sufficiently
accurate, they are several orders of magnitude faster than the rigorous discrete ordinate
solution. Nevertheless, as we will see later, the ideas behind these approximate techniques
can be reused in the framework of state-of-the-art computational codes for engineering
applications.
4.1. Two-Stream Approximation
A natural approximation based on the theory outlined in Section 2 is the two-stream
approximation. It can be derived by setting N to 1. Thus, we have one discrete ordinate per
hemisphere. Physically the problem (3) is reduced to the case in which we are interested
only in irradiances of upwelling and dowelling radiation. Despite that fact that the radiance
distribution across the polar hemispheres is not considered, the radiant energy is preserved.
Note that such an approximation was formulated by Schuster [20] and Schwarzschild
[21] in the framework of astrophysical radiative transfer. Mathematically, the feature of the
two-stream approximation is that the eigenvalue problem can be solved analytically, and
the solution can be expressed in a closed form, thereby yielding high computational speed.
8 V. Budak et al
4.2. Diffuse approximation
We present all functions in the RTE in the boundary value problem (3) as a spherical
function expansion for the case P1 of the spherical harmonics method:
( )
1
1
Lr (, ˆl ) = C00 () Y00 (ˆl ) + C1m () Y1m (ˆl ) = E0 () + 3E()ˆl , (16)
m =−1 4
( )
1
1
(, ˆl ) = s00 () Y00 (ˆl ) + s1m () Y1m (ˆl) s0 () + 3s()ˆl , (17)
m =−1 4
2l+ 1
x(ˆl ˆl ) = xl Pl (ˆl ˆl ) , (18)
l = 0 4
where 𝐸0 (𝜏) ≡ ∮ 𝐿𝑟 (𝜏, 𝒍̂)𝑑𝒍̂ , 𝐸(𝜏) ≡ ∮ 𝐿𝑟 (𝜏, 𝒍̂)𝒍̂𝑑𝒍̂ are spatial irradiance and light vector,
respectively, and similarly s0(), s() for the source function;
2k + 1 (l − n)! n
Ykm (ˆl ) = Pl () eim , Pln () are the associated Legendre polynomials.
4 (l + n)!
Ckm , skm are coefficients of radiance expansion and source functions by spherical functions.
After substituting (16)-(18) in RTE, we consistently integrate the resulting expression
by the full solid angle and with the weight l̂ . Taking into account the orthogonality of
spherical functions we obtain a connected system of two equations [22]:
(1 − ) E0 () + E () = s0 (),
1 (19)
E0 () + (1 − x1 )E () = s(),
3
that corresponds to the P1-approximation of the spherical harmonics method [22].
Under the assumption that the medium is homogeneous e=const, we can calculate the
divergence from the second equation of the system (19), express E () from the first and
substitute it into the second, which finally leads to the equation [22]:
1
E0 () − E0 () = () , (20)
a2
s0 () s()
where a = (1 − )(1 − x1 ) , () − + is the source function.
3(1 − ) 3a 2
Equation (20) admits an analytical solution.
Methods Calculating the Slab Radiance Factor 9
4.3. Invariant imbedding and matrix operator method
The invariant imbedding method was proposed by V. A. Ambartsumian in the first half of
the 20th century [23]. This method is based on the fact that the addition to the semi-infinite
medium of a layer with an arbitrary thickness and the same optical properties does not
change the reflected light intensity. As the layer can be of arbitrary thickness, we add an
infinitely thin layer, so that in this layer the multiple scattering collisions can be neglected.
The resulting equation is formulated for the reflection function. Note, that here we deal with
the one-point boundary value problem. In the discrete ordinate setting, the Ambartsumian
equation is transformed into matrix Riccati equation [24].
A consequence of the invariant imbedding principle is the so-called matrix operator
method. Considering Eq. (8) for two adjacent layers, a similar equation can be derived for
an efficient layer, which has the same optical properties as the system of two-layers [13].
Note that this principle can be used for the adding-doubling algorithm, in which the
reflection and transmission matrices are found for a very thin layer using a certain
approximation, and then by multi-fold doubling the corresponding entries for the desired
optical thickness can be obtained.
4.4. Method of iterations
The radiation transfer equation was first formulated by Chwolson [25] and Lommel [26] in
photometry for calculating the transmission and reflection of milk glasses. The equation
was expressed in an integral form, and the solution was proposed by the iteration method.
In the method of iterations, the diffuse radiance can be expressed as follows [27]:
Ns
L = Ls , (21)
s =0
where Ls refers to the s-fold scattered radiance, while Ns the maximum scattering order
considered. The expansion (21) is called Neumann series, and the computational technique
based on Eq. (21) is called successive order of scattering (SOS). The computations of Ls
can be organized iteratively. Let us consider the radiative transfer equation from Eq. (3),
L(, ˆl )
+ L(, ˆl ) = J (, ˆl ) (22)
in which the multiple scattering term is denoted by the source function J (, ˆl ) , i.e.
𝛬
𝐽(𝜏, 𝒍̂′ ) = 4𝜋 ∮ 𝐿(𝜏, 𝒍̂′ )𝑥(𝒍̂, 𝒍̂′ )𝑑𝒍̂′ . (23)
At the zero-th iteration we set
L0 (, ˆl ) = e− / 0 (ˆl − ˆl 0 ) (24)
10 V. Budak et al
and evaluate J 0 (, ˆl ) according to Eq. (23). Now the right hand side of Eq. (22) is known
and Eq. (22) can be integrated, yielding for the diffuse downwelling and upwelling radiation
−
− d
Ls (, ˆl ) = J s (, ˆl )e
(25)
0
and
− L −
− L −
d
Ls (, ˆl ) = Ls (0 , ˆl ) e + J s (, ˆl ) e
. (26)
( )
respectively. Then Ls , ˆl is used to find J s +1 (, ˆl ) from Eq. (23) as
𝛬
𝐽𝑠+1 (𝜏, 𝒍̂′ ) = 4𝜋 ∮ 𝐿𝑠 (𝜏, 𝒍̂′ )𝑥(𝒍̂, 𝒍̂′ )𝑑𝒍̂′ . (27)
The iterative scheme (24)-(27) converges to a ’true’ solution. However, if the medium
is transparent, the convergence is awfully slow. Also, note, that the at low number of
iterations the energy is not conserved.
4.5. Small-angle approximation
For sighting angles close to incident angles (i.e. 0 ), the small angle approximation
can be used to compute the transmitted radiance for the case of extended in the forward
direction phase functions [28]. The initial boundary value problem is then transformed into
𝜕𝐿(𝜏,𝒍̂) 𝛬
𝜇0 + 𝐿(𝜏, 𝒍̂) = ∮ 𝐿(𝜏, 𝒍̂′ )𝑥(𝒍̂, 𝒍̂′ )𝑑𝒍̂′ ,
{ 𝜕𝜏 4𝜋 (28)
𝐿(𝜏, 𝒍̂)|𝜏=0,𝜇>0 = 𝛿(𝒍̂ − 𝒍̂0 ).
Note that is substituted by 0 before the derivative. The advantage of the small-angle
radiative transfer equation is that it can be solved analytically. In particular, Eq. (28) can be
solved by using the spherical harmonics method, i.e. expanding the radiance into Fourier
cosine series and then Legendre series. Note that on the left hand side of Eq.(28) we have
0 dL / d instead of dL / d , and thus, all terms in Eq. (28) have the same expansion
kernels. The solution reads as follows
1 ( n − m )! − 0 (1−n ) m 1
LSAA (, , ) =
4 m = 0 n = m
( 2 − m0 )( 2n + 1)
( n + m )!
e Pn () Pnm ( 0 ) cos m, (29)
where m 0 is 1 if m=0 and zero, otherwise, n is the Legendre expansion coefficients of
the phase function, while Pnm are the associated Legendre polynomials. The small-angle
Methods Calculating the Slab Radiance Factor 11
solution gives a Dirac delta function at the upper boundary and a peaked function inside the
scattering medium. The amplitude of the peaked function decreases with increasing the
optical depth, while its width increases with increasing the optical depth. Fig.2 illustrates
the small-angle solution. As we can see, it agrees well with the discrete ordinate solution
for angles close to the incident angle.
Fig. 2. Comparison between the discrete ordinate solution, the small-angle approximation, and the
single scattering approximation. The computations are performed for the Henyey-Greenstein phase
function with the asymmetry parameter 0.9. The single scattering albedo is 0.99. The angle of
incidence is 60 degrees. (left) the optical thickness is 0.3. (right) The optical thickness is 2.0.
5 Synthetic iterations
Keeping in mind the advantages and limitations of the exact discrete ordinate approach and
approximate solution techniques, we can formulate a hybrid approach based on synthetic
iterations [29-30]. Originally such an approach was formulated in nuclear physics. The
computations of the radiance field are performed in two steps. At the first step, an
approximate solution is found. It should preserve the energy balance (i.e. the integral
quantity) while we are not interested in a detailed angular dependency of the radiance field.
At the second step, the solution is refined by using the iteration technique similar to
Eqs. (25)-(26). However, here at the first iteration we already have a good approximation,
while in the original formulation of iterative method, the initial solution based on the single
scattering approximation has a large error.
A computationally efficient method of synthetic iterations applied to Eq. (3) consists in
the following steps:
1. the anisotropic part of the radiance field (together with the singular component) is
expressed by using the small-angle approximation;
2. the remaining regular part of the radiance is computed in the two-stream
approximation, which preserves the energy balance;
3. associating the current solution with the source function in Eqs. (25)-(26), the solution
12 V. Budak et al
is refined for given viewing angles by applying iterative scheme (25)-(26). To perform
integration over t, we make use of expressions for the radiance inside the layer (see Eqs
(13)-(15)).
Note that the synthetic iteration method inherits advantages from both approximate and
discrete ordinate solution techniques.
6 Summary
Analytical discrete transfer theory [9] allowed us to show the relationship of all approximate
methods for determining the radiance factor, its accuracy, and applicability. The most
significant method for improving the accuracy of approximate solution methods is to isolate
the anisotropic part (4) from the solution based on the small-angle approximation
(subsection 4.5). The regular part of the solution can be determined by any approximation,
which provides an accurate definition of the reflection coefficient. To improve determining
the accuracy of the angular distribution of the radiance factor (BRDF), the method of
synthetic iterations (section 5) is the most effective. It seems that this approach solves all
the problems of determining the radiance factor of the real-world objects.
References
1. Kurt, M., Edwards, D.A.: Survey of BRDF models for computer graphics. Computer Graphics
(ACM) 43(2), 4 (2009)
2. Lai, Q., Liu, B., Zhao, J., Zhao, Z., Tan, J.: BRDF characteristics of different textured fabrics
in visible and near-infrared band. Optics Express, 28(3), P.3561 (2020)
3. The optics of Ibn Al-Haytham. Books I—III: On Direct Vision. Translated with introduction
and commentary by A.I. Sabra. The Warburg institute university of London, London (1989)
4. Fabry, C.: Introduction Générale à la Photométrie, Editions de la revue d’optique théoretique et
instrumentale. Paris (1927).
5. ILV: International lighting vocabulary, CIE S 017/E:2011. Commission Internationale de
L'Eclairage, Vienna (2020)
6. Pokrowski, G.I.: Zur Theorie der diffusen Lichtreflexioni. Zeit. f. Physik, 30, 66 (1924)
7. Torrance, K.E., Sparrow, E.M.: Theory for off-specular reflection from roughened surfaces.
JOSA, 57(9), 1105-1114 (1967).
8. Basov, A.Y., Budak, V.P.: Mathematical model of a surface radiance factor. CEUR Workshop
Proceedings, 29th International Conference on Computer Graphics and Vision, GraphiCon
2019; Bryansk; Russia; 23 - 26 September 2019, 2485, 43-47 (2019)
9. Afanas'ev, V.P., Basov, A.Y., Budak, V.P., Efremenko, D.S., Kokhanovsky, A.A.: Analysis of
the discrete theory of radiative transfer in the coupled "ocean-atmosphere" system: current
status, problems and development prospects. Journal of Marine Science and Engineering, 8(3),
202 (2020). https://doi.org/10.3390/jmse8030202
10. Chandrasekhar, S.: Radiative transfer. Dover publications, inc. New York (1950)
11. Milne, E.A: The reflection effect of the eclipse binaries. Mon. Not. Roy. Astrophys. Soc.,
LXXXVII, 43-49 (1926).
Methods Calculating the Slab Radiance Factor 13
12. Bohren, C.F., Huffman, D.R.: Absorption and Scattering of Light by Small Particles. Wiley
(Apr 1998). https://doi.org/10.1002/9783527618156
13. Budak, V.P., Klyuykov, D., Korkin, S.V.: Complete matrix solution of radiative transfer
equation for pile of horizontally homogeneous slabs. J Quant Spectrosc Radiat Transfer 112(7),
1141–1148 (2011). https://doi.org/10.1016/j.jqsrt.2010.08.028
14. Budak, V.P., Efremenko, D.S.; Shagalov, O.V.: Efficiency of algorithm for solution of vector
radiative transfer equation in turbid medium slab. J. Phys. Conf. Ser. 369, 012021 (2012).
https://doi.org/1088/1742-6596/369/1/012021
15. Karp, A.H.; Greenstadt, J.; Fillmore, J.A.: Radiative transfer through an arbitrarily thick,
scattering atmosphere. J. Quant. Spectrosc. Radiat. Transf. 24, 391–406 (1980). doi:
10.1016/0022-4073(80)90074-6
16. Efremenko, D.S., Molina García, V., Gimeno García, S, Doicu, A.: A review of the matrix-
exponential formalism in radiative transfer. Journal of Quantitative Spectroscopy and Radiative
Transfer 196, 17–45 (2017). https://doi.org/10.1016/j.jqsrt.2017.02.015
17. Plass, G.N., Kattawar, G.W., Catchings, F.E.: Matrix Operator Theory of Radiative Transfer.
Appl. Opt. 12, 314-326 (1973)
18. Kolmogorov, A.N.: Combinatorial foundations of information theory and the calculus of
probabilities. Russ. Math. Surv. 38, 29–40 (1983).
https://doi.org/10.1070/RM1983v038n04ABEH004203
19. Budak, V.P., Shagalov, O.V., Zheltov, V.S.: Numerical radiative transfer modeling in turbid
medium slab. Proc. of SPIE, 9292, 92920Y (2014).
20. Schuster, A.: The influence of radiation of the transmission of heat. Phil.Mag., 5, 243–257
(1903)
21. Schwarzschild, K.: Über das Gleichgewicht der Sonnenatmosphäre. Nachr. Konig. Gesel. der
Wiss., Gottingen, 1,41–53 (1906).
22. Davidson, B.: Neutron Transport Theory. Clarendon, Oxford, UK, 1957
23. Ambartsumian, V.A.: On the problem of diffuse reflection of light. J. Phys. USSR, 8, 65–75
(1944)
24. Flatau, P.J., Stephens, G.L.: On the Fundamental Solution of the Radiative Transfer Equation.
J. Geophys. Res., 93, 11037–11050 (1988)
25. Chwolson, O.D.: Grundzüge einer matimatischen Theorie der inneren Diffusion des Licht. Bull.
Acad. Imp. Sci. St. Petersbourg. 33, 221–256 (1889)
26. Lommel, E.: Die Photometric der diffusen Zuruckwerfung. Sitzber. Acad. Wissensch.
Munchen, 17, 95–124 (1887)
27. Lenoble ,J., Herman, M., Deuzé, J.L., Lafrance, B., Santer, R., Tanré, D.: A successive order of
scattering code for solving the vector equation of transfer in the earths atmosphere with aerosols.
Journal of Quantitative Spectroscopy and Radiative Transfer, 107(3), 479–507 (2007)
28. Kokhanovsky, A.A.: Small-angle approximations of the radiative transfer theory. Journal of
Physics D: Applied Physics, 30(20), 2837–2840 (1997)
29. Budak, V.P., Kaloshin, G.A., Shagalov, O.V.; Zheltov, V.S.: Numerical modeling of the
radiative transfer in a turbid medium using the synthetic iteration. Opt. Exp. 23, A829–A840
(2015). https://doi.org/10.1364/OE.23.00A829
30. Budak, V.P., Zheltov, V.S., Lubenchenko, A.V., Freidlin, K.S., Shagalov, O.V.: A fast and
accurate synthetic iteration-based algorithm for numerical simulation of radiative transfer in a
turbid medium. Atmos. Ocean. Opt. 30, 70–78 (2017).