Home

Search

Collections

Journals

About

Contact us

My IOPscience

An extended analytical approach for diffuse optical imaging

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Med. Biol. 60 5103 (http://iopscience.iop.org/0031-9155/60/13/5103) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 128.255.6.125 This content was downloaded on 18/06/2015 at 05:42

Please note that terms and conditions apply.

Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 5103–5121

Physics in Medicine & Biology doi:10.1088/0031-9155/60/13/5103

An extended analytical approach for diffuse optical imaging H Erkol1,2, F Nouizi1, M B Unlu2 and G Gulsen1 1

  Center for Functional Onco Imaging, Department of Radiological Sciences, University of California, Irvine, CA 92697, USA 2   Department of Physics, Bogazici University, Bebek, 34342, Istanbul, Turkey E-mail: [email protected] and [email protected] Received 29 March 2015, revised 15 April 2015 Accepted for publication 28 April 2015 Published 17 June 2015 Abstract

In this work, we introduce an analytical method to solve the diffusion equation  in a cylindrical geometry. This method is based on an integral approach to derive the Green’s function for specific boundary conditions. Using our approach, we obtain comprehensive analytical solutions with the Robin boundary condition for diffuse optical imaging in both two and three dimensions. The solutions are expressed in terms of the optical properties of tissue and the amplitude and position of the light source. Our method not only works well inside the tissue but provides very accurate results near the tissue boundaries as well. The results obtained by our method are first compared with those obtained by a conventional analytical method then validated using numerical simulations. Our new analytical method allows not only implementation of any boundary condition for a specific problem but also fast simulation of light propagation making it very suitable for iterative image reconstruction algorithms. Keywords: photon propagation, diffuse optical imaging, point source (Some figures may appear in colour only in the online journal) 1. Introduction Diffuse optical imaging (DOI) uses near-infrared (NIR) light to determine optical properties of tissue and has a wide range of diagnostic applications from breast cancer to brain imaging (Tromberg et al 2005, Tromberg 2005, Gibson et al 2005, Cerussi et al 2006, Cerussi et al 2007, Tromberg et al 2008). In DOI, three dimensional images are obtained using the boundary measurements of NIR light (Arridge 1999, Dehghani et al 1999, Dehghani et al 2000, Hielscher et al 2002, Pogue et al 2003, Pogue 2006, Unlu and Gulsen 2008, Unlu et al 0031-9155/15/135103+19$33.00  © 2015 Institute of Physics and Engineering in Medicine  Printed in the UK

5103

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

2008a, 2008b, Valim et al 2010, Valim et al 2013). To accomplish this, sources and detectors are placed around tissue and measurements are recorded in reflection or transmission mode. While light propagates, its interaction with biological tissue occurs due to absorption and elastic scattering. Modelling of light propagation in tissue can be achieved accurately by radiative transport equation (RTE) (Darne et al 2014). However, photon transport in turbid medium is generally approximated by the diffusion equation because the RTE is very hard to solve and scattering in tissue is remarkably larger compared to the absorption (Patterson 1989, Schmitt et al 1990, Patterson and Madsen 1991, Arridge et al 1992, Farrel et al 1992, Arridge et al 1993, Schweiger et al 1993, Boas et al 1994, Wang and Wu 2007). Photon density throughout the region can be obtained by utilizing analytical, Monte Carlo and numerical methods (such as the Finite Element Method and the Finite Difference Method) (Flock et al 1989, Arridge et al 1993, Das et al 1996). Analytical and Monte Carlo methods yield considerably accurate results compared to numerical methods. Furthermore, analytic methods provide the fastest computational time. Solving diffusion equation analytically is generally a challenging task. Nevertheless, analytical solutions have been presented for some definite geometries such as spherical or cylindrical as well as sufficiently large or infinite slabs (Patterson 1989, Schmitt et al 1990, Patterson and Madsen 1991, Arridge et al 1992, Farrel et al 1992, Arridge et al 1993, Schweiger et al 1993, Boas et al 1994, Pogue and Patterson 1994, Contini et al 1997, Kienle and Patterson 1997, Walker et al 1998, Martelli et al 2002, Cong et al 2004, Kienle 2005, Martelli et al 2007, Liemert and Kienle 2010, Zhang et al 2010). In addition, photon propagation, which is obtained by solving the diffusion equation, has also been widely studied for homogeneous sub-layers since some imaging techniques model tissue as a layered structure (Takatani and Graham 1979, Schmitt et al 1990, Dayan et al 1992, Okada et al 1997, Schweiger and Arridge 1997, Svaasand et al 1999, Kienle and Glanzmann 1999, Pogue et al 1999, Dehghani and Delpy 2000, Martelli et al 2002, Dehghani et al 2005, Yalavarthy et al 2006, Sikora et al 2006, Martelli et al 2007, Liemert and Kienle 2010a, 2010b). However, since analytical solutions for irregular geometries have not been obtained, numerical approximation methods such as the Finite Element Method (FEM) and the Finite Difference Method (FDM) are used (Arridge et al 1993, Schweiger et al 1993). Up to date, the Green’s function approach has been mainly used to solve the diffusion equation  in regular geometries (Patterson 1989, Arridge et al 1992). For example, Arridge et al (1992) obtained solutions for an infinite cylinder in which an infinitely long line source is placed. They also derived solutions of the diffusion equation for a point source utilizing the Green’s function technique in various regular geometries. Furthermore, Sikora et al (2006) used a series expansion approach to solve the diffusion equation for concentric spheres. In another work, Liemert and Kienle (2010b) obtained continuous wave (cw), frequency and time domain solutions for the diffusion equation via the Green’s function method with the extrapolated boundary conditions for a multiple layered finite cylinder. Moreover, Zhang et al (2010) presented a cw solution for a point source utilizing the extrapolated boundary conditions in cylindrical coordinates. In addition to those, Liemert and Kienle obtained detailed solutions for the diffusion equation  in a homogeneous turbid medium with a point source using various integral transformations (Liemert and Kienle 2010). In this paper, we obtain both two and three dimensional solutions for the diffusion equation analytically considering steady state (cw) case in a cylindrical medium. Here the light source is modeled by the Dirac δ function with a given strength. We present an integral approach to derive the Green’s function for the Robin boundary condition. Our method is indeed very flexible allowing implementation of any boundary condition (i.e. not limited with the Robin boundary condition). Our approach can also be applied to other regular geometries 5104

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

such as spherical. The main motivation of our study is to obtain solutions for the diffusion equation at the boundary making our method very suitable for DOI of homogeneous or nearly homogeneous media. To be able to validate our method, we first compare it with the analytical method presented by Arridge et al (1992). Since Arridge et al (1992) utilize the known Green’s function of the infinite medium with the zero boundary condition, it corresponds to a special case for our solution. In addition, we compare our results with the solutions obtained by the finite element method in both two and three dimensions. 2.  Theoretical method The photon propagation in tissue is described by the time dependent diffusion equation in time domain (Wang and Wu 2007, Nouizi et al 2011) ∂Φ(r, t ) + μa Φ(r, t ) − ∇ ⋅ [D∇Φ(r, t )] = S (r, t ) (1) c∂t

where Φ, c, D, μa, and S represent the photon density, the speed of light, the diffusion coefficient, the absorption coefficient and the source term, respectively. Here, the diffusion coefficient is given as D(r) = 1/(3[μa (r) + μs′(r)]) where μs′ is the reduced scattering coefficient. The light source can be considered a point source since it is very small. Hence, the source term Λ can be approximated by the Dirac delta function. As a result of this, we define S (r ) = c δ (r, r′) where r′ is the position of the light source and Λ is a dimensionless parameter determining the Λ

source strength. For calculational simplicity, we define γ = c . For a steady state case (continuous wave), the diffusion equation takes the following form μ γ 2 −∇ Φ(r ) + a Φ(r ) = δ (r, r′) (2) D D where D is taken as constant. 2.1.  Two dimensional case

In two dimensional cylindrical coordinates, except for the source position r  ≠  ri, the solution of the diffusion equation is ∞

Φ(r , θ ) = ∑ [amJm(k r ) + bmYm(k r )][cm cos(mθ ) + dm sin(mθ )] (3) m =−∞

where k = i

μa D

and i is the complex number, Jm and Ym are the first and second kind Bessel

functions, respectively. Here am, bm, cm, and dm are the differentiation constants. If the source is placed along the x axis (θ  =  0), the photon density is symmetric with respect to this axis. In other words, the constant dm in equation (3) becomes zero so that the angular dependency of the solution comes from cos (mθ) term leading to ∞

Φ(r , θ ) = ∑ [amJm(k r ) + bmYm(k r )] cos(mθ ). (4) m =−∞

Now we solve the diffusion equation for a point source using the properties of the Dirac δ function (figure 1). Firstly, we divide the region into two sub-regions according to the position of the source. If r  ⩽  ri, the solution of the radial part consists of only the first kind Bessel 5105

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

Figure 1.  The schematic showing the geometry of a homogeneous medium with a δ

function source in two dimensions.

function Jm(kr) since the second kind Bessel function goes to infinity at r  =  0. Hence, the solution reduces to ∞

Φ (r , θ ) = ∑ [bmJm(kr ) + cmYm(kr )] cos(mθ ) (6) m =−∞

r  ⩾  ri. To relate constants bm and cm, we apply the following Robin boundary condition at the boundary surface, r  =  R: ∂Φ(r , θ ) ∣r = R = 0 R, θ ) + 2ξD Φ( (7) ∂r

where ξ is a constant corresponding to refractive index mismatched between tissue and its surrounding medium (Arridge 1999, Gao et al 2010, Liu et al 2010). Now we utilize the properties of the Dirac delta function. The first property is that the solutions of the two regions are equal to each other at the source position. In other words, the expression for the equality of the photon density at r  =  ri is Φ (8) (r , θ )∣r = ri .

To represent the sudden change in the derivative of photon density at the source position, let us write solutions into the diffusion equation ⎤ ⎧1 d ⎡ d ⎫ 2 μ ⎢r Gm(kr )⎥ cos(mθ ) + 1 d [cos(mθ )]Gm(kr ) − a Gm(kr ) cos(mθ )⎬ = − γ δ (r, ri) ⎨ 2 2 ⎥⎦ r dr ⎢⎣ dr r dθ D D ⎭ m =−∞⎩ ∞











(9) where Gm(kr) is the following set of Bessel function solutions, Gm(kr)={Jm(kr), Ym(kr)}. Multiplying equation  (9) with r cos (m′θ) and then integrating between (r  =  ri  −  ϵ, θ  =  0) and (r  =  ri  +ϵ, θ  =  2π) yields 5106

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103



ri + ϵ

∑ ∫ r −ϵ i

m =−∞



=−

1 γ D

⎡ d ⎛ dG (kr ) ⎞⎤ ⎢ ⎜⎜r m ⎟⎥dr ⎢⎣ dr ⎝ dr ⎟⎠⎥⎦

ri + ϵ

∫r −ϵ ∫0 i



∫0



cos(mθ ) cos(m′θ )dθ

δ (r, ri) cos(m′θ )r dr dθ.

(10)

Here notice that except from the first term on the left hand side of equation (9), all the terms go to zero due to the fact that Φ itself is continuous. Next we consider ϵ  →  0+ and use the equality of the photon density at the source position given by equation (8). Writing δ (r − ri )δ (θ − θi ) δ(11) (r, ri) = r

and the orthogonality of cosine functions 2π

cos(mθ ) cos(nθ )dθ = πδmn ∫(12) 0

(if m, n  =  0, equation (12) is 2π) into equation (10) as ϵ  →  0+ and using the equality relation of Φ gives the expression for the abrupt change of

∂Φ : ∂r

γ dGm(kr ) dG (kr ) cos(mθi ). ∣r = ri + ϵ − m ∣r = ri − ϵ = − (13) dr dr πDri

Writing the Robin boundary condition for Φ> and substituting Φ> (r, θ) and Φ< (r, θ) into equations (8) and (13) gives ⎡ dJ (kr ) dY (kr ) ⎤ bmJm(k R ) + cmYm(k R ) = −2ξD⎢bm m + cm m (14) ⎥ ∣r = R , ⎣ dr dr ⎦ bmJm(k ri ) + cmYm(k ri ) = amJm(k ri ) (15)

and ⎡ dJm(kr ) γ dY (k r ) dJ (k r ) ⎤ + cm m − am m cos(mθi ). (16) ⎢bm ⎥ ∣r = ri = − ⎣ ⎦ πDri dr dr dr

Solving equations (14)–(16) for the differentiation constants am, bm and cm yields the photon density ⎧ amJm(kr ) cos(mθ ) if r ⩽ ri Φ( r, θ ) = ∑ ⎨ (17) ⎩ (bmJm(kr ) + cmYm(kr )) cos(mθ ) if  r ⩾ ri ∞

m =−∞

where k = i

μa , D

γ cos(mθi ) 2D(2DkξRJm − 1(kR ) + (R − 2Dmξ )Jm(kR )) (18) × (Jm(kri )(2DkξRYm − 1(kR ) + (R − 2Dmξ )Ym(kR)) am =

+Ym(kri )(−2DkξRJm − 1(kR ) − (R − 2Dmξ )Jm(kR ))), Dkξ(Ym − 1(kR ) − Ym + 1(kR )) + Ym(kR ) bm = γR 2 D (2DkξRJm − 1(kR ) + (R − 2Dmξ )Jm(kR)) (19) × Jm(kri ) cos(mθi ) 5107

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

and γ cm = − Jm(kri ) cos(mθi ). (20) 2D

These two dimensional results can be used for obtaining the photon density around the center when length of the cylindrical volume under investigation is very large. A special case corresponding to the zero boundary condition can be obtained by setting ξ  =  0 in the Robin boundary condition. Since this solution has been already presented in detail by Arridge et al (1992), this approach will allow us to validate our new method. In this case, am, bm and cm reduce to ⎤ γ ⎡ Jm(kri )Ym(kR ) am = − Ym(kri )⎥ cos(mθi ), (21) ⎢ ⎦ Jm(kR ) 2D ⎣ γ Jm(kri )Ym(kR ) cos(mθi ), bm = (22) 2D Jm(kR )

and γ cm = − Jm(kri ) cos(mθi ), (23) 2D

respectively. Now let us use the following relations between the Bessel functions (Jm and Ym) and the modified Bessel functions (Im and Km): Im(x ) = i−mJm(ix ) (24)

and π m+1 K i [Jm(ix ) + iYm(ix )]. m (x ) = (25) 2 μa D

Defining k = i

≡ iα and then substituting

Jm(iαr ) = i mIm(αr ) (26)

and 2 m+1 Y(27) Im(αr ) m(iαr ) = − m Km(αr ) + i πi

into equations (21)–(23) yields the photon density for r  >  ri given by equation (17)

∞ ⎫ 1 I (αr ) ⎧ Φ ∑ m i ⎨Km(αr )Im(αR) − Km(αR)Im(αr )⎬ cos(mθ ) (28) >(r , θ ) = πD m =−∞ Im(αR ) ⎩ ⎭ ⎪







by placing the source along the x axis (θi  =  0). A comprehensive comparison with the solution given by Arridge et al (1992), can be obtained by transforming the differentiation constants am, bm and cm to am′ , bm′ and cm′ and by using these new constants in equation (17), respectively. The new differentiation constants are am′ = am /( 8π ), bm′ = bm /( 8π ) and cm′ = cm /( 8π ), respectively. In other words, we take the scaled photon density Φscaled = Φ/( 8π ). Therefore, we obtain the exact solution as the Arridge et al's solution (Arridge et al 1992): Φ>(r , θ ) =



Im(αri )  {Km(αr )Im(αR ) 3 I (2π )2 D m =−∞ m(αR ) 1





5108

− Km(αR )Im(αr )} cos(mθ )

(29)

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

for r  >  ri. By interchanging r  →  ri, we can also obtain Φ (r,θ,z) and Φ< (r,θ,z) into equations (40) and (42) gives

⎛ ⎞ ⎛ 2 ⎞ ⎛π ⎛ ⎛π ⎛ 1 ⎞⎞ 2 1 ⎞⎞ bm, lJm⎜⎜ k 2 −⎜ ⎜l + ⎟⎟ R⎟⎟ + cm, lYm⎜⎜ k 2 − ⎜ ⎜l + ⎟⎟ R⎟⎟  ⎝L⎝ 2 ⎠⎠ ⎠ 2 ⎠⎠ ⎠ ⎝L⎝ ⎝ ⎝ ⎡ ⎛ ⎞ ⎛ ⎞⎤ ⎛π ⎛ ⎛π ⎛ 1 ⎞⎞2 1 ⎞⎞2 + 2ξD⎢bm, lJm′ ⎜⎜ k 2 −⎜ ⎜l + ⎟⎟ R⎟⎟ + cm, lYm′ ⎜⎜ k 2 −⎜ ⎜l + ⎟⎟ R⎟⎟⎥= 0, ⎢ 2 ⎠⎠ ⎠ 2 ⎠⎠ ⎠⎥⎦ ⎝L⎝ ⎝L⎝ ⎝ ⎝ ⎣  ⎛ ⎞ ⎛ ⎛π ⎛ ⎛π ⎛ 1 ⎞⎞ 2 bm, lJm⎜⎜ k 2 − ⎜ ⎜l + ⎟⎟ ri⎟⎟ +cm, lYm⎜⎜ k 2 − ⎜ ⎜l + ⎝ ⎝ ⎠ ⎠ ⎝L⎝ 2 L ⎝ ⎠ ⎝  ⎛ ⎛π ⎛ = am, lJm⎜⎜ k 2 − ⎜ ⎜l + ⎝L⎝ ⎝

(43)

⎞ 1 ⎞⎟⎟⎞2 ⎟ ri⎟ 2 ⎠⎠ ⎠

⎞ 1 ⎞⎟⎟⎞2 ⎟ ri⎟ 2 ⎠⎠ ⎠

(44)

and ⎛ ⎞ ⎛ ⎞ ⎛π ⎛ ⎛π ⎛ 1 ⎞⎞ 2 1 ⎞⎞ 2 bm, lJm′ ⎜⎜ k 2 − ⎜ ⎜l + ⎟⎟ ri⎟⎟ + cm, lYm′ ⎜⎜ k 2 − ⎜ ⎜l + ⎟⎟ ri⎟⎟ ⎝L⎝ ⎝L⎝ 2 ⎠⎠ 2 ⎠⎠ ⎝ ⎠ ⎝ ⎠  ⎛ ⎞ ⎛ ⎛ ⎞ ⎞ ⎛π ⎛ γ π 1 1 ⎞⎞2 −am, lJm′ ⎜⎜ k 2 − ⎜ ⎜l + ⎟⎟ ri⎟⎟ = − cos(mθi ) cos⎜⎜ ⎜⎜l + ⎟⎟zi⎟⎟ (45) ⎝L⎝ 2⎠ ⎠ πDLri 2 ⎠⎠ ⎝ ⎠ ⎝L⎝ 5111

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

where ′ stands for the differentiation with respect to r. Solving equations (43)–(45) for am, l, bm, l and cm, l yields the photon density

⎧ ⎛π ⎛ 1⎞ ⎞ if r ⩽ ri am, lJm(κr ) cos(mθ ) cos⎜ ⎜l + ⎟z⎟ ⎪ ⎪ ⎝L⎝ 2⎠ ⎠ Φ(r , θ , z ) = ∑ ⎨ ⎛π ⎛ 1⎞ ⎞ m, l =−∞⎪ (b J (κr ) + c Y (κr )) cos(mθ ) cos⎜ ⎜l + ⎟z ⎟ if  r ⩾ r m, l m i ⎪ m, l m ⎝ ⎝ 2⎠ ⎠ L ⎩  ∞

where κ =

(46)

⎛π ⎛ 1 ⎞⎞ 2 k 2 − ⎜ ⎜l + ⎟⎟ for the calculational simplicity and ⎝L⎝ 2 ⎠⎠ γ cos(mθi ) cos

( (l + )z ) π L

1 2

i

a m, l = 2DL (2DκξRJm − 1(κR ) + (R − 2Dmξ )Jm(κR )) (47) × (Jm(κri )(2DκξRYm − 1(κR ) + (R − 2Dmξ )Ym(κR )) +Ym(κri )(−2DκξRJm − 1(κR ) − (R − 2Dmξ )Jm(κR ))), bm, l = γR



Dκξ(Ym − 1(κR ) − Ym + 1(κR )) + Ym(κR ) 2DL (2DκξRJm − 1(κR ) + (R − 2Dmξ )Jm(κR ))

⎛π ⎛ 1⎞ ⎞ × Jm(κri ) cos(mθi ) cos⎜ ⎜l + ⎟zi⎟ ⎝L⎝ 2⎠ ⎠

(48)

and ⎛π ⎛ γ 1⎞ ⎞ cm, l = − Jm(κri ) cos(mθi ) cos⎜ ⎜l + ⎟zi⎟. (49) ⎝ ⎝ 2⎠ ⎠ 2DL L

Similar to the two dimensional case, we also obtain a solution corresponding to the zero boundary condition for the three dimensional case to be able to compare our method with Arridge et al's method and validate it. Setting ξ  =  0 in the Robin boundary condition, the differentiation constants reduce to ⎤ ⎛π ⎛ γ ⎡ Jm(κri )Ym(κR) 1⎞ ⎞ − Ym(κri )⎥ cos(mθi ) cos⎜ ⎜l + ⎟zi⎟, a m, l = (50) ⎢ ⎦ ⎝L⎝ 2⎠ ⎠ 2D ⎣ Jm(κR) ⎛π ⎛ γ Jm(κri )Ym(κR ) 1⎞ ⎞ bm, l = cos(mθi ) cos⎜ ⎜l + ⎟zi⎟, (51) ⎝L⎝ Jm(κR ) 2⎠ ⎠ 2D

and ⎛π ⎛ γ 1⎞ ⎞ cm, l = − Jm(κri ) cos(mθi ) cos⎜ ⎜l + ⎟zi⎟, (52) ⎝L⎝ 2⎠ ⎠ 2D

respectively. Now we define κ =

⎛π ⎛ 1 ⎞⎞ 2 k 2 − ⎜ ⎜l + ⎟⎟ ≡ iαl for the calculational simplicity. ⎝L⎝ 2 ⎠⎠

We substitute Jm(iαlr ) = i mIm(αlr ) (53) 5112

H Erkol et al

Phys. Med. Biol. 60 (2015) 5103

and 2 m+1 Y(54) Im(αlr ) m(iαlr ) = − m Km(αlr ) + i πi

into equations (50)–(52). Next we transform the differentiation constants am,l, bm,l and cm,l to a m, l

b

c

, m,l and m,l , or equivalently, we take the scaled solution Φscaled = Φ/( 8π ). Therefore, 8π 8π the solutions are 8π

Φ>(r , θ , z ) =

1 π (2π )1/2DL





∑ ∑

m =−∞ l =−∞

Im(αlri ) {Km(αlr )Im(αlR ) − Km(αlR )Im(αlr )} Im(αlR ) ⎛π ⎛ 1⎞ ⎞ × cos(mθ ) cos⎜ ⎜l + ⎟z⎟ ⎝L⎝ 2⎠ ⎠



(55)

and Φ  ri and r  

An extended analytical approach for diffuse optical imaging.

In this work, we introduce an analytical method to solve the diffusion equation in a cylindrical geometry. This method is based on an integral approac...
2MB Sizes 0 Downloads 10 Views