CW laser transilluminance in turbid media with cylindrical inclusions

Optik 122 (2011) 577–581

Contents lists available at ScienceDirect

Optik journal homepage: www.elsevier.de/ijleo

CW laser transilluminance in turbid media with cylindrical inclusions H.O. Di Rocco a,b , D.I. Iriarte a,b , M. Lester a,b , J.A. Pomarico a,b,∗ , H.F. Ranea-Sandoval a,b a b

IFAS - UNCPBA - Pinto 399, B7000GHG Tandil, Argentina CONICET,Argentina

a r t i c l e

i n f o

Article history: Received 15 September 2009 Accepted 26 March 2010

Keywords: Turbid media Diffusion approximation Inclusions Continuous wave

a b s t r a c t The present work analyzes the process of detection of small diffusive inclusions in turbid hosts. Experiments were carried out using a transillumination geometry and continuous wave laser radiation, considering cylindrical inclusions in different environments. A comparison between experimental data, theoretical approaches for the radiative transfer equation in the diffusion approximation, and numerical Monte Carlo simulations, is presented. © 2010 Elsevier GmbH. All rights reserved.

1. Introduction The study of the propagation of light through highly dispersive turbid media, such as biological tissues, is of essential interest in modern biomedical and technical studies, and it is thus a constantly growing field of research, as pointed out in numerous references [1]. It is well known that light in the Near Infrared (NIR) can penetrate several centimeters into tissues, due to the comparatively low absorption at these wavelengths. This is the main reason why NIR light can be used to explore tissues to detect inclusions, inhomogeneities, or alien bodies. The presence of these inhomogeneities in turbid media can, generally speaking, be detected either by time resolved, light modulation or continuous wave (CW) techniques. Moreover, the fact that optical techniques can be implemented for different wavelengths, allows using them to identify the substances that constitute the inclusions. This has been used, for example, to determine the oxygenation level of the blood in tissues, using the detection capabilities of a double wavelength instrumentation. The oxygen content in the blood can also be used to identify the character of tumors, besides localizing them within the organs [2]. The purpose of this paper is to analyze the detection of small inclusions in turbid hosts, and to discuss the procedures most commonly used to validate the models and the parameters involved. The experimental technique we have used is a CW, NIR-based transilluminance. We compare experimental results with both, theoretical models and Monte Carlo (MC) simulations. Samples consisted of slabs of turbid materials including turbid

∗ Corresponding author. Tel.: +54 2293 439660. E-mail address: [email protected] (J.A. Pomarico). 0030-4026/$ – see front matter © 2010 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2010.04.013

inhomogeneities of cylindrical shape. Theoretical calculations were performed within the Diffusion Approximation (DA) of the Radiative Transfer Equation (RTE). This work is mainly motivated by the fact that, despite some authors consider the CW method – which is by far less expensive than time resolved or light modulation techniques – as a lower resolution technique [3,4], we have found that it is, nevertheless, an interesting tool for the detection of inhomogeneities in diffusive media, as it was shown in several previous publications on the subject [5–8]. Its simplicity does not preclude it from attaining good resolution in the transillumination arrangement presented here. Experiments presented are carried out in the transillumination modality, in which source and detector are facing each other defining the optical axis (taken as the Z axis), and they are scanned across the X axis of a slab of turbid medium containing the desired inhomogeneity. Since CW is used, at each position, x, only the transmitted diffuse intensity, It (x), is registered. The paper is organized as follows: in the next section we present a brief description of the theoretical background involved. In

2. Theoretical considerations In this Section we shall give a brief description of the multiple scattering theory, emphasizing what is more relevant to the comparison to the experimental results. As it is well known, turbid media are described by three optical parameters, namely the refractive index, n, the absorption coefficient, a , and the reduced scattering coefficient, s . In a previous publication, we described a CW method for obtaining sufficiently precise values of both, the absorption coefficient and the reduced scattering coefficient, which are required to evaluate the light dis-

578

H.O. Di Rocco et al. / Optik 122 (2011) 577–581

tribution and/or propagation inside the turbid medium [9]. The refractive index needs to be determined from independent measurements. Other coefficient, called the diffusion coefficient, D, is sometimes used; it is related to s through D = 1/(3s + aa ) ∼ = 1/(3s ) when s >> a , which is the case of interest in biological tissues.

Km (x):





∞

scatt (rs , rd ) =

dp

B0 (p) K0



2  p2 + out d

0

× cos (pyd ) cos (m) Km





∞

+2

m=1 2  p2 + out d



