Markov Chain Investigation of Discretization Schemes and Computational Cost Reduction in Modeling Photon Multiple Scattering

Establishing fast and reversible photon multiple scattering algorithms remains a modeling challenge for optical diagnostics and noise reduction purposes, especially when the scattering happens within the intermediate scattering regime. Previous work has proposed and verified a Markov chain approach for modeling photon multiple scattering phenomena through turbid slabs. The fidelity of the Markov chain method has been verified through detailed comparison with Monte Carlo models. However, further improvement to the Markov chain method is still required to improve its performance in studying multiple scattering. The present research discussed the efficacy of non-uniform discretization schemes and analyzed errors induced by different schemes. The current work also proposed an iterative approach as an alternative to directly carrying out matrix inversion manipulations, which would significantly reduce the computational costs. The benefits of utilizing non-uniform discretization schemes and the iterative approach were confirmed and verified by comparing the results to a Monte Carlo simulation.


Introduction
Multiple scattering of photons propagating through turbid (optical opaque) medium has remained a challenge for many multi-disciplinary applications, such as biomedical tissue diagnostics and aerosol characterization, smoke and fog analysis, and spray.The reason that simplification assumptions usually do not apply well for such problems is that the multiple scattering phenomenon is intrinsically stochastic, particularly in the intermediate regime where the optical depth (OD) of the media is somewhere between two and 10.Many have attempted to interpret multiple scattering through modeling using Monte Carlo methods [1][2][3][4], or other approaches grounded on the random walk theory or empirical approximations [5][6][7][8].
Monte Carlo methods have been widely applied in multiple scattering numerical modeling because of the methods simple implementations and statistical outcomes.Monte Carlo algorithms provide statistical predictions for stochastic processes through repeated random sampling within the domain of interest.More specifically, for scattering simulations, the model continuously sends one photon into the scattering medium and tracks its trajectory by predicting its propagation and scattering events.If the photon experiences a scattering event, the new propagation angle will be calculated based on the scattering phase function, determined by the physics of the scattering process.This process ends up with the photons leaving the domain of interest.
Albeit that the Monte Carlo simulation can be readily implemented and allows simulations with a high fidelity, it suffers from a high computational cost for practical applications.Such a computational cost is primarily due to two reasons: (1) Monte Carlo algorithms require a great number of samples to guarantee the desired accuracy.Simulation error from Monte Carlo is inversely proportional to the square root of the number of samples N. Therefore, a sufficiently large number of samples, usually millions or billions of photons in this case, is required by the Monte Carlo model for multiple scattering problems to achieve meaningful predictions; (2) The total computational cost of the Monte Carlo simulation is also impacted by the cost of each evaluation, i.e., starting from when one photon enters the domain of interest and ending when the photon leaves the domain.The computational cost of the Monte Carlo algorithm is significantly higher with a larger OD because more scattering events are expected and more evaluations are required.
Given a large sample size and the relatively high computational cost in evaluating each photon, Monte Carlo simulations used in multiple scattering usually take hours for one single run with a single core CPU.The computational cost is less for forward problems as long as the optical properties of the turbid medium are known, but it is prohibitively high for inverse problems where the optical properties of the turbid medium are sought.Thousands or millions of cost function evaluations (multiple scattering calculations) are required for global optimization/reconstruction algorithms, such as simulated annealing (SA) algorithms and genetic algorithms (GA) in practical applications [9][10][11][12].Consequently, alternative algorithms and models, such as the matrix inversion methods [13] or iterative methods [14] are needed to enable multiple scattering diagnostics and reconstructions with more reasonable computational costs.With such computationally efficient methods and algorithms, more applications will be made possible, such as cancer diagnosis for skin tissue.Meanwhile, the capacity of many existing technologies will be further expanded, such as landscape imaging.
Over the past decades, the Markov process has been utilized for multi-disciplinary, probability-based investigations and applications, such as economic predictions [15], computer vision [16], production control [17], and gambling [18].There are also Markov chain related studies on the topics of diffusion [19] and radiative transfer [20].In previous work [21,22], the authors have introduced a Markov chain approach for modeling transmitted photon angular distribution (i.e., the probability density function of propagation angles of transmitted photons) via multiple scattering based on the Markov chain approximation.There has been some research using the Markov chain method to model the multiple scattering problems.For instance, Geddes et al. [23] have used the Markov chain radiative transfer to invert OII 83.4 nm dayglow profiles with multiple scattering considerations.Yuan et al. [24] utilized the Markov chain approach to model obliquely incident rays on the surface of participating media for atmospheric research.Although Markov chain methods in dealing with light propagations in multiple scattering medium can be directly derived by discretization of the generalized radiative transfer equations (RTE) [25], the modeling process can be independent from solving the ordinary differential equations (ODE) in typical RTE solution procedures [26].Such a feature would enable customized Markov chain approximations based on the needs of the investigation without solving the whole RTE system, and the authors' previous work using the Markov Chain has successfully built a precision connection between actual measureable data (angular distribution) and the desired optical properties of the medium, with less computational costs.
The discrete Markov chain method discretizes the domain of interests into multiple "states", and the transition probability from one state to another can be found for specific problems.A transition matrix T is then built with each entry T ij representing the probability of transition from state i to state j, so that the Markov chain approach converts the stochastic, random process problem into a mathematical problem involved with matrix manipulations.The Markov chain approximation approach for photon multiple scattering has been verified for both isotropic scattering and anisotropic scattering through the turbid slab and more details regarding Markov chain approximations for multiple scattering can be found in the authors' previous publications.
Despite the many advantages of Markov chain approximations over Monte Carlo simulations, the method suffers from a variety of computational complexities.For example, the "curse of dimensionality" [27].Here, the exponential increase in computational cost with higher dimensions occurs due to the size increase of the transition matrix (and the number of states).Techniques have been developed to improve performance at higher dimensions.Illustrations of these techniques include state reduction methods, such as the lump process established by Kemeny and Snell [28] and the GTH elimination algorithm (proposed by Grassmann, Taksar, and Heyman) [29].Methods for improving the performance of Markov chain approximations also include investigations into proper discretization of continuous Markov chain problems [30,31].
For multiple scattering problems, it is necessary to properly discretize the domain of interest to improve the efficiency of Markov chain methods, and thus reduce the computational cost.It is also of interest to find an effective way to discretize the zenith angle (photon propagation angle), so that the transmitted and reflected photon angular distributions can be simulated with improved fidelity.Such investigation is expected to improve the performance of the Markov chain method, for that many scattering applications are in Mie scattering regimes with high anisotropies (g), indicating that highly forward scatterings are desirable.Finer discretization in the forward direction could potentially better resolve the resultant angular distribution.Furthermore, there is usually a detection limit of the transmitted/reflected photon angles experimentally, for example, up to 580 m/rad (33.2 • ) [32,33].This issue also suggests an angular discretization that emphasizes finer resolution within the detection limit.
In the authors' previous works on Markov chain modeling for multiple scattering, uniform discretization for both the OD and zenith angle was used.Alternative discretization methods can be used to further improve the accuracy of the proposed Monte Carlo model with reduced or similar computation costs.As such, the goal of this paper was to investigate the impact of discretization on the Markov chain approximation of anisotropic multiple scattering through a turbid slab, and to recommend appropriate discretization schemes for applications.Furthermore, it is desirable to optimize the mathematical realization of the Markov chain model in the aspect of computational cost and accuracy.As such, methods to accelerate Markov chain approximations will also be discussed in this work.

Theory
The Markov chain model for multiple scattering in mediums has been thoroughly introduced in Li et al. [21,22], assuming that the medium consists of parallel layers (the optical properties within a specific layer is uniform) and the light propagation follows the Beer-Lambert law, where the extinction coefficient consists of absorption and scattering.For simplicity, only the Mie scattering was considered in this research.
The fidelity of the Markov Chain scheme in dealing with inhomogeneous medium, absorption, and anisotropic phase functions have been demonstrated.Therefore, the authors will only briefly introduce Markov chain theory and its application to photon multiple scattering.Basically, the Markov chain model for anisotropic scattering defines a transition state by the location z in the 1D turbid slab, where a scattering event takes place, and the direction/angle θ (zenith angle) in which this photon will propagate after the scattering event.Essentially, the multiple scattering problem can be considered as a continuous Markov chain problem.To solve the Markov process, the scattering medium is discretized to convert this problem from a continuous challenge into a finite Markov chain problem, so that matrix forms can be implemented.With the discretization approach, a transition matrix T can be defined with each entry T((z m , θ i ), (z n , θ j )), representing the probability of the photon transit from one state (z m , θ i ) to another (z n , θ j ).Absorbing states are incorporated to represent the probability of photons exiting the slab or being absorbed by the medium.More details on the Markov formation can be found in Li et al. [22].From photon scattering physics, the following relationships can be inferred for the multiple scattering problem, as in Li et al. [22]: α = arccos(cos θ i cos θ j cos ϕ + sin θ i sin θ j ) Whereas, in Equation ( 1), P(z m , z n , θ i ) is the probability of a photon with two consequent scattering events taking place at z m , z n and propagates in the direction of θ i , z m represents the location of the prior scattering event and z n represents the location of the latter scattering event.P(θ i , θ j , n) is the probability that a photon propagates in the direction of θ i , scattered with a certain phase function in the nth layer, then propagates in the new direction of θ j .In Equation (2), Σ e,z is the extinction coefficient at the location z.In Equation (3), ϕ is the azimuth angle and α is the included angle between the incident direction and the exiting direction.In Equation ( 4), Γ n (α) is the phase function (probability density function, PDF) in the nth layer.It is worth noting that Equation (2) represents the propagating probability from point z m to point z n .To solve the probability of the photon propagating from layer m to layer n, proper discretizations and integrations are necessary.To derive such a probability, the cumulative probability density function (CPDF) of the OD can be defined as C OD (z).The upper lower and upper bound of layer m and n are defined to be z ml , z mu , and z nl , z nu , respectively.With simple mathematical manipulations, the discretized form of Equation ( 2) can be readily established after integrating the probability over the thickness, twice for each layer: Furthermore, Equation (4) requires numerical integration instead of analytical solutions.To the authors' knowledge, a general analytical solution for any given phase function of Equation ( 4) is not readily available for analytical investigation purposes.Therefore, in the Markov chain calculations, ϕ is discretized (i.e., into 900 equal intervals at a resolution of 0.2 • ) to numerically solve for P (θ i , θ j , n).
In this study, various discretization schemes are utilized to analyze the effect of discretization.The schematic shown in Figure 1 depicts the comparison of uniform discretization (as in [21,22]) and non-uniform discretization which are proposed by this work.The upper half of Panel (a) represents the uniform discretization and the lower half shows the non-uniform discretization, and in both schemes, 180 intervals were created.For uniform discretizations, every internal is equally divided by a one-degree resolution.For non-uniform discretizations, different discretization schemes are tested, for example, finer meshes/intervals adopted for forward directions (α~0 • ) and coarse meshes/intervals for backward directions (α~180 • ).Panel (a) of Figure 1 is a polar plot demonstration of a Mie scattering phase function with a resolution of 1 degree.In this work, this specific phase function shown in Panel (a) of Figure 1 is defined as phase function 1 in this manuscript.For phase function 1, a uniform particle diameter of 1.2 µm, incident light wavelength of 532 nm, and ambient refractive index of 1.0 were selected.The refractive index of the particle was set to be 1.33 and unpolarized light was assumed.Different thickness discretization schemes are shown in Panel (b) of Figure 1.The left-hand side of Panel (b) depicts a uniform discretization in which each layer has the same thickness ∆z.Correspondingly, the right-hand side of Panel (b) shows a non-uniform discretization for which the thickness ∆z in each layer may vary.The authors believe that different emphasis should be given for different cases.For example, finer meshes near the incident plane of the turbid are desired if the reflected angular distribution is sought.However, the relevance of such discretization is subject to debate for transmitted photon angular distributions, which the authors will also investigate in this work.Two additional phase functions were considered in this work and the three phase functions are depicted in Figure 2 in Cartesian coordinates.Phase function 2 is defined based on a particle refractive index of 1.57, and everything else is the same as phase function 1.To create a more structurally complex phase function, phase function 3 is defined based on a uniform particle diameter of 10.0 μm, incident light wavelength of 650 nm, a particle refractive index of 1.57, and an ambient refractive index of 1.0.The phase functions are depicted in Figure 2 and the anisotropy of phase function 1, 2, and 3 are determined to be 0.8314, 0.6653, and 0.8691, respectively.For Markov chain approximations, the transition matrix calculations and the matrix inversions were carried out with MATLAB because of its excellent performance in handling matrix manipulations.An Intel Xeon X5650 CPU (12 threads, 2.66 GHz) was unitized for the MATLAB calculations.It is worth noting that MATLAB uses parallel computing for matrix calculations; therefore, when evaluating the computational costs, the multi-core parallelization effect should be considered as well.
To evaluate the accuracy of different discretization schemes, the Monte Carlo method was chosen as a reference, since it is a widely-accepted standard for multiple scattering simulations.The fidelity of the Monte Carlo codes utilized in this research have been compared and verified with other methods [5,34], which can be found in Li et al. [21,22].To compare the Markov chain simulations with the Monte Carlo simulation results, a fine angular discretization of the phase function was used in the Monte Caro simulation, i.e., 0.05°.In this study, 2.5 billion photons were sent for each Monte Carlo simulation.The program was written in C language and run with the same type of CPU (single Two additional phase functions were considered in this work and the three phase functions are depicted in Figure 2 in Cartesian coordinates.Phase function 2 is defined based on a particle refractive index of 1.57, and everything else is the same as phase function 1.To create a more structurally complex phase function, phase function 3 is defined based on a uniform particle diameter of 10.0 µm, incident light wavelength of 650 nm, a particle refractive index of 1.57, and an ambient refractive index of 1.0.The phase functions are depicted in Figure 2 and the anisotropy of phase function 1, 2, and 3 are determined to be 0.8314, 0.6653, and 0.8691, respectively.Two additional phase functions were considered in this work and the three phase functions are depicted in Figure 2 in Cartesian coordinates.Phase function 2 is defined based on a particle refractive index of 1.57, and everything else is the same as phase function 1.To create a more structurally complex phase function, phase function 3 is defined based on a uniform particle diameter of 10.0 μm, incident light wavelength of 650 nm, a particle refractive index of 1.57, and an ambient refractive index of 1.0.The phase functions are depicted in Figure 2 and the anisotropy of phase function 1, 2, and 3 are determined to be 0.8314, 0.6653, and 0.8691, respectively.For Markov chain approximations, the transition matrix calculations and the matrix inversions were carried out with MATLAB because of its excellent performance in handling matrix manipulations.An Intel Xeon X5650 CPU (12 threads, 2.66 GHz) was unitized for the MATLAB calculations.It is worth noting that MATLAB uses parallel computing for matrix calculations; therefore, when evaluating the computational costs, the multi-core parallelization effect should be considered as well.
To evaluate the accuracy of different discretization schemes, the Monte Carlo method was chosen as a reference, since it is a widely-accepted standard for multiple scattering simulations.The fidelity of the Monte Carlo codes utilized in this research have been compared and verified with other methods [5,34], which can be found in Li et al. [21,22].To compare the Markov chain simulations with the Monte Carlo simulation results, a fine angular discretization of the phase function was used in the Monte Caro simulation, i.e., 0.05°.In this study, 2.5 billion photons were sent for each Monte Carlo simulation.The program was written in C language and run with the same type of CPU (single For Markov chain approximations, the transition matrix calculations and the matrix inversions were carried out with MATLAB because of its excellent performance in handling matrix manipulations.An Intel Xeon X5650 CPU (12 threads, 2.66 GHz) was unitized for the MATLAB calculations.It is worth noting that MATLAB uses parallel computing for matrix calculations; therefore, when evaluating the computational costs, the multi-core parallelization effect should be considered as well.
To evaluate the accuracy of different discretization schemes, the Monte Carlo method was chosen as a reference, since it is a widely-accepted standard for multiple scattering simulations.The fidelity of the Monte Carlo codes utilized in this research have been compared and verified with other methods [5,34], which can be found in Li et al. [21,22].To compare the Markov chain simulations with the Monte Carlo simulation results, a fine angular discretization of the phase function was used in the Monte Caro simulation, i.e., 0.05 • .In this study, 2.5 billion photons were sent for each Monte Carlo simulation.The program was written in C language and run with the same type of CPU (single thread for each Monte Carlo simulation).The results obtained from the Monte Carlo simulation were then used as benchmarks.

Uniform Discretizations
In this section, Markov chain solutions with uniform discretizations were solved and compared with Monte Carlo simulations.It is worth noting that in the zenith angle discretizations (to be discussed in Section 3.3), both the incident angle θ 1 and new propagation angle θ 2 need to be discretized and the same discretization scheme must be applied to both angles.Such a requirement is posed by the nature of the Markov chain approximation, in that the transition states cannot be re-defined during the Markov process.Similarly, the resulting transmitted and reflected angular distributions tend to follow the same discretization that one sets at the beginning.Therefore, direct comparison among different discretization schemes will not offer sizably meaningful insights.Instead, the Markov calculations will be compared against one common benchmark with the results from the Monte Carlo simulations.Transmitted probabilities were normalized with a unit of 1 •−1 in this work to make sure that the summation of the probability over the range from 0 • to 180 • was unity minus the percentage of ballistic photons (transmitted photons without scattering, not discussed in this work).The same method was applied for all results presented in this work.
Figure 3 depicts Markov chain predictions of the transmitted photon angular distributions compared with those from the Monte Carlo simulations.In the simulations, uniform OD distribution was assumed, and the total OD was set to either 2 or 5. On the x axis, the transmitted angle from 0 • to 90 • indicates that the photons leave the turbid slab from the bottom exit plane, whilst the transmitted angle from 90 • to 180 • represents the photons that are reflected and leave the slab from the entrance plane.Since the range of the OD in the simulations was within the intermediate regime, the transmitted angular distribution exhibited both features from the phase function and multiple scattering averaging effect, as reflected by the complex angular distribution structures, as reflected in Figure 3. From the results, it is evident that Markov chain approximations agree well with Monte Carlo simulation results with satisfactory fidelity, for all the six cases considered in this work.thread for each Monte Carlo simulation).The results obtained from the Monte Carlo simulation were then used as benchmarks.

Uniform Discretizations
In this section, Markov chain solutions with uniform discretizations were solved and compared with Monte Carlo simulations.It is worth noting that in the zenith angle discretizations (to be discussed in Section 3.3), both the incident angle θ1 and new propagation angle θ2 need to be discretized and the same discretization scheme must be applied to both angles.Such a requirement is posed by the nature of the Markov chain approximation, in that the transition states cannot be re-defined during the Markov process.Similarly, the resulting transmitted and reflected angular distributions tend to follow the same discretization that one sets at the beginning.Therefore, direct comparison among different discretization schemes will not offer sizably meaningful insights.Instead, the Markov calculations will be compared against one common benchmark with the results from the Monte Carlo simulations.Transmitted probabilities were normalized with a unit of 1° in this work to make sure that the summation of the probability over the range from 0° to 180° was unity minus the percentage of ballistic photons (transmitted photons without scattering, not discussed in this work).The same method was applied for all results presented in this work.
Figure 3 depicts Markov chain predictions of the transmitted photon angular distributions compared with those from the Monte Carlo simulations.In the simulations, uniform OD distribution was assumed, and the total OD was set to either 2 or 5. On the x axis, the transmitted angle from 0° to 90° indicates that the photons leave the turbid slab from the bottom exit plane, whilst the transmitted angle from 90° to 180° represents the photons that are reflected and leave the slab from the entrance plane.Since the range of the OD in the simulations was within the intermediate regime, the transmitted angular distribution exhibited both features from the phase function and multiple scattering averaging effect, as reflected by the complex angular distribution structures, as reflected in Figure 3. From the results, it is evident that Markov chain approximations agree well with Monte Carlo simulation results with satisfactory fidelity, for all the six cases considered in this work.Figure 4 shows the absolute relative error of the Markov chain approximations compared to the Monte Carlo simulations for the results shown in Figure 3, to further study the performance of Markov chain approximations.Note that other algorithms, such as the random walk [7], adding doubling [34], and the imperialist competitive algorithm [35] are also frequently used to evaluate such computational tasks.Since we have incorporated such comparisons in our previous work [22], we will not elaborate on such comparisons in this work.The absolute relative error in this paper is defined as: we will not elaborate on such comparisons in this work.The absolute relative error in this paper is defined as: The absolute relative error as a function of the transmitted angle is plotted in logarithm scale in Figure 4. We used Red dash lines to indicate an error threshold of 1.0%, and blue dash lines for a threshold of 0.1%.It is worth noting that we did not use the arbitrarily chosen thresholds to define any regime, but rather to evaluate the performance of each discretization scheme as a reference.From the data, several observations could be made.Firstly, the absolute relative error for all the cases was generally quite small (mostly below 0.1% for phase function 1 and phase function 2, below 1% for phase function 3).The significant error from the phase function 3 case implied that a more complex phase function would significantly increase the error during calculations.Furthermore, the maximum error would increase with a higher OD.More specifically, when the OD increased from 2 to 5, the maximum error (occurring around the transmission angle of 90°) increased from 0.6% to 6.0% in this investigation.The absolute relative error as a function of the transmitted angle is plotted in logarithm scale in Figure 4. We used Red dash lines to indicate an error threshold of 1.0%, and blue dash lines for a threshold of 0.1%.It is worth noting that we did not use the arbitrarily chosen thresholds to define any regime, but rather to evaluate the performance of each discretization scheme as a reference.From the data, several observations could be made.Firstly, the absolute relative error for all the cases was generally quite small (mostly below 0.1% for phase function 1 and phase function 2, below 1% for phase function 3).The significant error from the phase function 3 case implied that a more complex phase function would significantly increase the error during calculations.Furthermore, the maximum error would increase with a higher OD.More specifically, when the OD increased from 2 to 5, the maximum error (occurring around the transmission angle of 90 • ) increased from 0.6% to 6.0% in this investigation.

Optical Depth Discretizations
To solve the two challenges mentioned above, which are: (1) Increased error around 90 • with higher OD and (2) increased error from modeling media with complex phase functions, various discretization investigations were performed.Before optimizing the discretization schemes, a theoretical explanation for interpreting the significant error near 90 • was necessary.We considered a single photon that was scattered, then transmitted without being scattered again (the last scattering event).For this photon, if it is assumed that the optical depth from the scattering plane to the transmitted plane is to be OD t , then the probability of this photon being transmitted can be expressed as: Additionally, if we assumed that a disturbance/error ε of OD t takes place during the simulation, the relative change caused by this disturbance can be written as: In Markov chain discretizations, the disturbance ε would have a maximum value equal to the optical depth of the specific layer.As shown in Equation ( 8), when ε << cos θ t , exp(−ε/cos θ t )~1.0, the OD discretization error is less significant.However, when θ t is close to 90 • , cos θ t ~0.0, and exp(−ε/cos θ t ) no longer yields a value close to 1.0, this would cause significant errors in this range.It is also worth noting that when cos θ t ~0.0, the value of P(OD t ) will also be ~0.0 except when OD t has a sufficiently small value.Therefore, the probability for photons that scatter far away from the exit plane with a θ t ~90 • , then this is too small to introduce significant errors in the angular distribution predictions.Such observations imply transmitted photons with θ t ~90 • are mostly likely from the last layer (transmitted photons) or the first layer (reflected photons).Therefore, care should be taken when modeling these two cases.
In this work, two different non-uniform discretization schemes were examined.In the uniform OD discretization, each layer had 1% of the total OD (100 turbid slabs assumed).Non-uniform Equation ( 9) emphasized on setting the OD value of the first layer and the last layer, which can be expressed by: Non-uniform Equation ( 10) assumed an exponential OD discretization to vary the whole discretization, which can be expressed as: where a and b are modeling constants.Various values of a and b were examined to find the optimal values to yield the best performance.Although it was expected that a relatively small a and large b (which yields finer resolution near the exit planes) can minimize the error near 90 • , such a proposition might not be completely true.In the simulations, we found that optimal values of a and b exist which could minimize the relative error near 90 • .The relative error near 90 • can be minimized near a = 0.022 for all cases, as can be seen in Figure 5; and b = 0.03 for OD = 5 cases, as reflected in Figure 6.With discretization scheme 2, the relative errors for OD = 2 cases were not successfully reduced with any value of b.For discretization scheme 1, the maximum relative error near 90 • dropped from 0.6% to for all cases, as can be seen in Figure 5; and b = 0.03 for OD = 5 cases, as reflected in Figure 6.With discretization scheme 2, the relative errors for OD = 2 cases were not successfully reduced with any value of b.For discretization scheme 1, the maximum relative error near 90° dropped from 0.6% to 0.3% for OD = 2 cases, and it dropped from 6% to 0.8% for OD = 5 cases.For discretization scheme 2, the error dropped from 6% to 0.3% for OD = 5 cases.These observations verified that proper discretization can help improve computational accuracy.Upon closer examination, it could be readily observed that with the optimal value of either a or b, the relative error was minimized when the OD in the first layer and last layer was close to 0.022 (b = 0.03 yields an OD close to 0.022 in the first and last layer).For example, changing the OD in these layers to 0.022 reduced the error from 6% to 0.8%, and adopting an exponential scheme further reduced the error from 0.8% to 0.3%.Further decreasing the OD in the two layers would increase the maximum error.Scheme 2 did not work for OD = 2 cases because in Scheme 2, the OD in the first and the last layer would always be less than 0.02, and therefore, the optimal value of 0.022 could not be reached.Other cases with higher OD, such as OD = 10 and OD = 15, were also examined and similar observations were made.The authors admit that the rationale behind this optimal value of OD = 0.022 for the first layer and the last layer is not fully understood as yet, although the importance of the first and last layer has been shown.Further analytical investigations regarding this optimal value/scheme will be incorporated in future works.
In existing experiments, transmitted photons in θt~90° could not be readily accessed.As such, the relative error in the detectable regimes was more significant for potential spectroscopy methods, as least at present.For the sake of discussion, this angular range was set to 0-30° (transmitted photons) and 150-180° (reflected photons).The averaged absolute relative error in this range was used to evaluate the performance of different OD discretization schemes.Figure 7 depicts the averaged error with different OD discretization schemes.As seen in Figure 7, averaged absolute relative errors for all cases were quite low, from 0.05% to 0.3% with various phase functions and different ODs.It is also worth noting that the Markov chain approximation error for phase function 3 was significantly higher than phase function 1 and phase function 2. Such a discrepancy was believed to be caused by the discretization of the phase function, since phase function 3 was more spatially complex than the other two.An appropriate zenith angle discretization scheme may be needed for challenging phase functions, as will be discussed in Section 3.3 of this work.Finally, it is evident that the required averaged error was not remarkably affected by OD discretization.This observation was also verified with numerical repeats with other boundary conditions.Although OD discretization seemed to have a limited impact in the angular range of 0-30° and 150-180°, it was still recommended to adopt the scheme in setting the OD of the first and last layer to a value of 0.022 to incorporate likely measurements near 90°.Upon closer examination, it could be readily observed that with the optimal value of either a or b, the relative error was minimized when the OD in the first layer and last layer was close to 0.022 (b = 0.03 yields an OD close to 0.022 in the first and last layer).For example, changing the OD in these layers to 0.022 reduced the error from 6% to 0.8%, and adopting an exponential scheme further reduced the error from 0.8% to 0.3%.Further decreasing the OD in the two layers would increase the maximum error.Scheme 2 did not work for OD = 2 cases because in Scheme 2, the OD in the first and the last layer would always be less than 0.02, and therefore, the optimal value of 0.022 could not be reached.Other cases with higher OD, such as OD = 10 and OD = 15, were also examined and similar observations were made.The authors admit that the rationale behind this optimal value of OD = 0.022 for the first layer and the last layer is not fully understood as yet, although the importance of the first and last layer has been shown.Further analytical investigations regarding this optimal value/scheme will be incorporated in future works.
In existing experiments, transmitted photons in θ t ~90 • could not be readily accessed.As such, the relative error in the detectable regimes was more significant for potential spectroscopy methods, as least at present.For the sake of discussion, this angular range was set to 0-30 • (transmitted photons) and 150-180 • (reflected photons).The averaged absolute relative error in this range was used to evaluate the performance of different OD discretization schemes.Figure 7 depicts the averaged error with different OD discretization schemes.As seen in Figure 7, averaged absolute relative errors for all cases were quite low, from 0.05% to 0.3% with various phase functions and different ODs.It is also worth noting that the Markov chain approximation error for phase function 3 was significantly higher than phase function 1 and phase function 2. Such a discrepancy was believed to be caused by the discretization of the phase function, since phase function 3 was more spatially complex than the other two.An appropriate zenith angle discretization scheme may be needed for challenging phase functions, as will be discussed in Section 3.3 of this work.Finally, it is evident that the required averaged error was not remarkably affected by OD discretization.This observation was also verified with numerical repeats with other boundary conditions.Although OD discretization seemed to have a limited impact in the angular range of 0-30 • and 150-180 • , it was still recommended to adopt the scheme in setting the OD of the first and last layer to a value of 0.022 to incorporate likely measurements near 90 • .

Zenith Angle Discretizations
In this section, different zenith angle discretization schemes were tested.As previously mentioned, the zenith angle of photon propagation trajectory (from 0° to 180°) was divided into 180 intervals.Three different discretization schemes were utilized: (1) Uniform discretization.In this approach, each interval covered a uniform 1°; (2) Forward emphasizing approach.In this scheme, a finer resolution was applied for a small zenith angle range, for example, from 0° to 30.0°, and a coarse resolution was used for the rest of the angular range; and (3) Forward and backward emphasizing approach.In this scheme, a finer resolution was applied for both the forward and backward zenith angles, for example, from 0° to 15.0° and from 165° to 180.0°, and a coarse resolution was used for the rest of the angular range.The OD was assumed to be uniform and the value of the OD was also set to be either 2 or 5.
Figure 8 depicts the absolute relative error of the Markov chain approximations compared with the Monte Carlo simulations with zenith angle discretization scheme 2. To show the performance of the zenith angle discretization scheme, the phase function 3 was selected for its complex angular structure.For zenith angle discretization scheme 2, the angular resolution was 0.25° for θ < 30° and 2.5° for θ > 30°.As such, 180 meshes were used to cover the angular range from 0° to 180°, which kept the computational cost the same as the first discretization scheme.Uniform OD discretization was utilized in Figure 8.It can be seen in Figure 8 that the absolute relative error for both OD = 2 and OD = 5 cases was reduced from 1% to 0.2% in the angular range of 0° to 30°.Therefore, a proper zenith angle discretization scheme can indeed improve the predictive performance in certain regimes.However, the trade-off is the sacrifice of prediction accuracy in the angular regime from 30° to 180°.As can be seen, the relative error near 90° and 180° was more significant.Furthermore, adopting a different OD discretization scheme could not remarkably reduce the peak error for the new zenith angle discretization scheme.a and b values recommended in Section 3.2 were still valid; however, the values only reduced the peak error from 6% to 4.5% (OD = 2) and 3% to 2% (OD = 5).This observation suggested that with a coarse zenith angle mesh, angular discretization was the dominating source of the computing error, especially for the regimes with coarse mesh.However, it should be noted that experimental measurement can only be achieved within a limited angular regime, i.e., from 0° to 30°.Improving the accuracy in this regime was critical for practical applications whilst preserving the accuracy in the range from 30° to 180° was not necessary.Therefore, a proper zenith angle discretization was meaningful for angular distribution-based algorithms and potential spectroscopy methods.

Zenith Angle Discretizations
In this section, different zenith angle discretization schemes were tested.As previously mentioned, the zenith angle of photon propagation trajectory (from 0 • to 180 • ) was divided into 180 intervals.Three different discretization schemes were utilized: (1) Uniform discretization.In this approach, each interval covered a uniform 1 • ; (2) Forward emphasizing approach.In this scheme, a finer resolution was applied for a small zenith angle range, for example, from 0 • to 30.0 • , and a coarse resolution was used for the rest of the angular range; and (3) Forward and backward emphasizing approach.In this scheme, a finer resolution was applied for both the forward and backward zenith angles, for example, from 0 • to 15.0 • and from 165 • to 180.0 • , and a coarse resolution was used for the rest of the angular range.The OD was assumed to be uniform and the value of the OD was also set to be either 2 or 5.
Figure 8 depicts the absolute relative error of the Markov chain approximations compared with the Monte Carlo simulations with zenith angle discretization scheme 2. To show the performance of the zenith angle discretization scheme, the phase function 3 was selected for its complex angular structure.For zenith angle discretization scheme 2, the angular resolution was 0.25 • for θ < 30 • and 2.5 • for θ > 30 • .As such, 180 meshes were used to cover the angular range from 0 • to 180 • , which kept the computational cost the same as the first discretization scheme.Uniform OD discretization was utilized in Figure 8.It can be seen in Figure 8 that the absolute relative error for both OD = 2 and OD = 5 cases was reduced from 1% to 0.2% in the angular range of 0 • to 30 • .Therefore, a proper zenith angle discretization scheme can indeed improve the predictive performance in certain regimes.However, the trade-off is the sacrifice of prediction accuracy in the angular regime from 30 • to 180 • .As can be seen, the relative error near 90 • and 180 • was more significant.Furthermore, adopting a different OD discretization scheme could not remarkably reduce the peak error for the new zenith angle discretization scheme.a and b values recommended in Section 3.2 were still valid; however, the values only reduced the peak error from 6% to 4.5% (OD = 2) and 3% to 2% (OD = 5).This observation suggested that with a coarse zenith angle mesh, angular discretization was the dominating source of the computing error, especially for the regimes with coarse mesh.However, it should be noted that experimental measurement can only be achieved within a limited angular regime, i.e., from 0 • to 30 • .Improving the accuracy in this regime was critical for practical applications whilst preserving the accuracy in the range from 30 • to 180 • was not necessary.Therefore, a proper zenith angle discretization was meaningful for angular distribution-based algorithms and potential spectroscopy methods.Figure 9 closely examines the performance of zenith angle discretization scheme 2 in the range from 0° to 30°.In this comparison, phase function 3 was selected and OD = 5 for Figure 9. Panel (a) of Figure 9 shows the transmitted angular distribution predictions for both uniform zenith angle distribution (scheme 1) and non-uniform zenith angle distribution (scheme 2).As seen, a uniform angular distribution is not sufficient to resolve the fine structure of the angular distribution predicted by the Monte Carlo simulation (which is more obvious near 3°).In comparison, non-uniform scheme 2 successfully resolved the profile of the transmitted angular distribution.Panel (b) of Figure 9 plots the absolute relative error together with the phase function 3. It is evident that for angular discretization scheme 1, the profile of the absolute relative error resembled the profile of the phase function, indicating that the error was mainly caused by insufficient resolution with regards to the phase function.In contrast, the error curve from angular distribution discretization 2 eliminated such similarities, indicating that the error caused by angular discretization was less apparent.Finally, Figure 10 examines the absolute relative error of Markov chain approximations compared with Monte Carlo simulations with zenith angle discretization scheme 3. Phase function 3 was selected and the angular resolution was 0.25° for θ < 15° and θ > 165°, and 2.5° for anywhere else.It can be seen in Figure 10 that the absolute relative error for both OD = 2 and OD = 5 cases was again reduced from 1% to 0.2% in the angular range of 0° to 15°.However, the peak error near 180° was not Figure 9 closely examines the performance of zenith angle discretization scheme 2 in the range from 0 • to 30 • .In this comparison, phase function 3 was selected and OD = 5 for Figure 9. Panel (a) of Figure 9 shows the transmitted angular distribution predictions for both uniform zenith angle distribution (scheme 1) and non-uniform zenith angle distribution (scheme 2).As seen, a uniform angular distribution is not sufficient to resolve the fine structure of the angular distribution predicted by the Monte Carlo simulation (which is more obvious near 3 • ).In comparison, non-uniform scheme 2 successfully resolved the profile of the transmitted angular distribution.Panel (b) of Figure 9 plots the absolute relative error together with the phase function 3. It is evident that for angular discretization scheme 1, the profile of the absolute relative error resembled the profile of the phase function, indicating that the error was mainly caused by insufficient resolution with regards to the phase function.In contrast, the error curve from angular distribution discretization 2 eliminated such similarities, indicating that the error caused by angular discretization was less apparent.Figure 9 closely examines the performance of zenith angle discretization scheme 2 in the range from 0° to 30°.In this comparison, phase function 3 was selected and OD = 5 for Figure 9. Panel (a) of Figure 9 shows the transmitted angular distribution predictions for both uniform zenith angle distribution (scheme 1) and non-uniform zenith angle distribution (scheme 2).As seen, a uniform angular distribution is not sufficient to resolve the fine structure of the angular distribution predicted by the Monte Carlo simulation (which is more obvious near 3°).In comparison, non-uniform scheme 2 successfully resolved the profile of the transmitted angular distribution.Panel (b) of Figure 9 plots the absolute relative error together with the phase function 3. It is evident that for angular discretization scheme 1, the profile of the absolute relative error resembled the profile of the phase function, indicating that the error was mainly caused by insufficient resolution with regards to the phase function.In contrast, the error curve from angular distribution discretization 2 eliminated such similarities, indicating that the error caused by angular discretization was less apparent.Finally, Figure 10 examines the absolute relative error of Markov chain approximations compared with Monte Carlo simulations with zenith angle discretization scheme 3. Phase function 3 was selected and the angular resolution was 0.25° for θ < 15° and θ > 165°, and 2.5° for anywhere else.It can be seen in Figure 10 that the absolute relative error for both OD = 2 and OD = 5 cases was again reduced from 1% to 0.2% in the angular range of 0° to 15°.However, the peak error near 180° was not seen in Figure 10 that the absolute relative error for both OD = 2 and OD = 5 cases was again reduced from 1% to 0.2% in the angular range of 0 • to 15 • .However, the peak error near 180 • was not reduced significantly.Two observations were made regarding this phenomenon: (1) The significant error of ~1% near 180 • only occured when phase function 3 was used, as can be seen in Figures 4-6, indicating this was an issue related to angular discretization; (2) The error near 0 • was remarkably reduced by discretization scheme 3 whilst the error near 180 • did not, which reflected that the error caused by coarse meshes in the range from 15 • to 165 • had been accumulated with multiple scattering effects.reduced significantly.Two observations were made regarding this phenomenon: (1) The significant error of ~1% near 180° only occured when phase function 3 was used, as can be seen in Figures 4-6, indicating this was an issue related to angular discretization; (2) The error near 0° was remarkably reduced by discretization scheme 3 whilst the error near 180° did not, which reflected that the error caused by coarse meshes in the range from 15° to 165° had been accumulated with multiple scattering effects.

Computational Cost Optimizations
One feature of the Markov chain approximation is that the computational cost remains relatively constant with different discretization schemes, phase functions, or different ODs; given that the size of transition matrix is fixed (which is determined by the number of meshes in the system).In the previous works, the computational cost for anisotropic multiple scattering was reported as ~8 min.Such a computational cost is obviously superior compared to many Monte Carlo simulations, especially with larger OD values, because with more scattering events, the computational cost for a single photon increases considerably.
Despite that, it is still desirable to further decrease the computational cost of the original Markov chain program, so that more computationally efficient reconstruction algorithms can be made possible.To analyze the computational costs, the approximation process was broken down into three parts.The three parts were: (1) Solving Equations ( 3) and ( 4), which builds the connection between the propagation angle before scattering (θ1) and after scattering (θ2); (2) Solving Equations ( 1) and ( 2), which incorporates the OD and part 1 to form the transition matrix; and (3) Perform a matrix inversion and obtain the results.For Markov chain methods, part (1) and part (3) consume most of the computational cost in the simulation, with each process consuming a computational time of more than 100 s.
To reduce the computational cost, for part (1), it is worth noting that P (θi, θj, n) is only dependent on the phase function in this specific layer where the scattering happens, which means P (θi, θj, n) can be calculated in advance and a 180 by 180 matrix can be stored to represent each phase function.When the transition matrix is needed during the calculation, the simulation can load the matrices directly with negligible computational costs.This is also an added benefit of the Markov chain approximations when compared to the Monte Carlo method, because the trajectory of each photon must be calculated in every evaluation in Monte Carlo simulations.
In part (3), the computational cost for matrix inversion is also significant, mostly due to the computational complexity of inversing a significantly large transition matrix (18,000 by 18,000).From the definition of the fundamental matrix:

Computational Cost Optimizations
One feature of the Markov chain approximation is that the computational cost remains relatively constant with different discretization schemes, phase functions, or different ODs; given that the size of transition matrix is fixed (which is determined by the number of meshes in the system).In the previous works, the computational cost for anisotropic multiple scattering was reported as ~8 min.Such a computational cost is obviously superior compared to many Monte Carlo simulations, especially with larger OD values, because with more scattering events, the computational cost for a single photon increases considerably.
Despite that, it is still desirable to further decrease the computational cost of the original Markov chain program, so that more computationally efficient reconstruction algorithms can be made possible.To analyze the computational costs, the approximation process was broken down into three parts.The three parts were: (1) Solving Equations ( 3) and ( 4), which builds the connection between the propagation angle before scattering (θ 1 ) and after scattering (θ 2 ); (2) Solving Equations ( 1) and ( 2), which incorporates the OD and part 1 to form the transition matrix; and (3) Perform a matrix inversion and obtain the results.For Markov chain methods, part (1) and part (3) consume most of the computational cost in the simulation, with each process consuming a computational time of more than 100 s.
To reduce the computational cost, for part (1), it is worth noting that P (θ i , θ j , n) is only dependent on the phase function in this specific layer where the scattering happens, which means P (θ i , θ j , n) can be calculated in advance and a 180 by 180 matrix can be stored to represent each phase function.When the transition matrix is needed during the calculation, the simulation can load the matrices directly with negligible computational costs.This is also an added benefit of the Markov chain approximations when compared to the Monte Carlo method, because the trajectory of each photon must be calculated in every evaluation in Monte Carlo simulations.
In part (3), the computational cost for matrix inversion is also significant, mostly due to the computational complexity of inversing a significantly large transition matrix (18,000 by 18,000).From the definition of the fundamental matrix: where N is the fundamental matrix, and T is the transition matrix obtaining the transition probabilities.
Note that the angular distribution can be expressed as (based on Markov chain Theory): where Q total is the sought angular distribution, R is the absorption matrix, and P is the initial scattering distribution.One way to reduce the computational cost of the Markov chain inversion is to use the expanded form of N to solve for N, instead of obtaining the inverse of (I − T) directly.If it is defined that: Q t stands for the angular distribution of photons that experienced exactly t times of scattering events (which is defined as the scattering order of this transmitted photon), before exiting the turbid slab.
Then an iterative relationship can be easily derived: Additionally, the iteration process is stopped when the number of photons transmitted with a scattering order of t is sufficiently small compared to the total number of photons transmitted.In this study, the ratio of Q t and Q total was obtained during each iteration, and the program would stop if the maximum value of the ratio vector was smaller than 10 −4 .
Figure 11 shows the computational cost breakdown for solving for the angular distribution, with and without the iterative approach.Simulations with various OD were performed and the Phase Function 2 was selected.As can be seen in Figure 11, when utilizing the iteration process, the computational cost for the inversion process is reduced from 147 s to 4 s (OD = 2) or 7 s (OD = 5).By incorporating this approach with the pre-calculated phase function P(θ i , θ j , n), the total computational cost for executing Markov chain approximations can be reduced from 463 s to 15 s (for OD = 2) and 18 s (for OD = 5).
where N is the fundamental matrix, and T is the transition matrix obtaining the transition probabilities.Note that the angular distribution can be expressed as (based on Markov chain Theory): where Qtotal is the sought angular distribution, R is the absorption matrix, and P′ is the initial scattering distribution.One way to reduce the computational cost of the Markov chain inversion is to use the expanded form of N to solve for N, instead of obtaining the inverse of (I − T) directly.If it is defined that: Qt stands for the angular distribution of photons that experienced exactly t times of scattering events (which is defined as the scattering order of this transmitted photon), before exiting the turbid slab.
Then an iterative relationship can be easily derived: Additionally, the iteration process is stopped when the number of photons transmitted with a scattering order of t is sufficiently small compared to the total number of photons transmitted.In this study, the ratio of Qt and Qtotal was obtained during each iteration, and the program would stop if the maximum value of the ratio vector was smaller than 10 −4 .
Figure 11 shows the computational cost breakdown for solving for the angular distribution, with and without the iterative approach.Simulations with various OD were performed and the Phase Function 2 was selected.As can be seen in Figure 11, when utilizing the iteration process, the computational cost for the inversion process is reduced from 147 s to 4 s (OD = 2) or 7 s (OD = 5).By incorporating this approach with the pre-calculated phase function P(θi, θj, n), the total computational cost for executing Markov chain approximations can be reduced from 463 s to 15 s (for OD = 2) and 18 s (for OD = 5).Figure 12 examines the stopping scattering order (SO) and Markov Chan computational costs as compared to the Monte Carlo simulations.The same CPU was utilized for Markov chain approximation and Monte Carlo simulations and 2.5 billion photons were sent in each case.Since MATLAB matrix manipulations involve parallel computing, the computational cost of the Monte Carlo simulation was scaled accordingly.From Figure 11, it is evident that with a higher OD, the   12 examines the stopping scattering order (SO) and Markov Chan computational costs as compared to the Monte Carlo simulations.The same CPU was utilized for Markov chain approximation and Monte Carlo simulations and 2.5 billion photons were sent in each case.Since MATLAB matrix manipulations involve parallel computing, the computational cost of the Monte Carlo simulation was scaled accordingly.From Figure 11, it is evident that with a higher OD, the computational cost of the modified Markov chain method was also higher because the stopping scattering order t was also higher.The stopping scattering order was roughly equal to the OD multiplied by 12, representing a quasi-linear relationship.The computational time for the modified Markov inversion process ranged from 4 s for OD = 2 cases to 41 s for OD = 30 cases, which was fast enough even when global optimization algorithms were utilized.In comparison, the Monte Carlo simulation cost 3000 s for OD = 2 cases and 61,000 s for OD = 30 cases, and such computational costs would inevitably discourage any meaningful inversion process.It is worth noting that varying the discretization could slightly change the stopping scattering order, but the relative difference was within 10% as the authors examined.Therefore, it was established that the computational cost estimation in Figure 12 is valid for various discretization schemes.
Appl.Sci.2018, 8, x FOR PEER REVIEW 15 of 18 computational cost of the modified Markov chain method was also higher because the stopping scattering order t was also higher.The stopping scattering order was roughly equal to the OD multiplied by 12, representing a quasi-linear relationship.The computational time for the modified Markov inversion process ranged from 4 s for OD = 2 cases to 41 s for OD = 30 cases, which was fast enough even when global optimization algorithms were utilized.In comparison, the Monte Carlo simulation cost 3000 s for OD = 2 cases and 61,000 s for OD = 30 cases, and such computational costs would inevitably discourage any meaningful inversion process.It is worth noting that varying the discretization could slightly change the stopping scattering order, but the relative difference was within 10% as the authors examined.Therefore, it was established that the computational cost estimation in Figure 12 is valid for various discretization schemes.Figure 13 depicts the absolute relative error obtained with the modified Markov inversion process.Uniform zenith angle discretization was used and the OD discretization scheme 1 was chosen.As such, the conditions in Figure 13 are the same as Figure 5, except for where Markov inversion process was chosen.As can be seen in all six cases, from Panel (a) to Panel (f), the absolute relative error with the modified Markov inversion process was at the same level as the original Markov inversion with a higher computational cost, therefore verifying the application of the modified Markov inversion process, as well as the recommended value of 10 −4 as the stopping criterion.Note that this recommended stopping criterion is subject to change depending on the desired accuracy, and the program can be stopped even earlier if high accuracy is not required.Furthermore, with the modified method, transmitted angular distribution at each scattering order can also be derived without additional computational costs, thus this method is more efficient than directly performing matrix inversions in solving multiple scattering problems.
Finally, it is worth mentioning that the modified Markov inversion process can also lead to the use of a greater transition matrix, thus enabling finer resolution for Markov chain process.In the original process, the fundamental matrix needs to be solved by inversing the transition matrix; therefore typically, the simulation must store the whole transition matrix for matrix inversion, which would cause RAM usage related problems (i.e., matrix size exceeds available RAM).However, with the modified inversion process, there is no need to register the complete matrix in the RAM.Vectors can be easily calculated from P (θi, θj, n) with far less computational cost.Although this approach could increase the computational cost for such calculations, this method validates the use of the Markov chain approximation with higher dimensions or finer resolutions.Figure 13 depicts the absolute relative error obtained with the modified Markov inversion process.Uniform zenith angle discretization was used and the OD discretization scheme 1 was chosen.As such, the conditions in Figure 13 are the same as Figure 5, except for where Markov inversion process was chosen.As can be seen in all six cases, from Panel (a) to Panel (f), the absolute relative error with the modified Markov inversion process was at the same level as the original Markov inversion with a higher computational cost, therefore verifying the application of the modified Markov inversion process, as well as the recommended value of 10 −4 as the stopping criterion.Note that this recommended stopping criterion is subject to change depending on the desired accuracy, and the program can be stopped even earlier if high accuracy is not required.Furthermore, with the modified method, transmitted angular distribution at each scattering order can also be derived without additional computational costs, thus this method is more efficient than directly performing matrix inversions in solving multiple scattering problems.
Finally, it is worth mentioning that the modified Markov inversion process can also lead to the use of a greater transition matrix, thus enabling finer resolution for Markov chain process.In the original process, the fundamental matrix needs to be solved by inversing the transition matrix; therefore typically, the simulation must store the whole transition matrix for matrix inversion, which would cause RAM usage related problems (i.e., matrix size exceeds available RAM).However, with the modified inversion process, there is no need to register the complete matrix in the RAM.Vectors can be easily calculated from P (θ i , θ j , n) with far less computational cost.Although this approach could increase the computational cost for such calculations, this method validates the use of the Markov chain approximation with higher dimensions or finer resolutions.

Conclusions
This research investigated many practical aspects in using the Markov chain approximation to solve photon multiple scattering through turbid slabs.Although uniform OD discretization and zenith angle discretization have already ensured an accurate prediction of the transmitted and reflected angular distribution, it has been shown in this work that non-uniform discretization can be utilized in Markov chain approximations to model multiple scattering, and the improvement in simulation performance was easily seen.
A modified matrix inversion scheme was also proposed in this work.By converting the matrix inversion process into an iterative process, the computational cost of solving a Markov chain inversion process was reduced from 180 s to 10 s.By pre-processing the computational tasks that are dependent only on the phase function in the specific layer and incorporating the iterative process, the total computational cost was successfully reduced from 8 min to about 30 s, which will greatly contribute to the development of reconstruction/inversion algorithms.
Overall, with both the non-uniform discretization scheme and the iterative process, the capacity of the proposed Markov chain method was improved remarkably, both in the aspects of accuracy and computational cost.The authors believed that proper OD discretization and zenith angle discretization were crucial for improving the performance of the proposed Markov chain method by adopting finer meshes or coarser meshes at locations with different levels of importance based on practical needs.As such, the computational cost is reduced with less total states and a similar level of accuracy.The iterative process further reduces the computational cost of Markov chain evaluations, and it enables Markov chain approximations with a much higher number of states.These methods in total open the door to

Conclusions
This research investigated many practical aspects in using the Markov chain approximation to solve photon multiple scattering through turbid slabs.Although uniform OD discretization and zenith angle discretization have already ensured an accurate prediction of the transmitted and reflected angular distribution, it has been shown in this work that non-uniform discretization can be utilized in Markov chain approximations to model multiple scattering, and the improvement in simulation performance was easily seen.
A modified matrix inversion scheme was also proposed in this work.By converting the matrix inversion process into an iterative process, the computational cost of solving a Markov chain inversion process was reduced from 180 s to 10 s.By pre-processing the computational tasks that are dependent only on the phase function in the specific layer and incorporating the iterative process, the total computational cost was successfully reduced from 8 min to about 30 s, which will greatly contribute to the development of reconstruction/inversion algorithms.
Overall, with both the non-uniform discretization scheme and the iterative process, the capacity of the proposed Markov chain method was improved remarkably, both in the aspects of accuracy and computational cost.The authors believed that proper OD discretization and zenith angle discretization were crucial for improving the performance of the proposed Markov chain method by adopting finer meshes or coarser meshes at locations with different levels of importance based on practical needs.As such, the computational cost is reduced with less total states and a similar level of accuracy.The iterative process further reduces the computational cost of Markov chain evaluations, and it enables Markov chain approximations with a much higher number of states.These methods in total open the door to accurate and computationally affordable inversion methods by reducing the computational costs of a single cost function evaluation.Such methods would lead to further improved diagnostics applications, such as spray diagnostics in the dense regime and bio-tissue characterization, which will be incorporated in future works.

Figure 2 .
Figure 2. Phase functions used in this study.

Figure 2 .
Figure 2. Phase functions used in this study.

Figure 2 .
Figure 2. Phase functions used in this study.

Figure 7 .
Figure 7. Averaged absolute relative error of Markov chain angular distribution predictions compared to Monte Carlo in different cases.Angular range was set to 0-30° (transmitted photons) and 150-180° (reflected photons).

Figure 7 .
Figure 7. Averaged absolute relative error of Markov chain angular distribution predictions compared to Monte Carlo in different cases.Angular range was set to 0-30 • (transmitted photons) and 150-180 • (reflected photons).

Figure 9 .
Figure 9. Solutions from Markov chain approximations compared to Monte Carlo simulations with zenith angle discretization scheme 2. Phase function 3 was selected and OD = 5.0.Panel (a).Transmitted angular distribution, Panel (b).Absolute relative error.

Figure 9 .
Figure 9. Solutions from Markov chain approximations compared to Monte Carlo simulations with zenith angle discretization scheme 2. Phase function 3 was selected and OD = 5.0.Panel (a).Transmitted angular distribution, Panel (b).Absolute relative error.

Figure 9 .
Figure 9. Solutions from Markov chain approximations compared to Monte Carlo simulations with zenith angle discretization scheme 2. Phase function 3 was selected and OD = 5.0.Panel (a).Transmitted angular distribution, Panel (b).Absolute relative error.Finally, Figure 10 examines the absolute relative error of Markov chain approximations compared with Monte Carlo simulations with zenith angle discretization scheme 3. Phase function 3 was selected and the angular resolution was 0.25 • for θ < 15 • and θ > 165 • , and 2.5 • for anywhere else.It can be

Figure 11 .
Figure 11.Histogram of computational costs for the original and improved Markov chain approximations.

Figure 11 .
Figure 11.Histogram of computational costs for the original and improved Markov chain approximations.

Figure 12 .
Figure 12.Modified Markov chain inversion computational cost and stopping scattering order as a function of the OD, in comparison to Monte Carlo simulations.

Figure 12 .
Figure 12.Modified Markov chain inversion computational cost and stopping scattering order as a function of the OD, in comparison to Monte Carlo simulations.