OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 18, Iss. 9 — Apr. 26, 2010
  • pp: 9456–9473
« Show journal navigation

Light diffusion in a turbid cylinder. I. Homogeneous case

André Liemert and Alwin Kienle  »View Author Affiliations


Optics Express, Vol. 18, Issue 9, pp. 9456-9473 (2010)
http://dx.doi.org/10.1364/OE.18.009456


View Full Text Article

Acrobat PDF (1101 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

This paper is the first of two dealing with light diffusion in a turbid cylinder. The diffusion equation was solved for a homogeneous finite cylinder that is illuminated at an arbitrary location. Three solutions were derived for an incident δ-light source in the steady-state, frequency, and time domains, respectively, applying different integral transformations. The performance of these solutions was compared with respect to accuracy and speed. Excellent agreement between the solutions, of which some are very fast (< 10ms), was found. Six of the nine solutions were extended to a circular flat beam which is incident onto the top side. Furthermore, the validity of the solutions was tested against Monte Carlo simulations.

© 2010 Optical Society of America

1. Introduction

Light propagation in turbid media, e.g. biological tissue, is described with different equations depending on the considered applications and the required accuracy. On the microscopic scale Maxwell equations are usually used [1

1. F. Voit, J. Schäfer, and A. Kienle, “Light Scattering by Multiple Spheres: Comparison between Maxwell Theory and Radiative-Transfer-Theory Calculations,” Opt. Lett. 34, 2593–2595 (2009). [CrossRef] [PubMed]

], whereas on the mesoscopic scale transport equation and on the macroscopic scale diffusion equation are normally applied [2

2. A. Ishimaru, Wave Propagation and Scattering in Random Media, (Academic Press, New York, 1978).

].

The advantage of using the diffusion equation is that analytical solutions can be obtained for a variety of geometries. The disadvantage is that due to its approximations its applicability is restricted and, therefore, has to be checked against more precise theories, e.g. transport theory. Analytical solutions of the diffusion equation were found, e.g. for infinite, semi-infinite, and slab geometries [3–5

3. M. S. Patterson, B. Chance, and B. C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

], further for a sphere [6

6. S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

, 7

7. B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

], for a parallel-epiped [8

8. A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

], and for a cylinder [6

6. S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

, 7

7. B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

, 9

9. A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]

].

Besides these homogeneous geometries solutions of the diffusion equation were also derived for inhomogeneous geometries. Especially, layered geometries aroused interest because many different parts of the body exhibit layered structures. In this study we present solutions of the diffusion equation for a homogeneous cylinder, whereas in the accompanying paper solutions for a cylinder having an arbitrary number of layers are given [10

10. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

].

In the field of biomedical optics the applicability of the derived solutions is twofold. First, in phantom experiments often cylindrical geometries are used, so that it is obviously required to apply the corresponding solutions. Second, for in-vivo experiments the cylindrical geometry can be used to describe approximately parts of the human body or of animals, e.g. the arm, the leg, or the finger. Of course, these examples are not exactly described by a cylinder but, definitely, the assumption of a cylindrical geometry is more exact than the mostly used semi-infinite geometry. We note that the derived equations can be used to obtain the corresponding solutions of the heat conduction equation, as the diffusion equation includes the heat equation as a special case. In this article we present three different solutions of the diffusion equation for a homogeneous finite cylinder in the steady-state, frequency, and time domains, respectively. For these nine solutions the assumed pencil beam is incident onto the cylinder barrel or onto the cylinder top (or bottom). The solutions were compared among each other for the different domains. We discuss how the solutions given in the literature so far are related to ours. In addition, we derived the solutions of the diffusion equations for a circular flat beam with a finite diameter incident onto the cylinder top (or bottom). Finally, the solutions were compared to results obtained by the Monte Carlo method, which is a numerical solution of the transport equation. In general, the derived solutions can be easily implemented using standard programming languages. In addition, we provide a freely available executable code on the internet [11].

2. Theory

2.1. Diffusion Theory

In this section we present different solutions of the diffusion equation for a finite homogeneous cylinder. In the first part the finite Hankel and cosine transforms are used to obtain an ordinary differential equation for the z-coordinate, see Fig. 1. The boundary conditions in z-direction are then implemented applying two different approaches in the frequency domain: use of an infinite series of point sources (version A) and use of hyperbolic sine functions (version B), and applying two approaches in the time domain: use of the finite sine transform (version A) and again use of an infinite series of point sources (version B). In addition, for these six cases we give also the solutions for a flat beam incident onto the cylinder top. In the second part we use the finite sine transform in z-direction and then solve the resulting modified Bessel differential equation both in the frequency and time domains. We note that the results in the steady-state domain can be easily obtained from those derived in the frequency domain by setting the modulation frequency ω = 0.