.

(4)

2.1. Diffusion theory The propagation of light in turbid media is best described by the Radiative Transfer Equation (RTE), which is the expression of the balance of energy inside a volume element of the scattering medium. An analytical model widely used to describe this phenomenon is the Diffusion Approximation (DA) to the RTE, that proved to be precise even in fairly complex systems but much simpler to implement. Diffusion theory describes light transport accurately only in the strong scattering regime, that is, where the reduced scattering coefficient is much larger than the absorption coefficient, s >> a . This is the situation occurring in biological tissues, and therefore the DA is suitable to describe the physical facts involved in these media. In an infinite, homogenous media, where both of the D and a are constants, the differential equation for a continuous source takes the form of a Helmholtz-like equation

∇ 2  (r) −

a S (r)  (r) = − D D

(1)

where S(r) is the source term; we assume a point source: S (r) = ı (r-rs ). Because in our case the inhomogeneities are cylinders, we need to expand the solution of Eq. (1) in cylindrical coordinates. As we will work in the transillumination modality, we arrange the source and the detector on the z axis, with the restriction ys = yd = 0; moreover, since source and detector face each other, it results that the angles in cylindrical coordinates satisfy: (s − d ) =  (see for example Ref. [10]). For a point source the solution to Eq. (1) is





inc (rs , rd ) = S



exp − rs -rd 





4 rs -rd 



,



p2 + 2 < ) and

p2 + 2 > , being < (> ) the smaller (larger) of s and d .

The propercombination satisfying   the boundary  conditions is the product Im



p2 + 2 < Km

S 0 (rs , rd ) = × 22 Dout K0







∞



p2 + 2 >



dp I0

and then

p2 + 2 <



0

p2 + 2 >





p2 + 2 < Km

∞

+2

cos (pyd ) cos (m) Im

m=1



p2 + 2 >

 ,

(3)

that must be compared with the expression presented by Walker et al [11]. Out of the cylinder, the flux is given by out (rs , rd ) = 0 (rs , rd ) + scatt (rs , rd ), where, for the fulfillment of the condition scatt (rs , rd ) → 0 when rd → ∞ , we must use the functions

∞

in (rs , rd ) =

dp



C0 (p) I0



p2

2 + in d

0 ∞

+2



 2 p2 + in d

cos (pyd ) cos (m) Im

.

m=1

(5) In order to find the coefficients Bm (p) and Cm (p) , the following conditions are used at the boundary between bulk and inclusion [12]:

 n 2 in [0 + scatt ]r=a = [in ]r=a nout    Dout

∂ ∂ 0 + scatt ∂r ∂r

= Din r=a

∂ in ∂r



;

(6)

r=a

The first of that is approximated but, as it is explained in [12], is valid for 1.3 < nout < 1.5, which is the interesting case for this work. The expressions for the coefficients Bm (p) and Cm (p) are long and privy of an evident sense, and therefore they are not written here. In Eq. (6) the lower index “in” indicates the inclusion whereas the lower index “out” indicates the bulk. The normal flux is Jn = J · nˆ = J + (r), where nˆ is the outward unit normal vector, i.e., non-diffusive media is assumed at the interface z = zout . Thus, at the exit boundary of the slab [12,13]:



(r)|S = −˛

