Numerical Study of Turbid Slab Optical Properties Reconstruction from Multiple Scattering Signals Using Time-Based Markov Chain Model

: The reconstruction of optical properties for opaque mediums is highly desired for medical, atmosphere and aerosol applications. However, the modeling and reconstruction process is highly related with multiple scattering phenomena, which elevates both the complexity and computational costs for such efforts. This work introduces a time-based Markov chain method, which uses a sparse transition matrix to represent the likelihood for a photon to transit in the turbid media. The accuracy of the time-based Markov chain model was veriﬁed against the forwarding calculations of the scattering-based Markov chain model and Monte Carlo simulations. Then, reconstruction was performed with backscattered photon angular distributions. Based on the characteristics of the sparse transition matrix, the optical properties (droplet diameters) could be obtained layer by layer with transmitted photon distributions at different time durations. It is shown that the time-based Markov chain model can reconstruct the optical properties of a turbid slab with satisfactory accuracy and lower computational costs.


Introduction
There is an increasing need for interpreting light scattering signals from turbid media, since scattering signals can contain crucial information about the optical characteristics of the media. The scattering process is usually complicated due to the occurrence of multiple scattering. Literally, multiple scattering means a photon or ray of light is scattered more than once within the medium before transmitting (i.e., leaving the scattering medium). Such chaotic motions can significantly impact optical diagnosis in, for instance, optical computed tomography techniques. Meanwhile, multiple scattering signals can be used to infer scattering phase function-related information by tracking the transmitting process of the photons. Therefore, it is expected that multiple scattering signals can be used for applications such as tissue diagnosis [1,2], atmosphere detection [3] and spray droplet size estimation [4], among others.
However, the inversion procedure for reconstruction of the turbid media can be complex when considering multiple scattering. When the optical depth (OD) of the medium (representing an expected number of scattering events for a single photon penetrating the medium) is within the range of 2-10, no satisfying approximation can be applied to obtain the solution of the turbid media for the scattering process. The most commonly seen methods to solve the multiple scattering problem are solving the radiative transfer equation (RTE) [5][6][7] or using Monte Carlo simulations [4,[8][9][10]. As an analytical method, RTE can be solved to obtain the exact solutions, while it is challenging and time-consuming unless strong assumptions are made [11][12][13][14]. The Monte Carlo simulation method is also popular in that the algorithm itself is relatively simple. The scattering event in a Monte Carlo simulation is modeled using a random number generator, which determines whether or not a photon will be scattered at a certain location and the new propagation direction. The disadvantage of Monte Carlo simulations is also notable. The Monte Carlo simulation is a stochastic algorithm in which a great number of photons need to be sent to warrant statistically meaningful results, especially when the signal level is expected to be low. Furthermore, the stochastic nature of this algorithm prevents any analytical inversion or reconstruction approach based on Monte Carlo simulations.
In previous works, the authors proposed a Markov chain-based model [15][16][17][18] to calculate the transmitted angular distribution through turbid slabs. In this model, multiple scattering was converted into a matrix (transition matrix), which represented the likelihood of a photon to transit from one state (location and propagation angle) to another. It has been demonstrated that this method can obtain results showing angular distribution (i.e., the probability distribution of transmitted photons as a function of transmitted angles) through infinitely large turbid slabs with the assumption that the optical properties within each parallel layer is uniform. It has been demonstrated that the Markov chain method could incorporate practical scattering phase functions and be executed with a much lower computational cost. In the recent work of the authors, the angular distribution results from photons scattered once and twice were used as inputs, and the optical properties of the turbid slab (droplet diameter) were reconstructed with good accuracy and noise resistance. The use of global optimization algorithms notably increases the computational costs to around 150 minutes, although such a capacity is new and has hardly been demonstrated before.
However, the previously introduced Markov chain method is still not perfect, since it relies on obtaining the angular distribution of the photons with an exact scattering order. Although possible, such a requirement is hard to achieve in practical experiments. For practical optical diagnostics, it is more feasible to obtain temporally resolved data, such as the spatial intensity distribution of transmitted photons or the angular distribution data as a function of time. These types of measurements have been proven to be practical in the existing literature [19]. In the previous Markov model, the scattering was set as the event which the algorithm would track. As a result, the scattering number (scattering order) of a photon will be preserved by the algorithm, but the total path length or the residual time cannot be kept in the current Markov model. In order to find temporal resolved angular distribution using the Markov chain model, the time or path length must be set as the event so that this information can be preserved as the photons are transmitted.
In this work, the authors propose a so-called time-based Markov chain model, which aims to address the issues with the scattering-based model studied before. Because of the nature of the newly introduced time-based Markov chain model, the inversion or reconstruction can be done without using global optimization algorithms. Rather, using the backscattered angular distribution, the optical properties of the turbid slab can be inferred one-by-one using linear algebra operations, which benefits from the matrix form of the time-based Markov chain model. The rest of this work is organized as follows. Section 2 discusses the fundamentals and assumptions of the time-based model. Section 3 demonstrates the numerical verification and reconstruction results using the time-based model. Finally, Section 4 concludes this investigation.

Mathematical Formation of the Time-Based Markov Model
The mathematical grounds of both the scattering-based Markov chain model and the time-based Markov chain model are similar. In the model, a transition matrix is built, and photon transmission results are given by matrix manipulations, such as matrix multiplication and inversion. The major difference between the two models is the definition of what an event is in the Markov formation. Such a difference in definition is depicted in Figure 1. Figure 1a demonstrates the light propagation procedure of the scattering-based model, and Figure 1b shows the same procedure for the time-based model, respectively. In scattering-based models, each photon scattering is considered an event. Therefore, the photon shown in the figure experienced two events before being transmitted. In this formation, the corresponding transition matrix is usually nonzero, since there is always a nonzero probability for a photon to propagate from layer i to layer j without scattering during the propagation. In comparison, the time-based model considers that a new event takes place once a certain amount of time has elapsed. For instance, if we define the time during which the photon propagates for one layer, then the photon will experience approximately n events before transmission if there are n layers in total. It is also worth noting that the propagation angle will affect the time for a photon to pass through one layer; we just used this illustration for demonstration purposes. In this model, the transition matrix can be established by mostly zero entries. For instance, if the set time interval t corresponds to the time that a photon propagates for a distance of the thickness of one layer, the photon cannot propagate for more than one layer during this time period, which means only those representing probabilities from the first layer to the second layer are not zero in the transition matrix. With such characteristics, it is possible to reconstruct the optical properties of the turbid slab using linear algebraic methods.
Although the time-based model is intrinsically similar to the scattering-based model, two assumptions are made to simplify the establishment of the time-based model. The first assumption is that the photon will be scattered at most once during one time period, and the scattering happens at the end of this time period. The error induced by this assumption can be reduced by finer discretization so that the likelihood of multiple scattering within one time period can be minimized. Assuming the OD in a layer is OD ∆z , the probability of multiple scattering can be readily found by Equation (1): In this expression, the second term on the right-hand side represents the probability that the photon transmitted a specific layer with an OD (with the value of OD ∆z ) without being scattered, and the third term represents the likelihood that this photon is scattered only once in this layer. It can be readily found that the error introduced by this assumption is 1.97%, with an OD ∆z = 0.2, and 9%, with an OD ∆z = 0.5. Since 100 layers are used in this current model, this error is acceptable, but as the OD in each layer grows, this error also increases and accumulates with each scattering, which should be taken into account in actual applications. This error can be reduced by using finer spatial meshes or taking a higher scattering order into consideration, which will make the algorithm more complex.
The second assumption is that in each layer, the photon starts its propagation at a random location z within this layer rather than a preset location. This assumption is also utilized in the scattering-based Markov chain model. However, the starting location for the scattering-based Markov chain model is always where scattering happens, which is a random location in the layer. In comparison, for the time-based model, if the photon is not scattered in a specific event, then the ending location of this photon is actually fixed rather than randomized. This error can be reduced with finer spatial meshes.
The governing equations for the time-based model are also similar to those under the scattering-based model, and the following equations are still applicable: In the equations, P(z m , z n , θ i ) is the probability for the photon that starts its propagation at z m with the propagating direction of θ i and finishes at z n during a specific time step; P(θ i , θ j , n) is the probability for the photon that propagates in the direction of θ i , snth layer, and then propagates in the new direction of θ j ; and Γ n (α) is the phase function in the nth layer. The implications of Equations (2)-(4) have been introduced elsewhere [15][16][17][18], and we will briefly introduce these equations. Equation (2) describes the transition probability T from one state (z m , θ i ) to another state (z n , θ j ). This probability is co-determined by the probability of the path length the photon propagates before being scattered again and the probability of the new scattering angle θ j after being scattered. The relationship between θ i and θ j can be established by the scattering phase function Γ n (α) with α, which is the included angle between the incident light and the scattered light. The phase function is determined by the imaging system characteristics. ϕ is the azimuth angle, which impacts the determination of the included angle α. To consider all ϕ in determining P(θ i , θ j , n), an integral is solved, as shown in Equation (4) over 2π.
In the time-based model, the expression of P(z m , z n , θ i ) is as follows: where l is the path length for each Markov event, t = l/c and c is the speed of light in the scattering medium. Equation (5) is readily understood with the assumptions stated before. As can be seen, when θ i is fixed, the direction of the photon is fixed (either forward or backward). Then, the photon can only propagate to either the current layer or the next layer (either backward or forward), and the transition probability can be easily found with Equation (5). Therefore, only z m and z n have nonzero probabilities, which means the transition matrix T under the time-based scheme is a sparse matrix. With such definitions, the transition matrix T can be obtained with each entry defined by Equations (2)-(5), and we have: and: where Q n is the angular distribution (the weight of the photons in each angle) of all photons that have transmitted during the time (n − 1) × t and n × t, Q total is the angular distribution of all transmitted photons, P is the initial photon status and R is the absorption matrix. With the structure above, the time-based Markov chain method can simulate the timeresolved angular distributions for transmitted photons, which is more accessible than counting the scattering order of the photons practically. More mathematical details can be found in [15][16][17][18].

Numerical Verification against Monte Carlo Simulations
In this section, we compare the performance of the time-based Markov chain model against the Monte Carlo simulation model. Although it would be more desirable to directly compare the results from the time-based model to the experiments, matching all the parameters in such experimental efforts can be challenging. Therefore, we used the Monte Carlo algorithm in [4] as the benchmark in this work, which was validated against a few experiments under multiple scattering conditions. In the Monte Carlo model, the trajectory and scattering events of each photon were determined by randomly generated numbers and then tracked and monitored by the algorithm. Under each condition, 100 million photons were sent. Results from the scattering-based model will also be presented, but the scattering-based Markov model cannot estimate temporally resolved angular distributions for transmitted photons, so we only compared the results of the time-based Markov chain model and the Monte Carlo model when applicable. Figure 2 compares the total angular distributions under different ODs when using different algorithms. The phase functions used were in agreement with those in our previous work [18]. The external light wavelength used in this simulation was set 532 nm, which is usually seen from Nd:YAG laser sources for diagnostics purposes. The scattering medium was water, and the ambient medium was air. The water droplet diameter ranged from 1.0 µm to 2.0 µm (a diameter resolution of 0.1µm was used). The Mie scattering theorem was adopted to calculate multiple scattering characteristics. In Figure 2, and in the rest of the simulation, a uniform OD distribution was used, and the phase function from the 1.0 µm water droplets was adopted. The turbid slab was assumed to be infinitely large, and the optical properties of the slab were uniform across the plane perpendicular to the thickness direction (z). All photons entered the slab from the same point and initially propagated in the z-direction. The turbid slab was discretized into 100 layers in the z-direction, and the propagation angle θ was discretized into 180 bins. Therefore, the dimensions of the transition matrix were 18,000 × 18,000, which was the same as in the previous work of the authors. For more detailed interpretation of the results shown in Figure 2, each data point represents the probability that the photon is transmitted in the propagation angle θ. For instance, the value at θ = 45 • reflects the percentage of photons transmitted with 45 • < θ < 46 • . Ballistic photons were excluded (photons without scattering) in the figure, and the total percentage plotted from θ = 0 • to θ = 180 • (including ballistic photons) equaled 1.0. In the Monte Carlo simulations, such probabilities could be obtained by counting individual photons statistically and then dividing them by the total number of photons simulated by the Monte Carlo model.
Although it was not the purpose of this work to compare the computational cost and accuracy, since the time-based model was mainly designed for reconstruction purposes, it is worth noting that the computational cost of the time-based model was around 5 min, and it took about 60 min for the Monte Carlo algorithm to simulate 100 million photons under OD = 1.0 conditions. Even for different methods which have been proposed to reduce the computational cost of Monte Carlo simulations, such as using cloud computing [8,9] or cell processors [20], hardware improvement could also accelerate the computation of the Markov chain model, especially when no computational optimization was done to the Markov chain time-based model at the time. In this work, we just demonstrated the capability of the proposed Markov time-based model and its inversion method. We will focus on the computational cost and accuracy in our future work.
As can be seen in Figure 2, the newly developed time-based model approximated the results given by the scattering-based model and Monte Carlo simulations well, although the resultant angular distribution could be complex and wavy. Furthermore, as we have discussed, the accuracy of the time-based algorithm suffered from a higher OD, and it is notable that the results under OD = 10 and OD = 20 cases were a bit off from the ground truth when compared with the OD = 2 and OD = 5 cases. Nevertheless, the proposed time-based model could still achieve a relative error of less than 1% compared with the Monte Carlo simulation results, which indicates that the proposed time-based model can satisfy the need of approximating Monte Carlo simulations.
To further demonstrate the error from the Markov chain approximation, we incorporated the relative error of each model, as shown in Figure 3. As we introduced in our previous work [17], the Markov chain model had the worst performance near 90 • , around which the signal was hard to capture by the sensors in the experiments. Therefore, in Figure 3, we just compared the results at 135 • (corresponding to the maximal point in 90 • -180 • ) to demonstrate the error of the time-based model, and we also confirmed that this error was representative of the overall error of the Markov chain models. As can be seen, the relative error of the scattering-based model increased as the OD increased, but in general, the relative error was less than 0.5%, which was in agreement with our previous research [16]. For the time-based model, the error was around 0.7% when OD = 2. Then, the error increased steadily and reached~1.6% when OD = 10. The error increased less remarkably afterward with higher ODs. Thus, although the time-based model would incur error from the two assumptions mentioned above, the overall accuracy was still acceptable for calculating the transmitted angular distribution. In order to examine the capacity of the time-based Markov chain model for distinguishing transmitted photons at different timings (i.e., propagating different path lengths inside the turbid media), Figure 4 compares the results obtained by Markov chain approximation and the Monte Carlo simulations, respectively. In Markov chain realization, the photon transmitted can be obtained via Equation (6). For instance, the transmitted angular distribution within t = 100∆z/c-t = 125∆z/c can be achieved by adding Equation (6), with n ranging from 100 to 125. For Monte Carlo simulation, as each photon is transmitted, the pathlength that this photon has propagated for can be recorded. Therefore, when all the photons have been simulated, the sought distribution within a range from t = 100∆z/c to t = 125∆z/c can be obtained by counting all photons that propagate a total path length of 100-125∆z before being transmitted. In this investigation, a total of 100 turbid layers were used, and the results in Figure 4 were obtained under a total OD of 5.0, suggesting that the OD in each layer was 0.05. The thickness of each layer was noted as ∆z. Figure 4a-h presents the photons transmitted during different periods of time. Since the total thickness of the turbid slab was 100∆z, no photon could be transmitted from the exit plane before t = 100∆z/c because 100∆z was the shortest path length possible. Therefore, in Figure 4a-d, no photons were seen within the range from 0 • to 90 • , and the transmitted photons in these panels were all photons reflected back to the entrance plane. Figure 4e-h demonstrates the transmitted photons during t = 100∆z/c and t = 200∆z/c (with ballistic photons removed from the results). As can be seen, after t = 100∆z/c, the fluctuating features associated with the Mie scattering phase function were not pronounced for the backscattered photons due to diffusion effects. Finally, the time-based Markov chain model could generally approximate the angular distribution of the transmitted photons that exited the turbid slab at different times, which could not be achieved by our previous scattering-based Markov chain model. Thus, the performance of the time-based Markov model is verified. The proposed method can approximate the angular distributions received by the sensors at different moments, which are easier to examine experimentally and provide more information for the reconstruction purpose explored in the following section.

Time-Based Markov Chain Multiple Scattering Inversion
The time-based model can utilize temporal resolved transmitted photons as reconstruction inputs, which is more desirable compared with finding the scattering order (i.e., how many times a specific photon has been scattered before transmission) of the transmitted photons experimentally. Furthermore, for the backscattered photons, there is a maximum penetrating distance during a specific time period. For instance, within t = 0~20∆z/c, the furthest distance a back-transmitted photon can reach is 10∆z, which requires this photon to propagate perpendicular to the entrance plane (θ = 0 • ), scattered at z = 10∆z, then return to the entrance plane with a transmitted angle of θ = 180 • . Therefore, the optical properties within z = 0~10∆z can be inferred by the transmitted angular distribution obtained during t = 0~20∆z/c. In a similar approach, the optical properties within z = 10~20∆z can be found by the transmitted angular distribution obtained during t = 20~40∆z/c, together with the known properties in z = 0~10∆z. In an iterative manner, the spatially resolved optical properties of the turbid slab can be reconstructed without brute force trial and error of all possible combinations.
The schematic of such a method has been depicted in Figure 5. The transition matrix to be solved was initially a zero matrix. With the angular distribution data during the first time period (Iter 1), the optical properties of the first set of layers (10 in this work) could be found by minimizing the difference of the simulated angular distribution and the measured angular distribution, assuming the optical properties in the first set of layers were uniform. The entries at other locations (for instance, in Iter 2) had no contribution to the measurements, so the entries in these locations could be set to 0. After fixing the optical properties in the first set of layers (in this work, the droplet diameter), the corresponding entry values (denoted by x in Figure 5) were updated and fixed in the iterations that followed. Next, for Iter 2, we performed the same procedures and updated the corresponding values (denoted by y), and so on so forth. Eventually, we solved for the final set of layers and updated the rest of the entries (denoted by z). By using this scheme and backscattered photons (backscattered photons were more readily accessible in the experiments), we successfully decoupled the connections between different layers through multiple scattering and substantially reduced the computational cost compared with using Monte Carlo simulation or brute force calculations. To examine such a capacity of the time-based model, we established a synthetic slab with a customized water diameter distribution for reconstruction. The numerical setup of this reconstruction effort is depicted in Figure 6. The synthetic phantom was in agreement with the one in our previous work [18] (i.e., the total OD of the turbid slab was 5.0). The slab was discretized into 10 layers with the same OD ∆z and OD. The diameter of the water droplets in layers 1-3 and 8-10 was 1.3 µm, and for layers 4-7 it was 1.6 µm. In this simulation, we had 11 candidate droplet diameter sizes (1.0-2.0 µm, with a step length of 0.1 µm). Therefore, there were 1110 different possible diameter spatial distributions. Solving for the true diameter distributions by brute force could be impossible in the sense of computational cost, even with the Markov chain models. However, the proposed reconstruction method reduced the cost function evaluations to only 11 (number of candidates) × 10 (number of layers) = 110, which substantially reduced the computational cost and made the reconstruction feasible. The reconstruction performance of the time-based Markov inversion method is presented in Figure 7. The ground truth data was obtained using the forwarding time-based Markov model, and we also added 5% relative noise to the ground truth data to investigate the robustness of the algorithm under noise. Since the time-based inversion scheme is a deterministic method, only one reconstruction was performed, in comparison with the 10 tests using optimization algorithms in the previous work [18]. The computational cost to perform each evaluation was around 25 minutes with an Intel i7-8700 CPU using MATLAB software, which could potentially be improved by using other programming platforms. In the two evaluations performed, both reconstructions faithfully rebuilt the synthetic phantom, so we will not be showing the reconstruction here. Instead, we show the relative error of each diameter in each iteration in Figure 7. In each iteration, the transmitted angular distribution was obtained with an assumed diameter size. Then, the relative error for each angle was found against the ground truth value. In this reconstruction, only backscattered photons were used, and we chose a range within 100-170 • , where the signal strength was higher. Then, we summarized all the absolute values of the relative error and found the droplet diameter with the lowest relative error summation. As seen, the cost function values in the first few layers were distinctively lower at 1.3 µm (~0), compared with the rest of the candidates (~70). As the sought layer moves further (away from the entrance plane), the relative error using false diameter candidates also decreases and reduces the sensitivity of the reconstruction. It is worth noting that this issue is not only associated with the proposed Markov chain method, but also impacts all algorithms utilizing backscattered photon angular distributions. The reason for this is that the sensitivities of different phase functions will be reduced, as the photons have been scattered many times, making it so that the spectral features of the angular distribution are not apparent anymore, which is also evident in Figure 5e-h in the range from 90-180 • , where the angular distribution was a quite smooth single-peak distribution. Furthermore, although the true particle diameter has been found with a 5% relative error, as shown in Figure 6b, the summarized error in layers 81-90 for the true diameter (1.3 µm) was 0.82, compared with 0.89 for a diameter of 1.4 µm. Such an observation means the capability to distinguish different phase functions is deteriorated with a higher measurement error. Therefore, an OD threshold exists beyond which the droplet diameter cannot be resolved by backscattered photons, and the time-based Markov method is also limited by such a constraint.
To address the issues mentioned above, photons transmitted from the exit plane (forward scattering photons) can also be utilized to examine the fidelity of the reconstruction with additional computational costs. Such efforts will be examined and implemented in future work. Furthermore, with the forward scattering process converted into a fast matrix evaluation problem, it is also possible to utilize end-to-end deep learning algorithms for reconstruction purposes, for which the training dataset will be provided by the time-based Markov model. Such a scheme is believed to be able to further enhance the current capacity of the proposed time-based Markov model.