In literature Arridge et al. presented a two-dimensional solution of the diffusion equation for a cylinder in the frequency domain solving the modified Bessel differential equation [6

6. S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

]. In addition, they gave the solution of the diffusion equation for a finite cylinder in the frequency and time domains, but only for the incident plane. Sassaroli et al. extended this solution to locations outside the incident plane and for the extrapolated boundary condition [9

9. A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]

]. Pogue and Patterson indicated the solution of the diffusion equation for a finite cylinder in the time domain, which corresponds to version B in this study, but they did not present a derivation of the solution [7

7. B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

].

Figure 1 shows the cylindrical geometry used in the simulations. The center of the coordinate system is the middle of the cylinder top. The radius and height of the cylinder are denoted as a and lz, respectively. The incident light beam position can be chosen arbitrarily, e.g. it can be incident onto the top (or bottom) and onto the barrel. It is assumed that the incident light beam can be represented by a point source in the turbid medium in direction of the beam at a distance of l/(μs + μa) from the location of (perpendicular) incidence at the cylinder barrel, see Fig. 1, where μa and μs are the absorption and reduced scattering coefficients of the homogeneous cylinder. The position of the point source is given in cylindrical coordinates r0 = (ρ 0, ϕ 0,z 0). In case of an oblique incident beam the Snel’s law has to be applied to determine the position of the point source. Extrapolated boundary conditions are used to describe the influence of the border between the turbid medium and the surrounding, characterized by the refractive indices, n 1 and n 0, respectively. As in earlier papers we use the fluence rate and the flux term for calculation of the reflectance in the steady-state domain [12

12. A. Kienle and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]

], whereas for the other solutions only the flux term (Fick’s law) is applied. The rationale for this is the better agreement to Monte Carlo simulations when using the fluence and flux terms for calculation of the spatially resolved reflectance [12

12. A. Kienle and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]

].

2.1.1. Solutions derived via the finite Hankel transform

We start with the diffusion equation in the frequency domain

ΔΦ(r,ω)(μaD+iωDc)Φ(r,ω)=1Dδ(rr0),
(1)
Fig. 1. Scheme of the cylinder used in the calculations. Exemplarily, the beam is incident perpendicular onto the cylinder top at r⃗ = (ρ 0, ϕ 0,0). The isotropic point source is located at r⃗ = (ρ 0, ϕ, z 0)

where Φ, D = 1/(3μs), c, and ω denote the fluence rate, the diffusion coefficient, the speed of light in the turbid medium and the angular frequency of the intensity modulated light, respectively. Equation (1) can be rewritten by using cylindrical coordinates as follows

2ρ2Φ+1ρρΦ+1ρ22ϕ2Φ+2z2Φ(μaD+iωDc)Φ=1δ(ρρ0)δ(ϕϕ0)δ(zz0),
(2)

where we use the abbreviation Φ = Φ(r⃗, ω). First, a cosine transform of the form

Φ(ϕ,m)=02πΦ(ϕ)cos(m(ϕϕ)),m=0,1,2,
(3)

is applied to Eq. (2) in order to eliminate the ϕ derivative. Thus, we obtain

2ρ2Φ+1ρρΦm2ρ2Φ+2z2Φ(μaD+iωDc)Φ=cos(m(ϕϕ0))δ(ρρ0)δ(zz0).
(4)

Equation (4) can be reduced to an ordinary differential equation by using the finite Hankel transform of mth order

Φ(sn)=0aρΦ(ρ)Jm(snρ)dρ,
(5)

where sn are the positive roots of the Bessel functions of first kind and mth order divided by a′, i.e. Jm(asn) = 0. The roots of the Bessel functions are determined using Halley’s method and applying the start values described by Abramowitz and Stegun [17

17. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, (Dover Publications, New York, 1971).

]. The upper limit of the integral a′ results from the extrapolated boundary condition Φ(ρ = a′) = 0, where we use a′ = a + zb. The extrapolation length zb is calculated with

zb=1+Reff1Reff2D,
(6)

where Reff represents the fraction of photons that is internally diffusely reflected at the boundary to the non-scattering surrounding. Reff was calculated by solving the integrals presented by Haskell et al [13

13. R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary Conditions for the Diffusion Equation in Radiative Transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

]. For the finite Hankel transform the relation

0aρ(2ρ2Φ+1ρρΦm2ρ2Φ)Jm(snρ)dρ=sn2Φ(sn)+asnJm+1(asn)Φ(ρ=a)
(7)

is valid [14

14. E. Meissel and G.B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics, (Bibliobazaar, 2008).

]. By applying the finite Hankel transform to Eq. (4) and by considering the boundary condition at ρ = a′, we obtain the following ordinary differential equation where Φ = Φ(sn, ϕ, m, z, ω)

2z2Φ(sn2+μaD+iωDc)Φ=1DJm(snρ0)cos(m(ϕϕ0))δ(zz0).
(8)