functions are the modified Bessel functions Im ( Km



(2)

being  = a /D. The incidence of a spherical wave on a cylinder indicates that we must find the representation of inc (rs , rd ) in cylindrical coordinates; that work was presented by Walker et al. [11], but, to our judgment, it contains some inaccuracies. It can be calculated by solving the Eq. (1) for r = / r s and the proper



Inside the cylinder the solution is given by an expression similar to the former, by now using the Im (x) , because the flux must be finite on the cylinder axis:

∂(r)  , ∂nˆ 

(7)

S

where ˛ takes into account the refractive index mismatch between bulk and air. We have also performed a perturbative analysis of the scattering following the lines of references [14–17], but the results were not accurate enough to compare with our experimental outcomes, since the values for the s and a we should use to fit the experimental ones were greater than those suitable for this approximation to be valid. 3. Experiment A sketch of the experimental setup for a transillumination arrangement is shown in Figs. 1 and 2 gives some details about the geometry of the sample. The samples with the inclusions were mounted on a computercontrolled single-axis Linear Translation Stage (LTS), capable of 1 m minimum step and that can travel up to 60 mm. As a source we used a pigtailed, power-stabilized, temperature controlled, NIR CW laser diode that emits at s = 830 nm with a nominal power of 12 mW at the fiber output. It was placed with the fiber output end in contact with the entrance face of the sample (at plane z = 0), being the illuminated area less than 0.25 mm2 . The recording device was a 768 × 512 pixels, 14 bits CCD camera manufactured by DTA Instruments, Italy, model CHROMA C3 . A Video Zoom Microscope (VZM) images a region of the opposite face of the sample (at plane z = S), which was approximately 1 mm2 . The optical axis is defined

H.O. Di Rocco et al. / Optik 122 (2011) 577–581

Fig. 1. Sketch of the top view of the geometry for the transillumination arrangement. A pigtailed CW Laser Diode (LD) ( = 830 nm) is placed with the output of the Optical Fiber (OF) in contact with the entrance face of the sample (plane z = 0). The sample (host of thickness S containing the inclusion) is placed on a linear translation stage (LTS) and can be moved in the X–Z plane along the X direction. A CCD camera with a Video Zoom Microscope (VZM) images an area, of about 1mm2 on the Optical Axis (OA) at the opposite face of the slab (plane z = S). The mean intensity at this spot is measured by the image processing system.

as the line joining the tip of the Optical Fiber (OF) and the center of the camera lens, and remained fixed during the experiments; it coincides with the Z axis as shown in Figs. 1 and 2. For the transillumination experiments, samples were moved in small steps along the x direction, perpendicular to the optical axis. The actual size of the steps was conveniently decided for each experiment. Once the point of observation was chosen, the transmitted mean intensity, It , at the image spot was recorded by the CCD using its factory-installed software. With the measurements at each step of the scanning, a transillumination intensity profile of a particular sample, It (x), was obtained. For better comparison between model outcomes and experiments, each experimental profile (as well as the corresponding MC simulation and theoretical result) was normalized to the intensity registered at a reference position far away from the inclusion, termed IRef in Figs. 3–5. We have studied both, solid and liquid hosts. They were slabs of thickness S, made of media of which optical parameters had been previously measured by means of a method we devised in reference [9]. The lateral extent of the slab, along X and Y axis (shown as H and W in Fig. 2) were considered to be large enough as to avoid

Fig. 2. Perspective view of the sample. The cylindrical inclusion is placed parallel to the Y axis inside the test sample at a distance zinh from the face in contact with the laser optical fiber. Light enters the slab at plane z = 0 and is collected at plane z = S.

579

Fig. 3. Dependence of the MC simulations on the set of selected optical parameters. The values of g given are the same for both, host and inclusion. The refractive index of the host is n1 , and n2 is that of the inclusion. The values considered for s are taken such that s = s (1 − g) = 10 cm−1 for the host and 18 cm−1 for the inclusion. Five cases are shown: (A) g = 0.6, n1 = n2 = 1.33, (B) g = 0.7, n1 = n2 = 1.33, (C) g = 0.8, n1 = n2 = 1.33, (D) g = 0.8, n1 = 1.4, n2 = 1.5, and (E) g = 0.85, n1 = n2 = 1.33. Please notice that, due to high computation time, MC was run only for positive values of the position, and assuming symmetry, transilluminance values were mirrored about the origin to improve visualisation.

the possibility of detecting photons which have reached these far points. Liquid samples were prepared by mixing whole milk, distilled water and diluted India-ink (see Section 5 for more details). The solid inclusion considered was a cylinder of Grilon® , while in the case of solid hosts, the inclusion was a cylindrical cavity drilled in the host and filled with a turbid solution of known optical parameters. 4. Monte Carlo simulations 4.1. General considerations Monte Carlo simulations were made for both, liquid and polymer slabs, with inhomogeneities in the configurations that best emulate the experimental conditions. The values of the optical parameters used in the simulations were obtained by the method described in reference [9]. Geometry is the same as shown in Fig. 1.

Fig. 4. Example of a solid inclusion in a turbid liquid: Grilon® (inh = 0.7 cm) in milky solution of optical parameters a = 0.08 cm−1 , s = 9.8 cm−1 , filling a tank of S = 3 cm. The experimental profile (solid circles) was taken with the inclusion at zinh = S / 2 = 1.50 cm. The result of this experiment is compared to the scattering theory results (line) and to the Monte Carlo simulation outcome (open circles). As in Fig. 3 MC was run only for positive values of the position, whereas experimental values were taken along the whole profile.

580

H.O. Di Rocco et al. / Optik 122 (2011) 577–581

At each new position of the photon it is decided if absorption takes place or if scattering process goes on by means of R4 . Absorption occurs if a (11) R4 < t and the photon is killed; otherwise, scattering proceeds. Fresnel reflections were taken into account at the input and output faces of the slab. Additionally, photons that reach lateral positions greater than 4S are killed. Fresnel conditions between the bulk and inhomogeneity are taken into account as indicated in Eq. (6). Complementary details of MC simulations for photon propagation in turbid media can be found in [18,20]. 4.2. About the dependence of MC outcomes on the parameters s , g and n Fig. 5. Liquid inclusion in turbid solid environment: a hole (inh = 0.5 cm) was drilled in a Grilon® slab of thickness S = 1.75 cm, and then filled with milky solution of optical parameters a = 0.03 cm−1 , s = 9.2 cm−1 . The inclusion was at zinh = S / 2 = 0.87 cm. Again, MC was run only for positive values of the position, whereas experimental values were taken along the whole profile.

Simulations were carried out by launching at the entrance face of the slab of turbid medium of thickness S, an infinitely narrow photon beam consisting of N0 photons incident normal to the slab. Incidence direction was assumed to be Z, coincident with the optical axis, and X and Y were the directions perpendicular to the incident beam. Photons were launched one at a time starting at positions (x, 0, 0), being x scanned from x = − (a / 2) to x = (a / 2), as in the experiments. For the transillumination configuration, transmitted photons were collected at the opposite face of the slab, in an area centered at coordinates (x, 0, S). The photons collected were those emerging from this area, within the numerical aperture of the optical system. Cylindrical inclusions were located inside the slab at coordinates (0, 0, zinh ) with its generatrix parallel to the faces of the slab in the Y direction. The basic flow diagram adopted in this work for the MC code follows the original idea from Wang et al. [18], and we shall give here only the general guidelines that describe the code. The propagation of each photon is determined by a set of random numbers {Rj , j = 1, 2, 3, 4} uniformly distributed between 0 and 1, the optical properties of the medium, s , a , and by the anisotropy factor g. Between successive scattering events a propagation step of length l given by: l=

− ln(R1 ) t

(8)

is assumed, being t = s + a . The new direction of propagation is determined by two additional random numbers, R2 and R3 , which produce the new propagation angles. The new azimuth angle, ␸, is considered to be uniformly distributed between 0 and 2. Thus: ϕ = 2R2 .

(9)

The deflection angle has to be sampled according to certain probability distribution for its cosine. The usual choice is the wellknown function first proposed by Henyey and Greenstein for stellar atmospheres [19]. In accordance with this, we take

cos =

⎧  ⎪ ⎨ 1 ⎪ ⎩

2g

 2

1+g −

1 − g2 1 − g + 2gR3

2R3 − 1





2  if

g= / 0

if

g = 0,

(10)

being g = cos which, for the case of biological tissue has a value between g = 0.7 and g = 0.9.

It is interesting to note, that while in the theoretical approach the reduced scattering coefficient, s = s (1 − g), is what matters, in MC simulations both, s and the anisotropy factor g are independent. That means that for MC it is not sufficient that the product s (1 − g) equals s , as is the case for the theory. The best values for s and g must be found independently from each other. This conclusion is supported by the results shown in Fig. 3, which corresponds to several MC runs, for different values of g (and thus for s ) while keeping the same values for s of host and inclusion. In particular, in the example of Fig. 3 we considered a typical case for tissue like phantoms, namely s = 10 cm−1 for the host and s = 18 cm−1 for the inclusion. For all simulations the thickness of the host was 3 cm and the inclusion was a long cylinder of diameter 0.7 cm located at the center of the host. The absorption coefficient was kept constant for both, host and inclusion, since it can be taken into account after the simulations, by multiplying the outcome by the factor exp (− a L), where L is the average total photon path. Because of the very long simulation time needed for MC calculations, we run only for 11 positions, each considering 108 photons. Taking into account the results of this figure, which shows a high sensitivity of the transilluminance profile on the set of parameters chosen, and the extremely long calculation time for each MC simulation, the selection of the right set is far from trivial if no other hint about the proper optical parameters is available. This shows the importance of having the solution for some simple theoretical cases, i.e. a very long cylinder in a homogeneous host, as is the case of this work, in order to constrain the range of the parameters to be used in the simulations. Moreover, the influence of the refractive index needs to be considered. In fact, starting from the fundamental solution for infinite media (see Eq. (2) in [6]), it is known that both, a and s are linked to n through the speed of light in the medium. This adds a further degree of freedom for the MC simulations, because a is divided by n while s is multiplied by this factor. As an example, we show in Fig. 3 two cases, both with Host = 50 cm−1 , Inclusion = 90 cm−1 , s s and g = 0.8 (both) but two sets of ns. The difference is striking, since one of them results in a single peaked profile, the other shows a central dip (see Fig. 3 caption). 5. Results and discussion We have investigated two cases, namely solid inclusions in liquid hosts and liquid inclusions in solid (plastic) hosts. All profiles shown correspond to transilluminance experiments. For the first case, we used small tanks made of thin glass walls to contain the liquid host. They were parallelepipeds of 30 cm × 25 cm × 3 cm filled with a solution of whole milk in distilled water, using India-ink (Rotring® black) as an absorber. The optical properties of

H.O. Di Rocco et al. / Optik 122 (2011) 577–581

the solution were measured and calibrated using the fitting procedure described in detail in [9]. A solution in a mixture (3 : 1 : 0.04) of water:milk:India-ink, diluted (1 : 1000) gives average values of s ∼ = 9.8 cm−1 and a ∼ = 0.08 cm−1 . In the following, we shall refer to this mixture (3 : 1 : 0.04) as the milky solution. Varying the concentration of ink and milk allows changing the value of a and s , respectively. Actual optical parameters were determined for every new mixed solution. The cylindrical inclusions were slim bars of diameters much smaller than its length made of a extruded polymer, namely Grilon® . We used a positioning device to properly place the inclusions inside the tanks with respect to the light entrance surface. This device has an accuracy of 1 m or better. Because of their geometry, the optical parameters for the polymer was determined from slices (Ref. [9]) of this material, but not from the actual inclusion. This is an important point, since machining of the bar can produce alterations (i.e. through temperature variations and dirt) in the optical properties. Moreover, local variations in composition of the bulk material can affect measurements, but as given in [9] good average values for Grilon are s = 4.2 cm−1 , a = 0.0016 cm−1 . Fig. 4 shows the results for a solid inclusion made of a Grilon® stick of diameter inh = 0.60 cm immersed in a 3.0 cm thick tank filled with milky solution. The inclusion was at a depth zinh = 1.5 cm, that is at the middle point of the tank. In this case both, theoretical and MC calculations, are fairly in accordance with experimental values for the above set of optical parameters. In the curves shown in Fig. 5, we present the reciprocal case of a hole of diameter inh = 0.50 cm drilled in a Grilon® slab of thickness S = 1.75 cm, and filled with other milky solution, which measured optical parameters were a = 0.03 cm−1 , s = 9.2 cm−1 and n = 1.4. Lateral dimensions of the host were: W = 5 cm and H = 4 cm (refer to Fig. 2 for geometry). The inclusion was at the middle point in the sample, that is zinh = 1.5 cm. Whereas the relative experimental minimum is ≈0.85, the theoretical value is about 0.86 and the corresponding MC value for it is near 0.83. 6. Conclusions We have made studies of CW transilluminance in turbid media in order to detect the inhomogeneities, and we have compared the outcomes of the experiments with MC simulations and Scattering Theory in the DA regime as validation procedures. The optical properties of the material used were taken as known from previous measurements. Unlike theory, that depends only on the absorption and reduced scattering coefficients, we have shown that MC simulations clearly depend on s and g separately; therefore several MC runs are needed in order to compare outcomes to theory. This is a drawback for MC simulations, and it is to be taken into account, since using the analytical model reduces the calculation time (which in case of MC simulations can last weeks) to less than a minute. As an example we have performed a couple of experiments to detect these inhomogeneities in different circumstances and environments. Experimental results compare qualitatively well with both, the Scattering theory and the MC simulations regarding the transillumination intensity modulation between the values at the inhomogeneity and the value of the background transillumination intensity. Theoretical approaches are only available for a few cases of high symmetry, being thus, not useful in general. It may seem that numerical simulations are the answer to this point, but this work shows that the set of possible combinations of optical

581

parameters is so high, that calculation time may become prohibitive. The proposed theoretical model for cylinders as well as the MC simulations fit fairly well the experimental data concerning the modulation height; however both, MC and theory, slightly overestimate the FWHM in all cases (Figs. 4 and 5). In conclusion, we have shown that the CW laser transillumination procedure has enough sensitivity to give reliable results in detecting inhomogeneities immersed in turbid media. These results validate the proposed theoretical model, which is useful due to the almost negligible time needed for computation. It can thus be used to find the most suitable set of optical parameters for a given combination host – inclusion by comparison of it outcomes with experimental data. Acknowledgements Financial support from CONICET under grant no. PIP6431 and from ANCYPT grants PICT039598 and PICT38058/05 are gratefully acknowledged. Dr. Iriarte acknowledges the Third World Academy of Science (TWAS) for equipment donation. References [1] T. Vo-Dinh (Ed.), Biomedical Photonics Handbook, CRC Press, Boca Raton, 2003, see specially Chapters 2, 3,4,21 and 22 and References therein. [2] F. Crespi, A. Bandera, M. Donini, C. Heidbreder, L. Rovati, Non-invasive in vivo infrared laser spectroscopy to analyze endogenous oxy. haemoglobin, deoxyhaemoglobin and blood volume in the rat CNS, J. Neurosci. Methods 145 (2005) 11–22. [3] G. Mitic, J. Kolzer, J. Otto, E. Plies, G. Solkner, W. Zinth, Time-gated transillumination of biological tissues and tissue-like phantoms, Appl. Opt. 33 (1994) 6699–6710. [4] K. Michielsen, H. De Raedt, N. García, Time-gated transillumination and reflection by biological tissues and tissue-like phantoms: simulations versus experiment, J. Opt. Soc. Am. A 14 (1997) 1867–1871. [5] H.O. Di Rocco, D.I. Iriarte, J.A. Pomarico, H.F. Ranea-Sandoval, Object detection in solid Phantoms by Speckle Contrast Measurement, J. Hologr. Speckle 2 (2005) 90–95. [6] H.O. Di Rocco, J.A. Pomarico, H.F. Ranea-Sandoval, Light transmittance in turbid media slabs used to model biological tissues; scaling laws, Optik – Int. J. Light Electron Opt. (2008), doi:10.1016/j.ijleo.2008.07.028. [7] D. Contini, H. Liszka, A. Sassaroli, G. Zaccanti, Imaging of highly turbid media by absorption method, Appl. Opt. 35 (1996) 2315–2324. [8] A.M. Siegel, J.J.A. Marota, D.A. Boas, Design and evaluation of a continuous-wave diffuse optical tomography system, Opt. Express. 4 (1999) 287–298. [9] H.O. Di Rocco, D.I. Iriarte, J.A. Pomarico, H.F. Ranea-Sandoval, R. Macdonald, J. Voigt, Determination of the optical parameters of slices of turbid media by diffuse CW Laser Light scattering profilometry, JQSRT 105 (2007) 68–83. [10] J. Jackson, Classical Electrodynamics, second ed., John Wiley & Sons, 1975. [11] S.A. Walker, D.A. Boas, E. Gratton, Photon density waves scattered from cylindrical inhomogeneities: theory and experiments, Appl. Opt. 37 (1998) 1935–1944. [12] Jorge Ripoll Lorenzo, Difusión de luz en medios turbios con aplicación biomédica, Ph.D. Thesis, Universidad Autónoma de Madrid, Spain, 2000. [13] R. Aronson, Boundary conditions for diffusion of light, J. Opt. Soc. Am. A 12 (1995) 2532–2539. [14] J.C. Paasschens, G.W.T. Hooft, Influence of boundaries on the imaging of object in turbid media, J. Opt. Soc. Am. A 15 (1998) 1797–1812. [15] M.R. Ostermeyer, S.L. Jacques, Perturbation theory for the use light transport in complex biological tissues, J. Opt. Soc. Am. A 14 (1997) 255–261. [16] S. Feng, F.-A. Zeng, B. Chance, Photon migration in the presence of a single defect: a perturbation analysis, Appl. Opt. 34 (1995) 3826–3837. [17] M.R. Ostermeyer, S.L. Jacques, Perturbation theory for optical diffusion theory: a general approach for absorbing and scattering objects in tissues, in: B. Chance, R.R. Alfano (Eds.), Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media: Theory, Human Studies, and Instrumentation, SPIE 2389, 1995, pp. 98–102. [18] L. Wang, S.L. Jacques, L. Zheng, MCML - Monte Carlo modeling of light transport in multi-layered tissues, Comput. Methods Programs 47 (1995) 131–146. [19] L.G. Henyey, J.L. Greenstein, Diffuse radiation in the Galaxy, Astrophys. J. 93 (1941) 70–83. [20] F. Martelli, D. Contini, A. Taddeucci, G. Zaccanti, Photon migration through a turbid slab described by a model based on diffusion approximation. II. Comparison with Monte Carlo results, Appl. Opt. 36 (1997) 4600–4611.