1. Introduction

Light can penetrate to several centimeters in biological tissue owing to the absence of chromophores that could absorb radiation in the near infrared region (NIR) [1]. Recently, information extracted from such photon transport studies have been used increasingly in medical diagnosis [1]. Notably, photon transport in biological media has been used for interesting applications such as imaging of premature infant brain [2], characterization of dentin and enamel [3], investigation of developing embryonic morphology [4] and noninvasive detection of glucose using spectroscopic methods [5]. However, due to strong scattering relative to absorption, imaging and non-invasive characterization in NIR becomes a very difficult task [1, 6]. Therefore, having a clear understanding and a reliable model of photon transport in turbid biological media plays a significant role in improving techniques used for medical diagnosis [1, 7, 8].

Researchers have used variety of numerical strategies for approximately solving standard photon transport equation in biological media such as Monte Carlo simulation techniques [8, 9], discrete ordinate methods [10], diffusion theory [8, 11], finite element methods [12] and various flux models [8] (Kubelka-Munk theory [6], four flux model [13], seven flux model [13]). In Monte Carlo techniques, photons are injected from a source either individually or as a bunch (or a packet) and traced through the scattering/absorptive medium following photon ray paths [9]. Each photon path incurs a random set of scattering and absorption events. A large number of photons is used to predict an average path. Even though Monte Carlo techniques can provide high accuracy result, due to high computation cost, they are used mainly for validating alternative, faster techniques. In the discrete ordinate method, the transport equation is discretized in different directions (or ordinates) and solved using accurate quadrature techniques [10]. However, similar to Monte Carlo techniques, prohibitive computational cost has reduced its usage in practice. The diffusion equation provides an approximation to the transport equation when the number of scattering events is high. Given that biological media is strongly scattering relative to absorption, diffusion equation has the ability to give very accurate results in practice [8, 11]. Apart from requiring a low computational effort, the diffusion equation can also be solved analytically [8]. Recently, finite element methods have also attracted much attention because of their ability to account for complex geometries and inhomogeneous media [12].

It is interesting to observe that up to now, these techniques have been mainly used for solving homogeneous media with constant refractive indices. However, most of the these techniques have the ability to account for refractive index discontinuities using Snell’s Law and Fresnel equations for plane boundaries [8, 9]. At a refractive index discontinuity, the incident photon beam or packet is resolved into reflected and refracted wave using Fresnel equations [14] and traced forward as two separate photon beams or packets with modified intensity or photon number. However, this technique fails to satisfactorily represent the photon transport in turbid media where refractive index continuously changes along the spatial dimensions [14]. This is mainly because when the refractive index is continuously changing along a certain direction, light rays propagating in that direction do not reflect any energy in opposite direction but tend to follow a curved path depending on the gradient of the refractive index along the path [14, 15]. Recently, Ferwerda [16] suggested a modification to standard radiative transport equation to account for media with isotropic, spatially varying continuous refractive index profiles. Khan and Jian [17] identified some redundant terms in Ferwerda’s derivation [16] and hence gave a modified, but equivalent form of equation. They also derived the analogous diffusion approximation for Ferwerda’s [16] formulation.

In this paper, we show that Ferwerda’s [16] (and hence Khan and Jian’s [17]) formulations do not satisfy energy conservation even in lossless, nonscattering media and fail to account for wavefront curvature. Using the principle of energy conservation [10, 14] and laws of geometrical optics [14, 15], we derive the photon transport equation for turbid biological media with spatially varying isotropic refractive index. We show that when the refractive index is constant, our result reduces to the standard radiative transfer equation [18, 16] and when the medium is lossless and free of scattering to the well known geometrical optics equations in refractive media [19, 15, 14].

This paper is organized as follows. In the Section 2, we track a photon packet propagating through isotropic, continuously spatially varying refractive index profile in scattering/ absorptive turbid biological media. In Section 2.1, we show that photon propagation can be viewed as a coordinate transformation or mapping through the turbid medium and hence the Jacobian [20] of the transformation plays a vital part in the characterization of the photon flow. In Section 2.2, using the definition of the Jacobian, we derive its differential properties when the coordinate transformation is governed by the laws of geometrical optics. Using the properties of the Jacobian for ray congruences [14], we derive two different but equivalent formulations of photon transport equation in turbid biological media. The first format given in Section 2.2 is based on the arc length parameter of the photon trajectory (or light path) and also demonstrates the conservative features of the result. The second formulation given in Section 2.3 is the standard form used in literature [18, 16, 17, 10] where Cartesian position coordinates and spherical propagation direction coordinates are used. In Section 3, we look at special cases where results are known from independent analysis. As expected, we show in Section 3.1 that the current formulation reduces to the well-known standard photon transport equation [18, 10] when the refractive index is homogeneous. In Section 3.2, we investigate the behavior of the photon transport equation for a lossless and non-scattering medium and show the results completely agrees with what is expected from geometrical optics analysis [19, 15, 14]. In Section 4, we conclude the paper after summarizing results and discussing further research directions.