We define G(sn, z, ω as the solution of

2z2G(sn,z,ω)α2G(sn,z,ω)=1Dδ(zz0),
(9)

where α 2 = s 2 n + μa/D + /(Dc). This equation is similar to the equation we derived solving the diffusion equation with lateral infinitely large extensions [15

15. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Nonin-vasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

]. In addition, we note that for the case of an N-layered cylinder the corresponding N-layered equation has to be solved, see accompanying paper [10

10. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

]. Equation (8) is solved by using the solution of Eq. (9) given below in version A and B. Thus, Φ can be obtained with

Φ(sn,ϕ,m,z,ω)=G(sn,z,ω)Jm(snρ0)cos(m(ϕϕ0)).
(10)

We note that Eqs. 3 and 5 together give the following transform

Φ(sn,ϕ,m)=02π0aρΦ(ρ,ϕ)Jm(snρ)cos(m(ϕϕ))dρdϕ,
(11)

with which Eq. (8) was obtained from Eq. (2). Next, an inversion formula of Eq. (11) is obtained by, first, deriving the inversion formula of Eq. (3) and then of Eq. (5). By summing both sides of Eq. (3) over all m ≥ 0, the following equation is obtained

m=0Φ(ϕ,m)=02πΦ(ϕ)m=0cos(m(ϕϕ))dϕ.
(12)

By applying the relation

m=0cos(m(ϕϕ))=12+πδ(ϕϕ)
(13)

we get

02πΦ(ϕ)[12+πδ(ϕϕ)]=Φ(ϕ,0)+πΦ(ϕ)
(14)

and, finally, we obtain

Φ(ϕ)=12πm=0(2δm,0)Φ(ϕ,m).
(15)

An inverse transform for the second transform, the Hankel transform (Eq. (5)), is found by using the completeness relation for the Bessel functions [14

14. E. Meissel and G.B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics, (Bibliobazaar, 2008).

]

n=1ρJm(snρ)Jm(snρ0)Jm+12(asn)=a22δ(ρρ0).
(16)

By applying Eq. (15) and Eq. (16) to Eq. (11), we find the inversion formula as follows

Φ(ρ,ϕ)=1πa2m=(2δm,0)n=1Φ(sn,ϕ,m)Jm(snρ)Jm+12(asn).
(17)

By using Eq. (10), the solution of Eq. (1) in the frequency domain is given by

Φ(r,ω)=1πa2m=0(2δm,0)cos()n=1G(sn,z,ω)Jm(snρ0)Jm(snρ)Jm+12(asn),
(18)

where φ = ϕ - ϕ 0. In order to simplify Eq. (18) the following recurrence relation is applied

Jm+1(asn)=2masnJm(asn)Jm1(asn)=Jm1(asn).
(19)

Thus, we obtain for the general case of a pencil beam incident onto an arbitrary location of the cylinder

Φ(r,ω)=1πa2m=cos()n=1G(sn,z,ω)Jm(snρ0)Jm(snρ)Jm+12(asn).
(20)

For the special case where the pencil beam is incident onto the center of the cylinder top, r⃗ = (0,0,0), Eq. (20) becomes

Φ(r,ω)=1πa2n=1G(sn,z,ω)J0(snρ)J12(asn).
(21)

As stated above, the light transmitted from the cylinder is calculated using Fick’s law. The results for the reflectance and transmittance from the top and bottom of the cylinder are easily obtained and, thus, are not given. For the case of the transmittance from the cylinder barrel the following equation

T(ϕ,z,ω)=DΦ(ρ,ϕ,z,ω)/ρ|ρ=a
(22)

is used to obtain

T(ϕ,z,ω)=Dπa2m=cos()
×n=1G(sn,z,ω)Jm(snρ0)(snJm+1(asn)maJm(asn))Jm+22(asn).
(23)

For the special case of central incidence we get

T(z,ω)=Dπa2n=1snG(sn,z,ω)J1(asn)J12(asn).
(24)

The results presented above were obtained by assuming an incident δ-source. For the derivation of the solution of a finite flat beam incident onto the center of the cylinder top we use the diffusion equation for the rotationally symmetric case (no ϕ-dependence)

2ρ2Φ+1ρρΦ+2z2Φ(μaD+iωDc)Φ=1DS(r,ω),
(25)

where S(r⃗, ω) is a rotationally symmetric source. For a flat beam we have

S(r,ω)=1πρw2circ(ρρw)δ(zz0),
(26)

where ρw is the radius of the flat beam and circ is defined as

circ(ρρw)={1,ρρw0,ρ>ρw.
(27)

By employing the finite Hankel transform of zero order, Eq. (25) becomes

2z2Φ(μaD+sn2+iωDc)Φ=1DS(sn,z,ω)=1DJ1(snρw)πsnρwδ(zz0).
(28)

By using the solution of Eq. (9), G(sn,z,ω), Eq. (28) is solved with

Φ(sn,z,ω)=J1(snρw)πsnρwG(sn,z,ω).
(29)

The inverse relation for the finite Hankel transform of zero order is

Φ(r,ω)=2a'2n=1Φ(sn,z,ω)J0(snρ)J12(sna).
(30)

Thus, the solution of the diffusion equation for an incident flat beam is given by

Φ(r,ω)=2πa2ρwn=1G(sn,z,ω)J0(snρ)J1(snρw)snJ12(sna).
(31)

We note that Eq. (31) is also used for solving the diffusion equation for N-layered turbid media, see the accompanying paper [10

10. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

]. The final step for obtaining the complete solutions is the derivation of the Green’s function, G, see Eq. (9). The boundary conditions in z-direction for the homogeneous cylinder are given by

Φ(ρ,ϕ,z=zb,ω)=0
Φ(ρ,ϕ,z=lz+zb,ω)=0.
(32)

Below we present two solutions in the frequency and two solutions in the time domains.

Solutions in the frequency domain

Version A

We start with the ordinary differential equation given in Eq. (9). It can be solved with the method of image sources where the exponential functions are summed to satisfy the boundary conditions for the z-direction given in Eq. (32). The solution is

G(sn,z,ω)=12k=exp(αzz1k)exp(αzz2k),
(33)

where

z1k=2klz+4kzb+z0
z2k=2klz+(4k2)zbz0
α=μaD+sn2+iωDc.
(34)

Thus, by inserting Eq. (33) in Eq. (20) we obtain the fluence rate for the pencil beam incident at an arbitrary location onto the cylinder as

Φ(r,ω)=12πDa2m=cos()n=1Jm(snρ0)Jm(snρ)αJm+12(asn)
×k=exp(αzz1k)exp(αzz2k).
(35)

For the case when the pencil beam is incident at the center of the cylinder Eq. (35) is simplified to

Φ(r,ω)=12πDa2n=1J0(snρ)αJ12(asn)k=exp(αzz1k)exp(αzz2k).
(36)

Version B

There is another possibility to solve Eq. (9) for the boundaries given in Eq. (32). By applying the method described in reference [15

15. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Nonin-vasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

] the Green’s function for a one-dimensional layer is given in the interval 0 ≤ z < lz as

G(sn,z,ω)=exp(αzz0)exp(α(z+z0+2zb))2sinh[z0+zb]
×sinh[α(z+zb)]exp(α(lz+2zb))sinh[α(lz+2zb)].
(37)

We note that for the evaluation of Eq. (37) no approximation is needed, in contrast to the solutions presented earlier [15

15. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Nonin-vasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

]. In general, the solution calculated with Eq. (37) is faster than the solution obtained with Eq. (33), if lz is small.

Finally, Eq. (37) is inserted in Eq. (20) to obtain the fluence rate for a pencil beam incident at an arbitrary position, in Eq. (21) for a pencil beam incident in the middle of the cylinder top, or in equation Eq. (31) for an incident flat beam.

Solutions in the time domain

Version A

For the derivation of the solution of the fluence rate in the time domain we start with the ordinary differential equation [Eq. (8)]. The boundary conditions given in Eq. (32) are fulfilled by using a finite sine transform:

Φ(λk)=z1z2Φ(z)sin(λk(zz1))dz,λk=/(z2z1),k=1,2,,
(38)

where z 1 = -zb and z 2 = lz+zb. By applying Eq. (38) and the following relation

z1z22z2Φ(z)sin(λk(zz1))dz=λk[Φ(z1)(1)kΦ(z2)]λk2Φ(λk)
(39)

we obtain from Eq. (8)

(sn2+λk2+μaD+iωDc)Φ=1DJm(snρ0)sin(λk(z0+zb))cos(m(ϕϕ0)).
(40)

Thus, the fluence rate in the frequency domain is given by

Φ(sn,ϕ,m,λk,ω)=cJm(snρ0)sin(λk(z0+zb))cos(m(ϕϕ0))Dc(sn2+λk2)+μac.
(41)

The fluence rate in the time domain is obtained from the solutions in the frequency domain using the following Fourier transform

exp(at)ε(t)exp(iωt)dt=1a+,a>0,
(42)

where ε(t) describes the Heaviside step function. By employing Eq. (42), Eq. (41) becomes

Φ(sn,ϕ,m,λk,t)=cJm(snρ0)sin(λk(z0+zb))cos(m(ϕϕ0))
×exp(μact)exp(Dcλk2t)exp(Dcsn2t),t>0.
(43)

We use the completeness relation for the sine function to obtain an inversion formula of Eq. (38)

Φ(z)=2z2z1k=1Φ(λk)sin(λk(zz1)).
(44)

Finally, the inversion formula given in Eq. (20) is used to obtain the time resolved fluence rate in real space

Φ(r,t)=2πa2(lz+2zb)exp(μact)k=1exp(Dcλk2t)sin(λk(z0+zb))sin(λk(z+zb))
×m=cos()n=1exp(Dcsn2t)Jm(snρ)Jm(snρ)Jm+12(asn).
(45)

For a pencil beam incident onto the center of the cylinder top Eq. (45) is simplified to

Φ(r,t)=2cπa2(lz+2zb)exp(μact)k=1exp(Dcλk2t)sin(λk(z0+zb))sin(λk(z+zb))
×n=1exp(Dcsn2t)J0(snρ)J12(asn).
(46)

Version B

We also present a second possibility to solve Eq. (8). Considering again the boundary conditions given in Eq. (32) the series of point sources already applied for version A of the frequency solutions is used, see Eq. (33), where α is now α=μa/D+sn2+/Dc. By employing the Fourier transform

1πtexp(at)exp(b24t)ε(t)exp(iωt)dt=exp(ba+)a+,a>0
(47)

the solution for the homogeneous cylinder in the time domain is for an arbitrary incident pencil beam

Φ(r,t)=12πa2cDπtexp(μact)k=exp((zz1k)24Dct)exp((zz2k)24Dct)
×m=cos()n=1exp(Dcsn2t)Jm(snρ0)Jm(snρ)Jm+12(asn).
(48)

For the pencil beam incident at the middle of the cylinder top we get

Φ(r,t)=12πa'2cDπtexp(μact)k=exp((zz1k)24Dct)exp((zz2k)24Dct)
n=1exp(Dcsn2t)J0(snρ0)J12(a'sn).
(49)

2.1.2. Solutions derived via modified Bessel differential equation

In this section we derive the solutions of the homogeneous diffusion equation for an incident pencil beam by solving the modified Bessel differential equation. We start in the frequency domain and, then, give the time domain solution.

The modified Bessel differential equation is obtained by neglecting the z-dependency in Eq. (2)

2ρ2Φ+1ρρΦ+1ρ22ϕ2α2Φ=1Dδ(ρρ0)δ(ϕϕ0),
(50)

where α=μa/D+/Dc. This differential equation can be solved considering the radial boundary condition Φ(ρ = a′, ω) = 0. The solution can be derived by using a formula given in reference [6

6. S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

] as

Φ(ρ,φ,ω)=12πD[K0(αρ2+ρ022ρρ0cosφ).
m=cos()Im(αρ0)Km(αa)Im(αρ)Im(αa)].
(51)

Φ(r,ω)=14πDk=exp(μeffr1)r1exp(μeffr2)r21πD(lz+2zb)
×k=1sin(λk(z0+zb))sin(λk(z+zb))m=cos()Im(αρ0)Km(αa)Im(αρ)Im(αa),
(52)

where

r1=ρ2+ρ022ρρ0cosφ+(zz1k)2
r2=ρ2+ρ022ρρ0cosφ+(zz2k)2
μeff=μa/D+/Dc.
(53)

The transmitted light from the cylinder barrel is found by using the ρ derivative for the modified Bessel functions, see Eq. (22). The result is

T(φ,z,ω)=aρ0cosφ4πk=(μeff+1/r1)exp(μeffr1)r12(μeff+1/r2)exp(μeffr2)r22
+1π(lz+2zb)k=1sin(λk(z0+zb))sin(λk(z+zb))
×m=cos()Im(αρ0)Km(αa')[αIm+1(αρ)+m/aIm(αa)]Im(αa').
(54)

For a pencil beam incident onto the middle of the cylinder top the fluence rate is

Φ(r,ω)=14πDk=exp(μeffr1)r1exp(μeffr2)r2
1πD(lz+2zb)k=1sin(λk(z0+zb))sin(λk(z+zb))K0(αa')I0(αρ)I0(αa'),
(55)

and the transmittance is obtained with

T(z,ω)=a4πk=(μeff+1/r1)exp(μeffr1)r12(μeff+1/r2)exp(μeffr2)r22
1π(lz+2zb)k=1αsin(λk(z0+zb))sin(λk(z+zb))K0(αa)I1(αa)I0(αa).
(56)

For the solution in the time domain, the frequency domain solution, Φ(r⃗, ω), is used and split into real and imaginary parts [17

17. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, (Dover Publications, New York, 1971).

]. Then, the inverse Fourier transform is calculated by using the FFT algorithm.

Alternatively, we used the ‘damping’ rule of the Laplace transform for separation of that part in Eq. (52) which contains the z-coordinate [18

18. A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

]. The result is

Φ(r,t)=14πDπDctk=exp((zz1k)24Dct)exp((zz2k)24Dct)
×[12texp(μact)exp(ρ2ρ022ρρ0cosρ4Dct)
12πm=cos()Im(μeffρ0)Km(μeffa)Im(μeffρ)Im(μeffa)exp(iωt)dω],
(57)

and the transmittance is obtained with

T(φ,z,t)=14ππDctk=exp((zz1k)24Dct)exp((zz2k)24Dct)
×[aρ0cosφ4Dct2exp(μact)exp(a2ρ022aρ0cosρ4Dct)
+12πm=cos()Im(μeffρ0)Km(μeffa)Im(μeffa)Im(μeffa)exp(iωt)dω],
(58)

where Im(μeffa) = ∂Im(μeffρ)/∂ρρ=a = μeff I m+1(μeffa)+m/aIm(μeffa).

The integral is solved by applying the inverse discrete Fourier transform. The calculations using Eq. (58) are approximately 2 orders of magnitude faster than those calculated without using the ‘damping’ rule.

2.2. Monte Carlo Simulations

The solutions of the diffusion equation are compared to results obtained by Monte Carlo simulations. The Monte Carlo program we developed is able to handle different geometries. For the simulations presented in this study a cylindrical geometry is employed. The light is incident perpendicular to the cylinder top both as a pencil light beam and as a circular flat beam having a finite radius. The refractive index of the cylinder and the surrounding medium is set to n 0 = n 1 = 1.0. The Henyey-Greenstein phase function is used for calculating the scattering angles using an anisotropy factor of g = 0.8. 107 photons were used in Monte Carlo simulations.

3. Results

In this section we, first, compare the different solutions of the diffusion equation derived in the steady-state, frequency, and time domains among themselves. Then, we compare the derived solutions of the diffusion theory with Monte Carlo simulations.

3.1. Comparison of the different solutions of the diffusion equation

For the first comparison the cylinder is illuminated by a pencil beam. It is incident onto the cylinder barrel perpendicular to the barrel at a depth of z = 8mm, see inset of Fig. 2. The steady-state transmittance is calculated around the cylinder barrel at a depth of z = 10mm. The height and radius of the cylinder are lz = 20mm and a = 5mm, respectively. The reduced scattering coefficient is μs = 0.9mm-1, whereas the absorption coefficient is varied (μa = 0.005,0.01,0.015mm-1). For the calculations shown in Fig. 2 we used Eq. (54) for ω = 0. As expected the transmitted light decreases up to the opposite point relative to the incident beam (at ϕ = π) and the transmittance is smaller for increasing absorption coefficients.

Next, we compare the different solutions of the steady-state diffusion equation for the cylinder shown in Fig. 2. The relative difference of two of the three solutions [(Eqs. (35)–(54)]/[Eq. (35)] is shown in Fig. 3. The relative differences are in the range of 10-10, thus, the agreement is excellent. The relative differences between the two steady-state solutions calculated via the Hankel transform (version A and version B) are even smaller (data not shown).

Subsequently, the solutions in the time domain are compared. Figure 4 shows the time resolved transmittance from a homogeneous cylinder having a radius of a = 9mm and a height of lz = 12mm. The cylinder barrel is illuminated at z 0 = 9mm and the transmittance is detected at z = 9mm on the opposite site of the cylinder (ϕ = π) (black curve) and at the middle of the cylinder bottom (red curve).

Fig. 2. Steady-state radial transmittance around the cylinder barrel at z = 10mm and ρ = a. The cylinder is illuminated perpendicular to the cylinder barrel at z 0 = 8mm, see inset. The optical parameters are μs = 0.9mm-1, μa = 0.005 (black),0.01 (green),0.015-1(red)mm , n 0 = 1.0, and n 1 = 1.4. The geometrical data of the cylinder are a = 5mm and lz = 20mm.
Fig. 3. Relative difference of two steady-state solutions for the cylinder considered in Fig. 2.
Fig. 4. Time resolved transmittance from a cylinder that is illuminated at the cylinder barrel at z 0 = 9mm. The transmittance is detected in z-direction at z = lz, ρ = 0mm (red curve) and in radial direction at z = 9mm, ρ = a (black curve). The optical and geometrical parameters are μs = 0.8mm-1, μa = 0.03mm-1, n 0 = 1.0, n 1 = 1.4, a = 9mm, and lz = 12mm. The dashed curves show the transmittance from a homogeneous cube.

Figure 4 shows the results calculated with Eq. (48) (solid curves). The relative difference between the three solutions in the time domain is smaller than 10-10 (data not shown). For comparison, in Fig. 4 also the time resolved transmittance from a rectangular parallelepiped is shown for both detection positions (dashed curves) [8

8. A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

]. The lengths of the edges of the parallelepiped in x-, y-, and z-directions are 9mm, 9mm, and 12mm, respectively. As expected, the results diverge for longer times when the detected photons have experienced the region where the two geometries (cylinder and parallelepiped) are different.

For frequency domain measurements an intensity modulated beam, which consists of a continuous and a sinusoidally oscillating part, is incident onto the turbid medium. In the following we consider only the oscillating part. Figure 5 shows the results of the oscillating part (ω 0 = 2π · 500MHz) for a pencil beam that is incident onto the cylinder barrel computed with the solution obtained by solving the modified Bessel differential equation. The light transmitted from the barrel at an angle of ϕ = 3π/4 at the same height is shown. The transmittance is computed for different cylinder radii a = 5mm (black curve), a= 6mm (green curve), a = 7mm (red curve).

In order to compare the solution shown in Fig. 5 to the solution obtained with the hyperbolic functions (version B) in the frequency domain, we calculated the transmittance by using the following time harmonic signal

S(r,t)=δ(rr0)cos(ω0t).
(59)

The transmittance is given by

T(r,t)=Re{T(r,ω=ω0)eiω0t}.
(60)

The disagreement between the two solutions is small. Thus, we give the numerical values for the above curve having a radius of a = 6mm:

Fig. 5. The transmittance of a sinusoidally oscillating source incident at z = 12mm, ρ = a and detected at the same height at ϕ = 3π/4. The optical and geometrical parameters are μs = 0.7mm-1, μa = 0.02mm-1, n 0 = 1.0, n 1 = 1.4, a = 5,6,7mm, lz = 22mm, z 0 = 11mm. The modulation frequency is ω 0 = 2π · 500MHz.
T(r,ω0)=104mm2×{2.57749063871.2540621701i2.57749063661.2540621700i.
(61)

The upper complex number is calculated with the solution derived by solving the modified Bessel differential equation and the lower number with the hyperbolic functions. Thus, the relative differences are in the order of 10-10.

3.2. Comparison with Monte Carlo simulations

The solutions of the diffusion equation are, first, compared to Monte Carlo simulations for an incident pencil beam, and, then, for an incident flat beam. Figure 6 and Fig. 7 show the spatially resolved reflectance and transmittance from the cylinder top and bottom, respectively. The radius of the cylinders is a = 14mm, whereas the height of the cylinders is lz = 5mm (blue curves) and lz = 10mm (red curves). The optical properties are μs = 1.3mm-1 and μa = 0.008mm-1.

In general, the results of the Monte Carlo simulations (circles) agree well with those obtained from the diffusion theory (solid curves) both for the spatially resolved reflectance and transmittance. For small distances of the spatially resolved reflectance the differences are larger due to the well-known break down of the diffusion theory. At large distances the influence of the border on the light propagation can be seen causing a stronger decrease in the remitted and transmitted intensities. As expected, the increase of the height of the cylinder results in an increase of the remitted light especially at large distances from the source, whereas the transmittance is decreased at small distances and increased at large distances.

Figure 8 shows the spatially resolved reflectance from a homogeneous cylinder that is illuminated in the middle of the top side by a flat beam having radii of ρw = 0mm (blue curve), ρw = 5mm (green curve), and ρw = 10mm (red curve). We note that the beam with ρw = 0mm corresponds to a pencil beam. The optical properties of the cylinder are μs = 1.2mm-1 and μa = 0.01mm-1. The radius and height are a = 20mm and lz = 5mm, respectively.

Fig. 6. Comparison of the spatially resolved reflectance calculated with the solution of the diffusion equation and with the Monte Carlo method for a turbid cylinder illuminated with a pencil beam in the middle of the cylinder top. The optical and geometrical parameters of the cylinders are μs = 1.3mm-1, μa = 0.008mm-1, a = 14mm, n 1 = n 0 = 1.0, lz = 5 and 10mm.
Fig. 7. Comparison of the spatially resolved transmittance calculated with the solution of the diffusion equation and with Monte Carlo simulations for a turbid cylinder illuminated with a pencil beam in the middle of the cylinder top. The optical and geometrical parameters of the cylinders are μs = 1.3mm-1, μa = 0.008mm-1, a = 14mm, n 1 = n 0 = 1.0, lz = 5 and 10mm.
Fig. 8. Spatially resolved reflectance from a homogeneous cylinder that is illuminated by a flat beam with ρw = 0,5,10mm. The optical and geometrical parameters are μs = 1.2mm-1, μa = 0.01mm-1, a = 20mm, l2 = 5mm, and n 1 = n 0 = 1.0.

In general, Fig. 8 shows good agreement between diffusion theory and Monte Carlo simulations. Interestingly, the break down of the diffusion equation at small distances can only be seen for ρw = 0mm. For large beam diameters the agreement is even good for small distances from the source.

4. Discussion

Solutions of the diffusion equation for a homogeneous cylinder were derived in the steady-state, frequency, and time domains. Comparison of the reflectance and trans-mittance calculated with these solutions showed excellent agreement among themselves in each domain. The speed for calculation of the solutions, however, are different. In the time domain we found for the general case of a beam incident at an arbitrary position that the solution obtained via solving the modified Bessel differential equation is about hundred times faster than the other two solutions. For example, the time needed for obtaining the time resolved transmittance using Eq. (58) is in the order of 10ms on a state-to-the-art computer programmed in Pascal (Delphi).

The derived solutions permit the calculation of the fluence rate for a pencil beam that is incident at an arbitrary position onto the surface of the cylinder. In addition, we derived the solutions for a circular flat beam incident onto the plane surfaces of the turbid cylinder for all three domains.

The solutions of the diffusion equation were compared to Monte Carlo simulations. Good agreement was found as expected from results calculated for turbid media having other geometries. Interestingly, quite a good agreement was obtained for the spatially resolved reflectance at small distances from the center of the incident circular beam having finite beam diameters, although it is expected that diffusion theory breaks down in this case. Probably, the distance dependent errors of the solution of the diffusion equation for an incident pencil beam (at small distances diffusion theory gives larger and at very small distances smaller values than Monte Carlo simulations) cancel out for finite beam diameters, compare Fig. 8.

Acknowledgement

We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076).

References and links

1.

F. Voit, J. Schäfer, and A. Kienle, “Light Scattering by Multiple Spheres: Comparison between Maxwell Theory and Radiative-Transfer-Theory Calculations,” Opt. Lett. 34, 2593–2595 (2009). [CrossRef] [PubMed]

2.

A. Ishimaru, Wave Propagation and Scattering in Random Media, (Academic Press, New York, 1978).

3.

M. S. Patterson, B. Chance, and B. C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

4.

T. J. Farrell, M. S. Patterson, and B. C. Wilson, “A Diffusion Theory Model of Spatially Resolved, Steady-State Diffuse Reflectance for the Noninvasive Determination of Tissue Optical Properties in Vivo,” Med. Phys. 19, 879–888 (1992). [CrossRef] [PubMed]

5.

D. Contini, F. Martelli, and G. Zaccanti, “Photon Migration through a Turbid Slab Described by a Model Based on Diffusion Approximation. I.Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

6.

S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

7.

B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

8.

A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

9.

A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]

10.

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

11.

http://www.uni-ulm.de/ilm/index.php?id=10020200.

12.

A. Kienle and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]

13.

R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary Conditions for the Diffusion Equation in Radiative Transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

14.

E. Meissel and G.B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics, (Bibliobazaar, 2008).

15.

A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Nonin-vasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

16.

H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, (Clarendon, Oxford, 1959).

17.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, (Dover Publications, New York, 1971).

18.

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

OCIS Codes
(170.0170) Medical optics and biotechnology : Medical optics and biotechnology
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(170.5280) Medical optics and biotechnology : Photon migration

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: February 12, 2010
Revised Manuscript: April 7, 2010
Manuscript Accepted: April 14, 2010
Published: April 21, 2010

Virtual Issues
Vol. 5, Iss. 9 Virtual Journal for Biomedical Optics

Citation
André Liemert and Alwin Kienle, "Light diffusion in a turbid cylinder. I. Homogeneous case," Opt. Express 18, 9456-9473 (2010)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-9-9456


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. F. Voit, J. Schäfer, and A. Kienle, “Light Scattering by Multiple Spheres: Comparison between Maxwell Theory and Radiative-Transfer-Theory Calculations,” Opt. Lett. 34, 2593–2595 (2009). [CrossRef] [PubMed]
  2. A. Ishimaru, Wave Propagation and Scattering in Random Media, (Academic Press, New York, 1978).
  3. M. S. Patterson, B. Chance, and B. C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]
  4. T. J. Farrell, M. S. Patterson, and B. C. Wilson, “A Diffusion Theory Model of Spatially Resolved, Steady-State Diffuse Reflectance for the Noninvasive Determination of Tissue Optical Properties in Vivo,” Med. Phys. 19, 879–888 (1992). [CrossRef] [PubMed]
  5. D. Contini, F. Martelli, and G. Zaccanti, “Photon Migration through a Turbid Slab Described by a Model Based on Diffusion Approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]
  6. S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Path length in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]
  7. B. W. Pogue, and M. S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]
  8. A. Kienle, “Light Diffusion through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]
  9. A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]
  10. A. Liemert, and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).
  11. http://www.uni-ulm.de/ilm/index.php?id=10020200.
  12. A. Kienle, and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]
  13. R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary Conditions for the Diffusion Equation in Radiative Transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
  14. E. Meissel, and G. B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics, (Bibliobazaar, 2008).
  15. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]
  16. H. S. Carslaw, and J. C. Jaeger, Conduction of Heat in Solids, (Clarendon, Oxford, 1959).
  17. M. Abramowitz, and I. A. Stegun, Handbook of Mathematical Functions, (Dover Publications, New York, 1971).
  18. A. Liemert, and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited