1. Introduction

Understanding practical multiple scattering phenomenon has been long considered a challenging task, especially for the intermediate regime where the optical depth (OD) ranges between 1 and 10 [1]. However, there is a significant need to use multiple scattering measurements for reconstruction purposes, for instance, performing tissue imaging or rebuilding droplet diameter spatial distribution. A robust method to model the phenomenon would allow light scattering diagnostic techniques to become useful for variety of applications where they are currently not practical. For example, optical methods could be more commonly used to examine dense liquid sprays used in diesel engine combustion, if errors induced by multiple scattering can be properly addressed [2].

Typical solutions for the multiple scattering problem include solving a set of radiative transfer equations (RTEs) [3] or using stochastic methods such as Monte Carlo simulations [1, 4, 5]. These methods can yield accurate predictions of multiple scatterings with preset parameters, such as scattering phase functions, OD, incident light angle and strength distribution, etc. However, these methods would consume a considerable amount of computational resources. Further, their complex or stochastic nature has prohibited any successful attempts at multiple scattering inversions.

Beyond solving RTEs and Monte Carlo simulations, there have also been many studies focused on developing reduced models for multiple scattering processes. Such efforts include Random Walk theorem [6, 7], empirical predictions [8], and adding-doubling method [9], etc. These methods can usually provide simplified, analytical expressions for predictions of crucial observations such as total transmission, average cosine, and traveled path length distributions. However, reduced models also have drawbacks as most utilize simple phase functions as a function of g (anisotropy factor) or assuming isotropic scattering. Such simplifications, while giving close results for some of the averaged observations, cannot distinguish similar phase functions, for instance, phase functions with the same g but completely different probability density functions (PDFs). They also suffer from limitations including the inability to deal with anisotropic scattering, absorbing mediums, and/or non-uniform OD/phase function distributions.

Such needs and limitations motivate the authors to use improved models to simulate multiple scattering processes and for multiple scattering reconstructions. Previous work by the authors [10] suggested the prospective of utilizing discrete Markov Chain approximation to model photon multiple scattering. We proved that the performance of the Markov Chain model was superior to Monte Carlo simulation in solving scattering through turbid slabs. However, the previous work was limited in that scattering was assumed to be isotropic, which was not true for many cases.

Extending the previous study, this work investigates Markov Chain approximation with an anisotropic scattering assumption. More specifically, this work focuses on analyzing the transmitted/reflected photon angular distribution with the developed Markov Chain model. Although sharing some common purposes, the Markov Chain model has many additional capabilities compared with Monte Carlo simulations. The semi-analytical form of the Markov Chain model can promote more thorough parameter studies for multiple scattering processes. Furthermore, the matrix form of Markov Chain model will enable advanced inversion algorithms based on matrix manipulations or more effective stochastic optimization algorithms, as will be discussed in the following sections.

2. Markov Chain approximation for anisotropic scattering

The Markov Chain method, which is used to infer statistical predictions by utilizing the transition probability from one state to another, has found numerous applications in various fields [11]. By defining the transition matrix P, in which each entry represents the probability that the random process transits from one state to another, the Markov Chain theorem predicts how the random process will evolve. Markov Chain approximation requires the random process to be “memoryless”, which is; the probability distribution of the next state depends solely on the current state. Photon multiple scattering is a process that satisfies the requirement of Markov Chain theorem, since the new propagation direction and distance of the photon before the next scattering event is only dependent on the current scattering event. Successful Markov Chain modeling for photon propagation through isotropic turbid slabs has been demonstrated in a previous work [10] where more details regarding Markov Chain approximation can be found.

To properly address anisotropic scattering, one additional category of state has to be considered in the Markov Chain model. The schematic of the Markov Chain process for anisotropic scattering is shown in Fig. 1. In the schematic, the turbid slab is discretized into multiple layers with the same thickness Δz and we assume the optical properties within each laser are uniform. m and n are arbitrary layer indexes and z_m and z_n represents the distance between the layer and the incident plane for layer m and n, respectively. We assume two consecutive scattering events take place in the mth and nth layer. θ_i and θ_j are photon propagation angles (zenith angles) and φ is the azimuth angle. It is assumed that φ = 0 degree for layer m due to symmetry. α is defined as the angle between the two propagation vectors.

 

