1. Introduction

Multiple scattering (due to Mie scattering) represents a fundamental problem with applications in a wide spectrum of practical systems, ranging from imaging through biological tissues [1], remote sensing in the atmosphere [2], and development of laser diagnostics for dense sprays [3]. Consequently, this problem has attracted a considerable amount of research efforts from different disciplines. In many of these applications, the angular distribution (AD) of the multiply scattered photons plays an important role. For example, the AD plays a critical role in designing the imaging system involving multiple scattering, such as the selection of the lens [4], estimation of the detection limit and resolution [5, 6], and the application of angular filtering [7]. Furthermore, practical applications typically involve an optimization facet, which requires the AD to be computed a great many times. Therefore, it is highly desirable to have a simple (preferably closed-form) method to calculate the AD efficiently.

Fundamentally, the multiple scattering problem requires one to solve the wave equation in a system of a large number of scatterers, which is extremely difficult [8]. Under certain assumptions, the problem can be simplified and the so-called radiative transfer equation (RTE) can be derived [9]. However, even the solution of the simplified RTE is difficult under most non-ideal conditions (e.g., realistic geometry). Thus various approximate solution techniques, both theoretical and numerical, have been developed [4, 8, 10–12]. The theoretical techniques seek analytical solutions to the RTE under certain approximations [13]. For instance, under large optical thickness, the photon diffusion model [4, 10] approximates the multiple scattering problem as a diffusion problem. Another example, the random walk model [11, 14] approximates the multiple scattering problem as a random walk process. These models result in analytical expressions that can be used to predict the spatial and temporal distribution of multiply scattered photons. However, these analytical expressions can be regarded as semi-closed-form (at lease from the perspective of practical applications). For example, these expressions typically involve special functions and summation of infinite series [10, 11, 15], which complicates implementation. Furthermore, much of the analytical work did not target the AD specifically, and therefore application of the previous results to calculate the AD requires extension work, sometimes considerable extension work.

The numerical techniques solve the RTE computationally and a major limitation is the computation cost. The Monte Carlo (MC) technique provides a notable example of the numerical technique. The MC technique has been demonstrated as a powerful tool, capable of solving multiple scattering problems under realistic conditions and resolving quantities that are difficult to study experimentally [12, 16–18]. However, the MC technique encounters considerable challenges in terms of computational cost, especially at relatively large optical depths. The challenge is especially acute for calculating the AD, because resolving the AD requires receiving sufficient photon counts in a small angular range out of the transmitted photons, which are already a small fraction of incident photons at large optical depth (illustrated by the results in Section 3).

The goal of this paper is therefore to develop a closed-form method to calculate the AD of multiply scattered photons. Closed-form expressions are derived to calculate the AD of photons transmitted from isotropic turbid slabs, based on which, the average cosine and fraction of photons within a given acceptance angle (numerical aperture) can be readily calculated. The results obtained from the method were compared against MC data, and good agreement is achieved over a wide range of optical depths. The method is simple to implement and highly efficient, suitable for application in a range of practical applications. Section 2 describes the method. Section 3 reports the results and discusses the extension of the current method to anisotropic scatterers. Finally, Section 4 summarizes the paper.

2. Description of the method

Figure 1 illustrates the schematic of the multiple scattering problem and the MC configuration used in this work. As shown in Panel (a), incident photons enter an infinitely large turbid slab as a pencil beam. The slab, with thickness Z, contains isotropic scatterers. Note that isotropic slabs in this work mean the scatterers are isotropic and the distribution of the scatterers in the slab is also isotropic (i.e., uniform). All length scales in this work are expressed in units of the mean scattering length, i.e., ${\Sigma }_S^{-1}$, where ${\Sigma }_S^{}$(mm-1) represents the scattering coefficient of the scattering media. The plane on which the incident photons enter is named the Incident Plane (IP), and the opposite plane named the Exit Plane (EP). After a certain number of scattering events, the incident photons exit either via the IP (as reflected photons) or via the EP (as transmitted photons). Here we focus on the AD of the transmitted photons, and extension of the method to calculate the AD of the reflected photons is straightforward.

 

Fig. 1 Schematic of the problem and the configuration of MC simulation.

Download Full Size | PDF

The MC simulation in this work approximates the infinitely large slab by a rectangle slab, whose width and height are both 100 × of its thickness as shown in Panel (a). This approximation is justified by the fact that for all cases simulated, the fraction of photons leaked via the sides is less than 10⁻⁶. This work used a standard MC model as described in [16], and the implementation of the model largely follows [12]. The MC model involves the following four major steps: 1) it sends incident photons into the IP one by one, 2) it determines the scattering direction and pathlength of the incident photon stochastically to calculate the location of the photon after a scattering event, 3) it repeats step 2) until the photon exits the slab, and 4) it registers the location and orientation of each photon that exits. The MC model terminates after a sufficient number of exit photons have been registered to determine the AD with acceptable statistical accuracy, which corresponds to a number of incident photons on the order of 10⁷ in this work. The free scattering pathlength in the MC model was taken to be exponentially distributed, with the mean free pathlength equal to${\Sigma }_S^{-1}$.

A Cartesian coordinate system is defined as shown in Panel (b), where the z-axis is perpendicular to the EP. The origin is defined to be on the EP, and the positive of the z-axis points towards the IP. The scattering angle of a transmitted photon (θ) is defined as the angle formed relative to the negative z-axis after its last scattering event (illustrated by photon 1). Under these definitions, the AD function (Q_T(θ)) of the transmitted photons is:

With the pathlength assumed to be exponentially distributed, it can be shown that:

 

Fig. 2 Panel (a). Comparison of P against MC data. Panels (b) and (c). Comparison of Q_z against MC data under various OD.

Download Full Size | PDF

The rigorous derivation is of Q_z is more complicated. Here we derive Q_z under some simplifications, and then briefly discuss the rigorous derivation; though the results derived under these simplifications are already in reasonable agreement with MC data. The first simplification is to assume that the transmitted photons are always scattered at θ = 0 by the last scattering event (i.e., they all exit along the negative z-axis as photon 2. This appears to be a dramatic simplification, however, its validity on isotropic scatterers have been well demonstrated in previous work where a random walker on a lattice is used to model multiple scattering [11, 14]. Also note that this simplification is only invoked for the derivation of Q_z (one step in the derivation of the AD), and does not imply that the AD is a delta-function. The second simplification assumes that all transmitted photons exit exactly on the EP. With these simplifications, Q_z could be understood as the distribution of photons with respect to z after one scattering event along the positive z-axis if the photons were injected into the slab on the EP. With this understanding, because the photon pathlength is exponentially distributed, Q_z is:

The constant, C in Eq. (1), depends on how the AD is normalized. If the AD is normalized relative to the transmitted photons, i.e. ${\displaystyle {\int }_0^{\pi /2}{Q}_T(\theta )\sin \theta d\theta }=1$, then C = 1. One could alternatively normalize Q_T(θ) relative to the incident photons such that ${\displaystyle {\int }_0^{\pi /2}{Q}_T(\theta )\sin \theta d\theta }$ = fraction of photons that are transmitted (i.e., transmittance). In this case, it can be shown that [11]:

3. Results and discussions

Figure 3 shows the AD of the transmitted photons calculated using Eq. (5) compared to MC data. As shown in Fig. 3, the method developed above agrees with MC data across a wide range of OD. Two interesting observations can also be made: 1) the AD is peaked towards θ = 90° at small OD (Panel (a)), and 2) after OD becomes larger than 5, the AD becomes insensitive to OD (Panel (b)). The first observation can be explained by noting that at small OD (e.g., 0.1), photons undergo fewer scattering events to be transmitted and consequently are more likely to exit at large angles (if they do exit). The second observation can be explained by noting that at large OD, photons undergo a large number of scattering events and their directions become completely randomized, resulting in the AD’s insensitivity.

 

Fig. 3 AD of transmitted photons calculated using Eq. (5) compared to MC data. The symbols represent the MC data (with ballistic photons excluded), the solid lines the calculation using Eq. (5).

Download Full Size | PDF

With the validity of Eq. (5) confirmed, other properties of interest can be readily calculated. For example, the average cosine of the transmitted photons can be expressed as:

 

Fig. 4 Panel (a). Comparison of the average cosine at various OD. Panel (b). Fraction of transmitted photons within given acceptance angles. Panel (c). Error relative to the MC data.

Download Full Size | PDF

Another quantity of interest to the design of imaging optics is the fraction of photons that can be collected. Such fraction can be calculated by${\displaystyle {\int }_0^{\theta a}{Q}_T(\theta )\sin \theta d\theta }$, where θ_a is the acceptance angle defined by the imaging system. For a simple lens, θ_a is determined by the diameter of the lens and the distance between the lens and the EP. Example calculations are shown in Panel (b) and (c) of Fig. 4. As can be seen, the calculation from Eq. (5) agrees well with the MC data (within −1 to 4%). At a θ_a of 1.5°, the fraction of transmitted photons within this angle varies between 4 × 10⁻⁴ to 8 × 10⁻⁴, i.e., on average, one photon exits within θ < 1.5° out of ~1250 to 2500 transmitted photons. Only a fraction of the incident photons are transmitted. Therefore, as mentioned earlier, the MC techniques can be prohibitively expensive when applied simulate the AD of multiply scattered photons, especially at relatively large OD; and the method developed above can overcome such computational limitation.

Lastly, we discuss the simplifications made in the derivation of Q_z, and the extension of the method to anisotropic scatterers. As shown in Fig. 2 (Panels (b) and (c)), the two simplifications made in the derivation of Eq. (4) results in more slow decay of Q_z compared to the MC data. Such discrepancy can be addressed by either empirically fitting the MC data to obtain Q_z, or by removing the simplifications. The former approach is straightforward and our results have shown that an empirically fitted Q_z can improve the agreement for calculations shown in Figs. 3 and 4. The latter approach is closely related to the extension to anisotropic scatterers. Equation (5) is fundamental and applicable for anisotropic scatterers. However, derivation of P and Q_z need to be modified because of the directionality of the scattering for anisotropic scatterers. It has been shown previously that certain scaling laws based on the asymmetric factor of the scatterers can be used to extend isotropic models [11, 14] to anisotropic scatterers, and we plan to apply these scaling laws in the derivation of the AD.

4. Summary

In summary, this paper develops a method for calculating the AD of multiply scattered photons through isotropic turbid slabs. A closed-form method has been developed and validated using MC data over a wide range of conditions. Based on this method, various quantities of interest to multiple scattering can be readily calculated to aid the design of imaging optics. Calculating the AD using numerical techniques encounters significant computational limitations at relatively large OD. These limitations become especially acute when one needs to optimize the design of the imaging optics iteratively, or to solve the inverse scattering problem. The method developed here, due to its simplicity in implementation, can overcome these limitations; and we therefore expect it to be a valuable tool to study multiple scattering.