2. The photon transport equation with spatially varying isotropic refractive index

To derive the photon transport equation (TPE) with spatially varying refractive index, we track the evolution of an infinitesimally small convex elementary scattering volume along a infinitesimally small ray tube surrounding a central ray in phase space. Figure 1 shows the evolution of initial photon packet with phase space volume, $\delta V_{S_0}$ at arc length s=s ₀. The arc length s is measure from an arbitrary but fixed reference point (relative to the laboratory frame) of the central light ray.

To proceed further we introduce the rectangular cartesian coordinate system shown in Fig. 2. Positions along the right ray are given by the vector, r=x x̂+y ŷ+z ẑ where x̂, ŷ and ẑ are unit vectors along x, y and z directions, respectively. We also assume that tangent vector at position r is given by $\widehat{\Omega }$ ≡Ω_x x̂+Ω_y ŷ+Ωz ẑ with polar angle, θ with the ẑ direction and azimuthal angle, φ with the x̂ direction as shown in Fig. 2.

 

Fig. 1. Transport of an infinitesimally small photon packet with phase space volume, $V_{S_0}$, along an infinitesimally small ray tube surrounding a central light ray.

Download Full Size | PDF

 

Fig. 2. Details of the stationary and moving coordinate systems

Download Full Size | PDF

The unit tangent vector, $\widehat{\Omega }$ , can be then written as:

where φ∊[0,2π), θ∊[0, pi] and µ=cos(θ) giving the range µ∊[-1,1]. Suppose ℒ(r,$\widehat{\Omega }$ ,t)≡ℒ(x,y,z,µ,φ,t) is the radiance of the photon flow [10, 18], which describes the amount of energy in the photon beam at position r, flowing per second through a unit area perpendicular to the tangential unit vector $\widehat{\Omega }$ . Therefore, the integral

is proportional to the change in photon number inside the phase space volume V per unit time. Applying this result to the case given in Fig. 1, we track an infinitesimally small photon packet at s=s ₀ with phase space volume, $\delta V_{S_0}$ along the trajectory of a light ray. We basically follow the approach taken by Ferwerda [16] in his derivation of photon transport equation with spatially varying, isotropic refractive index media. However, our approach is systematic and rigorous, leading to a different result.

We observe that the photon number in this initial volume can grow due to incoming photons from distributed sources along the light path and photons scattering to the volume from outside, whereas photon number in this initial volume can reduce due to attenuation and scattering of photons to outside. Application of energy conservation principle [10, 18, 16] (i.e. photon number conservation in this case due to tracking of a single wavelength) to this moving volume leads, for a small Δs length along the trajectory, the following result:

where µ_a≡µ_a(x,y,z) is the absorption coefficient per unit length and µ_s≡µ_s(x,y,z) is the scattering coefficient per unit length. ε≡ε(x,y,z,µ,φ,t) represents the source distribution per unit volume [16, 18]. The probability of scattering a photon moving in direction (µ̀,φ̀) to (µ,φ) is give by the normalized phase function, f (µ,φ,µ̀,φ̀) [18, 16]. Noting that this infinitesimal photon packet move along ray congruences [14], the photon packet volume sizes, δV _s+Δs and δV_s, at s and s+Δs, respectively, can be seen as a mapping of initial volume $\delta V_{S_0}$ along the ray congruences [14, 15]. In multidimensional space, such mapping happens through the Jacobian, J≡J(x,y,z,µ,φ) [20]. Using this relation, we could get the following expression for the left hand side of Eq.(3):

where (x ₀,y ₀,z ₀,µ ₀,φ ₀) is the coordinate variables (x,y,z,µ,φ) at s=s ₀. However, for a medium with spatially varying refractive index, n≡n(x,y,z), we have the relation, ds/dt=c/n, where c is the speed of light in vacuum [14]. Using, this, we could get the following relation for the integrand of right hand side of Eq.(4):

Now to evaluate Eq.(4), we need to evaluate differential of Jacobian with respect to arc length at a certain time instant (i.e. $\frac{\partial J}{\partial s}$).

2.1. Evaluation of term $\frac{\partial J}{\partial s}$

Suppose |x| represents the determinant of matrix x. Then we could write the Jacobian J as [20]:

where (x ₀,y ₀,z ₀,µ ₀,φ ₀) is the coordinate variables (x,y,z,µ,φ) at s=s ₀. To evaluate this we introduce the permutation function, σ, of the integer set {1,2, …,n} [21]. Transposition (i, j) of σ is an operation which interchanges i and j but leaves other numbers fixed. It can be shown that every permutation mapping of σ is a combination of transpositions [21]. It can be shown easily that if one representation of a permutation by transposition involves an even number of transpositions, then all representations involving transpositions result in even number of transpositions [21]. An analogous statement applies to the combination of odd numbers of transpositions. Using this property, we denote the sign of a permutation σ, denoted sgn(σ) is +1 if σ is a combination of even number of transpositions and is -1 if σ is a combination of an odd number of transpositions. Using these definitions, we could write J in Eq.(6) as a sum of permutations of σ={1,2,3,4,5} [21]:

If ϑ₀ is a one-to-one, onto map from index set {1,2,3,4,5} to variable set (x ₀,y ₀, z ₀,µ ₀,φ ₀) and ϑ is a similar map to variable set (x,y, z,µ,φ) then we could write Eq.(7) in the following expanded differential form:

Differentiating Eq.(8) respect to s when t is held constant, we get:

Careful inspection of each of the terms under summation index v in Eq.(9) shows that only a single term out of five expressions contributes to the final sum due to repeats of identical rows. Thus, we have the following result after some tedious but straightforward calculations:

To write this equation in compact form we could use the definitions of divergence operator in a rectangular coordinate system. For the cartesian system (x̂, ŷ, ẑ) the divergence operator of the vector function V=V_x x̂+V_y ŷ+V_z ẑ has the form [22]:

Using this definition, we could write Eq.(10) as:

However, $\widehat{\Omega }$ denotes a unit vector in the light propagation direction of the ray’s trajectory. It is a well known result in geometric optics that the divergence of $\widehat{\Omega }$ is given as the sum of principal radii of curvatures R ₁(s) and R ₂(s) of the geometrical wave-fronts [19, 14, 15]:

It is also well known from geometric optics that the light ray trajectory is given by [14, 15]:

where n≡n(x,y, z) is the spatially varying refractive index of the propagation medium. After carrying out the differentiation with respect to s, Eq.(14) can be rearranged to get

Because $\widehat{\Omega }$ is a unit vector (i.e. $\widehat{\Omega }$ ·$\widehat{\Omega }$ =1), it is straightforward to show that ∂ $\widehat{\Omega }$ /∂ s is perpendicular to the vector $\widehat{\Omega }$ [20, 22]. Therefore we introduce following unit vectors, µ̂ and φ̂ in the plane perpendicular to the unit vector $\widehat{\Omega }$ :

Differentiating Eq.(1) with respect to s, it possible to derive the following result:

Substituting Eq.(18) to Eq.(15) and projecting results along µ̂ and φ̂, we derive the following expressions for the rate of change of µ and φ along the light trajectory:

The unit vectors µ̂ and φ̂ in Eq.(16) and Eq.(17), respectively have the following rate changes along µ and φ directions.

Differentiating Eq.(19) and Eq.(20) with respect to µ and φ, respectively, we get the following results:

Substitution of these results and Eq.(13)to Eq.(12), we finally get the following compact equation for the differential of Jacobian J along the ray trajectory:

Noting that $\frac{dn}{ds}={\nabla }_{\mathbf{r}}n·\hat{\Omega }$, Eq.(26) can be written in the following equivalent form:

2.2. Photon transport equation in light path trajectory arc length format

We use the results in previous section to find the photon transport equation along a light path using arc length variable, s. Substitution of Eq.(26) and Eq.(5) in Eq.(4), we get the following relation, integrated over the initial phase space volume, $\delta V_{S_0}$:

Recognizing that jacobian, J maps the infinitesimal volume $\delta V_{S_0}$ to δV_s, we could rewrite Eq.(27) as:

Using this result in Eq.(3) and taking the limit as Δs→0, we get the integral equation:

This result is valid for any arbitrary infinitesimal phase space volume, δV_s. Therefore, the integrand of the above equation must be identically equal to zero [22].

