A closed-form method for calculating the angular distribution of multiply scattered photons through isotropic turbid slabs

This paper develops a method for calculating the angular distribution (AD) of multiply scattered photons through isotropic turbid slabs. Extension to anisotropic scattering is also discussed. Previous studies have recognized that the AD of multiply scattered photons is critical for many applications, such as the design of imaging optics and estimation of image quality. This paper therefore develops a closed-from method that can accurately calculate the AD over a wide range of conditions. Other virtues of the method include its simplicity in implementation and its prospective for extension to anisotropic scattering. ©2011 Optical Society of America OCIS codes: (290.4210) Multiple scattering; (290.4020) Mie theory; (290.7050) Turbid media. References and links 1. J. C. Hebden, D. J. Hall, M. Firbank, and D. T. Delpy, “Time-Resolved Optical Imaging of a Solid TissueEquivalent Phantom,” Appl. Opt. 34(34), 8038–8047 (1995). 2. R. M. Measures, Laser Remote Sensing: Fundamentals and applications (Krieger Publishing Company, 1992). 3. M. A. Linne, M. Paciaroni, J. R. Gord, and T. R. Meyer, “Ballistic Imaging of The Liquid Core for a Steady Jet in Crossflow,” Appl. Opt. 44(31), 6627–6634 (2005). 4. J. A. Moon, P. R. Battle, M. Bashkansky, R. Mahon, M. D. Duncan, and J. Reintjes, “Achievable Spatial Resolution of Time-Resolved Transillumination Imaging Systems Which Utilize Multiply Scattered Light,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 53(1), 1142–1155 (1996). 5. J. A. Moon and J. Reintjes, “Image Resolution by Use of Multiply Scattered Light,” Opt. Lett. 19(8), 521–523 (1994). 6. J. A. Moon, R. Mahon, M. D. Duncan, and J. Reintjes, “Resolution Limits for Imaging through Turbid Media with Diffuse Light,” Opt. Lett. 18(19), 1591–1593 (1993). 7. M. Paciaroni and M. Linne, “Single-Shot, Two-Dimensional Ballistic Imaging Through Scattering Media,” Appl. Opt. 43(26), 5100–5109 (2004). 8. A. A. Kokhanovsky, “Analytical Solutions of Multiple Light Scattering Problems: A Review,” Meas. Sci. Technol. 13(3), 233–240 (2002). 9. A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering (Prentice Hall, Englewood Cliffs, New Jersey, 1991). 10. D. Contini, F. Martelli, and G. Zaccanti, “Photon Migration Through a Turbid Slab Described by a Model Based on Diffusion Approximation. 2. Theory,” Appl. Opt. 36(19), 4587–4599 (1997). 11. A. H. Gandjbakhche, G. H. Weiss, R. F. Bonner, and R. Nossal, “Photon Path-Length Distributions for Transmission through Optically Turbid Slabs,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 48(2), 810–818 (1993). 12. E. Berrocal, D. L. Sedarsky, M. E. Paciaroni, I. V. Meglinski, and M. A. Linne, “Laser Light Scattering in Turbid Media Part I: Experimental and Simulated Results for The Spatial Intensity Distribution,” Opt. Express 15(17), 10649–10665 (2007). 13. M. D. King and Harshvardhan, “Comparative Accuracy of Selected Multiple Scattering Approximations,” J. Atmos. Sci. 43(8), 784–801 (1986). 14. R. F. Bonner, R. Nossal, S. Havlin, and G. H. Weiss, “Model for Photon Migration in Turbid Biological Media,” J. Opt. Soc. Am. A 4(3), 423–432 (1987). 15. K.-N. Liou, An Introduction to Atmospheric Radiation (Academic Press Ltd., 2002). 16. I. M. Sobol, The Monte Carlo Method (The University of Chicago Press, 1967). 17. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo Programs of Polarized Light Transport into Scattering Media: Part I,” Opt. Express 13(12), 4420–4438 (2005). #151545 $15.00 USD Received 22 Jul 2011; revised 12 Oct 2011; accepted 1 Nov 2011; published 10 Nov 2011 (C) 2011 OSA 21 November 2011 / Vol. 19, No. 24 / OPTICS EXPRESS 23932 18. E. Berrocal, D. L. Sedarsky, M. E. Paciaroni, I. V. Meglinski, and M. A. Linne, “Laser Light Scattering in Turbid Media Part II: Spatial and Temporal Analysis of Individual Scattering Orders via Monte Carlo Simulation,” Opt. Express 17(16), 13792–13809 (2009). 19. Handbook of Mathematical Functions, with Formulas, Graphs, and Mathematical Tables (Dover, 1965).