Fig. 1 Schematic of multiple scattering modeled by Markov Chain model.

Download Full Size | PDF

In order to solve the underlying multiple scattering problem, the transition matrix must be built. In this work we consider number of scattering events for each photon as the “event” in Markov Chain theorem and define (z, θ) as Markov Chain states. The goal then is to find the transition probability for photons to propagate from state (z_m, θ_i) to state (z_n, θ_j) between two adjacent scattering events. The following relationships be readily derived:

3. Markov Chain approximation computation

In order to verify the fidelity of the proposed Markov Chain method for anisotropic scattering calculations, we performed computations with different scattering medium properties. The results were then compared with Monte Carlo simulations. In this study we only used Monte Carlo simulation as the reference, however such comparison will still substantially verify the Markov Chain model since the same Monte Carlo codes have been verified and validated by other numerical methods and experiments [1, 5, 10]. Figure 2 represents the phase functions (probability density function, PDF) used in this study. The phase functions were calculated assuming Mie scattering with a resolution of 1 degree and a uniform particle diameter of 1.2 μm, incident light wavelength of 532 nm, and ambient refractive index of 1.0. The only difference between the two phase functions is that for the first phase function, scattering particle refractive index was set to be 1.33 and for the second phase function it was 1.57, which yields an anisotropy factor g of 0.7649 and 0.6653, respectively. Unpolarised light was assumed. In this study, the turbid slab was discretized into 100 layers and the propagation angle was discretized into 180 intervals with a resolution of 1 degree, which makes P a 18,000 by 18,000 matrix.

 

Fig. 2 Phase functions (PDF) selected for this study.

Download Full Size | PDF

Before we proceed with non-uniform turbid slabs, we examine the performance of Markov Chain approximation with other methods in an approach similar with [12]. To be consistent with previous work [10], OD was selected to be 4.242. Cases with and without absorption were both considered. As can be seen in Table 1, random walk theory yielded the worst performance. Adding-doubling method predicted the total transmission at an error ~1% and Markov Chain approximation about 0.1% with Monte Carlo simulation as the benchmark, thus verifying the proposed Markov Chain model. Figure 3 exhibits the comparison between Markov Chain approximation and Monte Carlo simulation for transmitted/reflected photon angular distribution. In this simulation, the thickness of the slab was set to be 1 cm and the extinction coefficient was set to be 10 cm⁻¹. We assumed no absorption, thus the total OD was 10.0. We considered two different cases as follows: In case 1, the scattering phase function in layers 1-25 and layers 76-100 follows phase function 2, and the phase function in layers 26-75 follow phase function 1; while for case 2 the phase functions were swapped. The Monte Carlo simulation sent 500 million photons. For simplicity, the transmitted/reflected angular distributions were presented in the same plots, with θ between 0 to 90 degrees for transmitted photons and 90 to 180 degrees for reflected photons. It took Monte Carlo 8 hours and Markov Chain about 10 minutes to perform the same simulation. From this, it is evident that the computational capacity of Markov Chain approximation is satisfactory compared to typical methods. It can be observed that Markov Chain approximations generated accurate transmitted/reflected angular distributions for both cases. Furthermore, the simulations show that the arrangement of turbid slab layers could change the resultant angular distributions. As seen, for θ within 0 to 90 degrees (transmitted photons), the angular distribution predictions for both cases are quite close. The reason for this phenomenon is that since the total OD is relatively large, the transmitted photons have scattered sufficient times so that the scatterings average out the impact of the phase function. In comparison, θ angular distributions within 90 to 180 degrees (reflected photons) for different phase function distributions have a substantial difference, indicating reflected photon angular distribution can be used as the input for diagnostics purpose.

Table 1. Comparison of Total Transmittance with Different Methods

View Table

 

Fig. 3 Transmitted and reflected photons with different phase function distributions.

Download Full Size | PDF

Figure 4 shows the comparison between Markov Chain approximation and MC simulation for transmitted/reflected photon exiting the domain of interest that went through exactly 1, 2, and 3 scattering events, respectively. As can be observed, Markov Chain approximation matches Monte Carlo simulation angular distribution well in that Markov Chain approximation successfully captured all the details of the angular distribution profiles, both in shape and magnitude. It is worth noting that Markov Chain approximation agrees with Monte Carlo simulation with any given number of scattering events and only scattering events up to 3 were plotted for simplicity. Furthermore, Markov Chain approximation has successfully predicted the scattering events with very simple matrix forms and no matrix inversion calculation is required, comparing to total transmitted/reflected photon solutions. Such features could make the reconstruction process less computational expensive and enable fast reconstructions with either deterministic algorithms or stochastic algorithms.

 

Fig. 4 Transmitted/reflected photon angular distribution for photons exiting with different times of scatterings. Markov Chain theorem angular distribution predictions for # of scattering = 1, 2, and 3 are P’ × R, P’ × P × R, and P’ × P² × R, respectively.

Download Full Size | PDF

Finally, Fig. 5 shows Markov Chain approximations with a scattering and absorbing medium. The scattering coefficient was set to be 9.5 cm⁻¹ and the absorption coefficient was set to be 0.5 cm⁻¹, both uniform across the turbid slab. Same phase function combinations were used as in Fig. 3. As evident in Fig. 5, the Markov Chain approximation faithfully reproduced the same transmitted/reflected angular distributions as Monte Carlo simulation did. It is worth noting that for all the simulations presented in this work, the maximum absolute relative error is mostly less than 1% under the current simulation conditions, which proves the capability of Markov Chain approximation in handling one-dimensional anisotropic scattering problems. Furthermore, absorbing medium changed the intensity and profile of the angular distribution for reflected photons, as evident by the transmission peak around θ = 170 degrees for case 1 (same extinction coefficients were used as in Fig. 3), therefore transmitted angular distribution can be used to infer absorbing coefficient too.

 

Fig. 5 Markov Chain approximation for an absorbing medium. The scattering coefficient and absorption coefficient were set to be 9.5 cm⁻¹ and 0.5 cm⁻¹, respectively.

Download Full Size | PDF

4. Discussion and conclusion

This paper has presented a novel method of utilizing a Markov Chain approximation to simulate photon multiple scattering through a turbid slab with anisotropic scattering assumptions. The purposed Markov Chain approximation is considered superior to Monte Carlo simulation in handling such problems because Markov Chain method reformats the problem into matrix form, allowing feasible input parameter studies and inversion studies for this complex problem. To the authors’ knowledge, the purposed Markov Chain method is one of the few methods that yields accurate predictions with minimum simplifications besides Monte Carlo simulations and solving radiative transfer equations.

Furthermore, as mentioned earlier, Markov Chain approximation provides a viable approach in reconstructing the optical properties of the turbid slab, such as optical depth, phase function, and their distributions. Traditionally Markov Chain method has been used to predict the steady state (in this study transmitted/reflected photons) with a known transition state. While recently there have been attempts of reconstructing the transition state with the steady state as priori, given the transition state is constrained in some way [13]. With the matrix form of Markov Chain approximation, one could easily rewrite the cost function as Q_simualtion-Q_measurement, then minimize the cost function to inverse desired quantities in a fashion similar as in [14].

It is also worth noting that there have been some successful experimental investigations [15] in reflected photon angular distribution measurements. In this previous study, the reflected photon angular distributions were recorded as a function of time. Although the distance (time) that a photon has traveled is not directly correlated to how many times the photon has scattered before the photon exit, it is likely that such relationship can be built to enable reconstructions in obtaining medium optical properties. Another potential solution to incorporate the experimental data is to improve the current Markov Chain approximation, so that time (distance) is also considered as a factor in building the transition matrix. With this method, proper optimization algorithms can be selected for practically feasible reconstructions, a capability that has been long desired by the optical research community.