where we used the following equivalent relations in function arguments: (x,y,z)→r; (µ,φ)→$\widehat{\Omega }$ and (µ̀,φ̀)→$\stackrel{`}{\Omega }$ . Eq.(30) is the main result of this paper. It is important to note that both r and $\widehat{\Omega }$ are functions of s with functional relation given by Eq.(14). Therefore, once the spatially dependent, isotropic refractive index relation, n(r) is given, photon trajectory can be constructed by solving Eq.(14). This result can then be used to calculate the eikonal, 𝒮 using the relation [14, 15]:

Once the eikonal, 𝒮 is known, principal radii of curvatures, R ₁(s) and R ₂(s) can be calculated using standard techniques in differential geometry [23].

2.3. Photon transport equation in (r,$\widehat{\Omega }$ ) format

It is also useful to express the photon transport equation given in Eq.(30) using (r,$\widehat{\Omega }$ ) variables because this is the standard form that researchers have reported results in the past [16, 17, 18]. Moreover, this format is much suited for numerical calculations [10]. Noting that ℒ(r,$\widehat{\Omega }$ )/n ²(r) is a function of both r and $\widehat{\Omega }$ , we could write:

where we used the definition for gradient operator in the moving orthogonal coordinate system ($\widehat{\Omega }$ ,µ̂,φ̂) [22]:

Using the Eq.(15) for ∂ $\widehat{\Omega }$ /∂s and noting the identities µ̂·φ̂=0, Ω·$\widehat{\nabla }$ _{$\widehat{\Omega }$}=0 and ∂/∂s=$\widehat{\Omega }$ ·∇_r, we could write Eq.(32) in following equivalent form:

Substitution of this result to Eq.(30) gives the photon transport equation in (r,$\widehat{\Omega }$ )format as:

It is instructive compare this result with Ferwerda [16](and so Khang and Jiang [17]). It is easy show that their [16, 17] calculation of ∇_r·Ω is incorrect because it depends only on differentials of the refractive index and hence does not take wavefront curvature into account. It is a well known fact that even in constant refractive index media, depending on the source configuration, wave fronts can have finite radius of curvature [14, 15]. However, Ferwerda’s ∇_r·Ω term reduces to zero in such cases. Also, his result does not have the term (-2/n(r))$\widehat{\Omega }$ ·∇_r n(r). This term is essential for making the Eq.(35) conserve energy.

3. Special cases

To confirm the validity of proposed photon transport equation given in Eq.(30) or Eq.(35), we investigate its behavior for some special cases where expected outcome is known from independent analysis.

3.1. Photon transport equation in constant refractive index media

It is instructive to see the behavior of Eq.(35) in a constant refractive index media (i.e. n(r)=n ₀) where wavefront curvatures are assumed to be plane. Noting that for such medium, ∇_r n(r)=0 and both principal radii of curvatures, R ₁(s) and R ₂(s) become infinite, we get:

This is the well know photon transport equation for homogeneous turbid media (with constant refractive index) [10, 18, 16].

3.2. Lossless non-scattering medium - steady state case

In a lossless, non-scattering medium, the photon transport equation in Eq.(30) simplifies to:

If the radiance of the photon beam has a value of ℒ=ℒ ₀ at s=s ₀ corresponding to r=r ₀, then Eq.(37) gives:

This result agrees with the geometric optics result given by Born and Wolf [14] and also Kline [19, 15] for wavefronts with finite radius of curvature. Especially when the wavefronts are plane, both principal radii of curvatures, R ₁(s) and R ₂(s) become infinite [14, 15]. Therefore, the integral of Eq.(38) becomes zero, giving the following result:

This result agrees with Ishimaru’s result [18] on plane wave photon transport in turbid media. It is interesting to note that none of the results given in this section can be derived from the photon transport equation reported by Ferwerda [16] or it subsequent modified form by Khan and Jiang [17].

4. Conclusions

In this paper, using the principle of energy conservation and laws of geometrical optics, we derived the photon transport equation for turbid biological media with spatially varying isotropic refractive index. We showed that photon propagation through the medium can be viewed as a coordinate transformation or mapping through the turbid medium and hence the Jacobian [20] of the transformation plays vital part in the characterization of the photon flow. Using the properties of the Jacobian for ray congruences [14], we presented two equivalent representation of the transport equation. The first format is based on the arc length parameter of the photon trajectory (or the light path) and also demonstrates the conservative features of the result. The second formulation is the standard form used in literature [18, 16, 17, 10] where cartesian position coordinates and spherical propagation direction coordinates are used. Recently, Ferwerda [16] and Khan and Jian [17] made attempts to derive such an equation. However we showed that Ferwerda’s [16] (and hence Khan and Jian’s [17]) formulations do not satisfy energy conservation and fail to account for wavefront curvature. We showed that when the refractive index is constant, our result reduces to the standard radiative transfer equation [18, 16] and when the medium is lossless and free of scattering to the well known geometrical optics equations in refractive media [19, 15, 14].