Introduction
Multiple scattering (due to Mie scattering) represents a fundamental problem with applications in a wide spectrum of practical systems, ranging from imaging through biological tissues [1], remote sensing in the atmosphere [2], and development of laser diagnostics for dense sprays [3].Consequently, this problem has attracted a considerable amount of research efforts from different disciplines.In many of these applications, the angular distribution (AD) of the multiply scattered photons plays an important role.For example, the AD plays a critical role in designing the imaging system involving multiple scattering, such as the selection of the lens [4], estimation of the detection limit and resolution [5,6], and the application of angular filtering [7].Furthermore, practical applications typically involve an optimization facet, which requires the AD to be computed a great many times.Therefore, it is highly desirable to have a simple (preferably closed-form) method to calculate the AD efficiently.
Fundamentally, the multiple scattering problem requires one to solve the wave equation in a system of a large number of scatterers, which is extremely difficult [8].Under certain assumptions, the problem can be simplified and the so-called radiative transfer equation (RTE) can be derived [9].However, even the solution of the simplified RTE is difficult under most non-ideal conditions (e.g., realistic geometry).Thus various approximate solution techniques, both theoretical and numerical, have been developed [4,8,[10][11][12].The theoretical techniques seek analytical solutions to the RTE under certain approximations [13].For instance, under large optical thickness, the photon diffusion model [4,10] approximates the multiple scattering problem as a diffusion problem.Another example, the random walk model [11,14] approximates the multiple scattering problem as a random walk process.These models result in analytical expressions that can be used to predict the spatial and temporal distribution of multiply scattered photons.However, these analytical expressions can be regarded as semi-closed-form (at lease from the perspective of practical applications).For example, these expressions typically involve special functions and summation of infinite series [10,11,15], which complicates implementation.Furthermore, much of the analytical work did not target the AD specifically, and therefore application of the previous results to calculate the AD requires extension work, sometimes considerable extension work.
The numerical techniques solve the RTE computationally and a major limitation is the computation cost.The Monte Carlo (MC) technique provides a notable example of the numerical technique.The MC technique has been demonstrated as a powerful tool, capable of solving multiple scattering problems under realistic conditions and resolving quantities that are difficult to study experimentally [12,[16][17][18].However, the MC technique encounters considerable challenges in terms of computational cost, especially at relatively large optical depths.The challenge is especially acute for calculating the AD, because resolving the AD requires receiving sufficient photon counts in a small angular range out of the transmitted photons, which are already a small fraction of incident photons at large optical depth (illustrated by the results in Section 3).
The goal of this paper is therefore to develop a closed-form method to calculate the AD of multiply scattered photons.Closed-form expressions are derived to calculate the AD of photons transmitted from isotropic turbid slabs, based on which, the average cosine and fraction of photons within a given acceptance angle (numerical aperture) can be readily calculated.The results obtained from the method were compared against MC data, and good agreement is achieved over a wide range of optical depths.The method is simple to implement and highly efficient, suitable for application in a range of practical applications.Section 2 describes the method.Section 3 reports the results and discusses the extension of the current method to anisotropic scatterers.Finally, Section 4 summarizes the paper.

Description of the method
Figure 1 illustrates the schematic of the multiple scattering problem and the MC configuration used in this work.As shown in Panel (a), incident photons enter an infinitely large turbid slab as a pencil beam.The slab, with thickness Z, contains isotropic scatterers.Note that isotropic slabs in this work mean the scatterers are isotropic and the distribution of the scatterers in the slab is also isotropic (i.e., uniform).All length scales in this work are expressed in units of the mean scattering length, i.e., 1 S − Σ , where S Σ (mm-1) represents the scattering coefficient of the scattering media.The plane on which the incident photons enter is named the Incident Plane (IP), and the opposite plane named the Exit Plane (EP).After a certain number of scattering events, the incident photons exit either via the IP (as reflected photons) or via the EP (as transmitted photons).Here we focus on the AD of the transmitted photons, and extension of the method to calculate the AD of the reflected photons is straightforward.
The MC simulation in this work approximates the infinitely large slab by a rectangle slab, whose width and height are both 100 × of its thickness as shown in Panel (a).This approximation is justified by the fact that for all cases simulated, the fraction of photons leaked via the sides is less than 10 −6 .This work used a standard MC model as described in [16], and the implementation of the model largely follows [12].The MC model involves the following four major steps: 1) it sends incident photons into the IP one by one, 2) it determines the scattering direction and pathlength of the incident photon stochastically to calculate the location of the photon after a scattering event, 3) it repeats step 2) until the photon exits the slab, and 4) it registers the location and orientation of each photon that exits.The MC model terminates after a sufficient number of exit photons have been registered to determine the AD with acceptable statistical accuracy, which corresponds to a number of incident photons on the order of 10 7 in this work.The free scattering pathlength in the MC model was taken to be exponentially distributed, with the mean free pathlength equal to 1 S − Σ .A Cartesian coordinate system is defined as shown in Panel (b), where the z-axis is perpendicular to the EP.The origin is defined to be on the EP, and the positive of the z-axis points towards the IP.The scattering angle of a transmitted photon (θ) is defined as the angle formed relative to the negative z-axis after its last scattering event (illustrated by photon 1).Under these definitions, the AD function (Q T (θ)) of the transmitted photons is: where P is the AD of the transmitted photons at z, Qz the fraction of the transmitted photons that undergo their last scattering event at z, ( )sin  The rigorous derivation is of Q z is more complicated.Here we derive Q z under some simplifications, and then briefly discuss the rigorous derivation; though the results derived under these simplifications are already in reasonable agreement with MC data.The first simplification is to assume that the transmitted photons are always scattered at θ = 0 by the last scattering event (i.e., they all exit along the negative z-axis as photon 2. This appears to be a dramatic simplification, however, its validity on isotropic scatterers have been well demonstrated in previous work where a random walker on a lattice is used to model multiple scattering [11,14].Also note that this simplification is only invoked for the derivation of Q z (one step in the derivation of the AD), and does not imply that the AD is a delta-function.The second simplification assumes that all transmitted photons exit exactly on the EP.With these simplifications, Q z could be understood as the distribution of photons with respect to z after one scattering event along the positive z-axis if the photons were injected into the slab on the EP.With this understanding, because the photon pathlength is exponentially distributed, Q z is: where B is the normalization factor such that In this case, it can be shown that [11]: Substitution of Eqs. ( 2) through (4) into Eq.( 1) yields the final expression for Q T : (1/ cos 1) It can be seen that the above method can be extended to calculate the AD of reflected photons by re-deriving Q z .Our preliminary study suggests that, under similar simplifications, the Q z for reflected photons should be similar to Eq. ( 3) above.A thorough treatment of this extension is undergoing.

Results and discussions
Figure 3 shows the AD of the transmitted photons calculated using Eq. ( 5) compared to MC data.As shown in Fig. 3, the method developed above agrees with MC data across a wide range of OD.Two interesting observations can also be made: 1) the AD is peaked towards θ = 90° at small OD (Panel (a)), and 2) after OD becomes larger than 5, the AD becomes insensitive to OD (Panel (b)).The first observation can be explained by noting that at small OD (e.g., 0.1), photons undergo fewer scattering events to be transmitted and consequently are more likely to exit at large angles (if they do exit).The second observation can be explained by noting that at large OD, photons undergo a large number of scattering events and their directions become completely randomized, resulting in the AD's insensitivity.
Fig. 3. AD of transmitted photons calculated using Eq. ( 5) compared to MC data.The symbols represent the MC data (with ballistic photons excluded), the solid lines the calculation using Eq. ( 5).
With the validity of Eq. ( 5) confirmed, other properties of interest can be readily calculated.For example, the average cosine of the transmitted photons can be expressed as:  As can be seen, the calculation from Eq. ( 5) agrees well with the MC data (within −1 to 4%).At a θ a of 1.5°, the fraction of transmitted photons within this angle varies between 4 × 10 −4 to 8 × 10 −4 , i.e., on average, one photon exits within θ < 1.5° out of ~1250 to 2500 transmitted photons.Only a fraction of the incident photons are transmitted.Therefore, as mentioned earlier, the MC techniques can be prohibitively expensive when applied simulate the AD of multiply scattered photons, especially at relatively large OD; and the method developed above can overcome such computational limitation.Lastly, we discuss the simplifications made in the derivation of Q z , and the extension of the method to anisotropic scatterers.As shown in Fig. 2 (Panels (b) and (c)), the two simplifications made in the derivation of Eq. ( 4) results in more slow decay of Q z compared to the MC data.Such discrepancy can be addressed by either empirically fitting the MC data to obtain Q z , or by removing the simplifications.The former approach is straightforward and our results have shown that an empirically fitted Q z can improve the agreement for calculations shown in Figs. 3 and 4. The latter approach is closely related to the extension to anisotropic scatterers.Equation ( 5) is fundamental and applicable for anisotropic scatterers.However, derivation of P and Q z need to be modified because of the directionality of the scattering for anisotropic scatterers.It has been shown previously that certain scaling laws based on the asymmetric factor of the scatterers can be used to extend isotropic models [11,14] to anisotropic scatterers, and we plan to apply these scaling laws in the derivation of the AD.

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

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

.
the fraction of the transmitted photons that is transmitted within (θ, θ + dθ), and C is a It will be shown later that C = 1 for isotropic scatterers.Intuitively, Eq. (1) divides the slab into a series of infinitely thin slices (each with thickness #151545 -$15.00USD Received 22 Jul 2011; revised 12 Oct 2011; accepted 1 Nov 2011; published 10 Nov 2011 (C) 2011 OSA dz), and calculates the AD of all transmitted photons by averaging the AD of the photons that undergo their last scattering event at each slice.With the pathlength assumed to be exponentially distributed, it can be shown that: can be explained using photon [REMOVED EQ FIELD] shown in Panel (b) of Fig. 1.Suppose the photon undergoes its last scattering event at an angle of θ.In order for it to be transmitted, its free scattering pathlength needs to be longer than z/cosθ as shown, which has a probability of exp(-z/cosθ), because the free scattering pathlength is exponentially distributed.Normalization of this probability leads to Eq. (2).Equation (2) is validated by MC data, with an example set of comparison shown in Panel (b) of Fig. 2.

Fig. 2 .
Fig. 2. Panel (a).Comparison of P against MC data.Panels (b) and (c).Comparison of Qz against MC data under various OD.
in Panels (b) and (c) of Fig.2., Eq. (3) agrees reasonably well with the MC data, and the agreement improves with the optical depth (OD) of the sample.The constant, C in Eq. (1), depends on how the AD is normalized.If the AD is normalized relative to the transmitted photons, i.e.C = 1.One could alternatively normalize Q T (θ) relative to the incident photons such that that are transmitted (i.e., transmittance).
Example results calculated using Eq.(6) are shown in Panel (a), Fig.4, in comparison to MC data.Again, good agreement (within −0.5 to 1.5%) is observed in a wide range of OD and the insensitivity with respect to OD is observed for OD larger than ~5.

Fig. 4 .
Fig. 4. Panel (a).Comparison of the average cosine at various OD.Panel (b).Fraction of transmitted photons within given acceptance angles.Panel (c).Error relative to the MC data.Another quantity of interest to the design of imaging optics is the fraction of photons that can be collected.Such fraction can be calculated by 0