1. Introduction

Optical techniques in biomedical applications such as optical tomography and light-aided sensing of substances have been receiving tremendous interest recently [1]. These techniques of sensing substances in tissue or blood require foreign structures to be embedded/implanted in tissue in order to condition optical signals. Therefore, having a detailed understanding of how light interacts with tissue is required. Owing to its ability to accurately represent light propagation through tissue, wave propagation through biological tissue is modeled using the photon transport equation (PTE) [2, 3]. However, the interaction of electromagnetic energy with embedded objects can be best studied using Maxwell’s equations. Therefore, in order to model wave propagation through tissue with implanted foreign structures, a mapping of the PTE to Maxwell’s equations is required.

In this paper we develop a technique to map the photon transport equation to Maxwell’s equations using phase-retrieval techniques. To illustrate the applicability and accuracy of our method we analyze laser pulse propagation through a metal slit in biological tissue. Light propagation in scattering and absorbing media such as biological tissue is modeled by the photon transport equation (PTE) [2, 3] which is written in terms of the magnitude of the intensity but not the phase. Light propagation through a slit in a metal screen is described by Maxwell’s equations [4, 5], which take into consideration both the magnitude as well as the phase of the electric and magnetic fields together with the vectorial character of these fields. However, to the authors’ knowledge no work has been reported which addresses the problem of coupling these two sets of equations. Therefore, we introduce a novel strategy to couple these two sets of equations at the tissue - metal screen interface by retrieving the phase information from the intensity profile.

This paper is organized in five sections. In section 2 we introduce the formulation of the proposed technique for coupling the PTE to Maxwell’s equations using a phase-retrieval technique. In section 3 we address the problem of modeling wave propagation through a slit in a metal screen implanted in tissue. Section 4 provides the simulation results for a composite slab of tissue layer and a metal screen with a slit, a discussion on those results and possible extensions of the proposed technique. We conclude in section 5 summarizing the key features and advantages of the proposed technique.

2. Formulation: construction of phase information from a radiance profile

This section presents how we map the PTE to Maxwell’s equations by constructing the phase information using the radiance profile obtained by solving the PTE. In section 2.1 we discuss how to obtain the radiance profile by solving the PTE. In section 2.2 we present the derivation of the transport-of-intensity equation which is used for phase retrieval. In section 2.3 we show how to construct the phase profile from the radiance profile obtained by solving the PTE.

2.1. Modeling light propagation in tissue

For modeling light propagation through biological tissue, we use the photon transport equation (PTE) [3]:

where I _PTE (z,u,ϕ,t) is the radiance (units:W.m ^-2.sr ^-1.Hz ^-1), (z,θ,ϕ) are the standard spherical coordinates, u=cosθ, t is the time variable, σ _t and σ _s are attenuation and scattering coefficients, respectively, and σ _t=σ _s+σ _a, where σ_a is the absorption coefficient. The speed of light in the medium is denoted by v, P(u ^′,ϕ ^′;u,ϕ) is the phase function and F (z,u,ϕ,t) refers to the source term.

The Laguerre Runge-Kutta-Fehlberg (LRKF) method [3] can be used to solve Eq.(1) for I _PTE. The LRKF method reduces the transient PTE, given by Eq.(1), to an ordinary differential equation of only one independent variable, z, and solves it numerically. In the LRKF algorithm the discrete ordinate method [6] is used to discretize the azimuthal and zenith angles (ϕ and θ), a Laguerre expansion [3] is used to represent the time dependence and the Runge-Kutta- Fehlberg (RKF) method [7] is used to solve the one-variable ordinary differential equation in z.

2.2. Derivation of the transport-of-intensity equation for phase construction

As described in detail in [8], for a static (i.e. the electrical permittivity and the magnetic permeability are independent of time), non-magnetic (i.e. with constant permeability) medium without any current or charge densities inside, and with scatterers which slowly vary over length scales comparable to the wavelength of the incident radiation, the Maxwell equations can be reduced to [8]

where ε is the electric permittivity of the medium, µ₀ is the permeability of free space, $\nabla =\left(\frac{\partial }{\partial x}\mathbf{i}+\frac{\partial }{\partial y}\mathbf{j}+\frac{\partial }{\partial z}\mathbf{k}\right)$, t represents time and E and H denote the electric field and the magnetic field, respectively. From Eq.(2) and Eq.(3), since there is no mixing between any of the components of the electric and the magnetic field vectors, we can move on to a scalar theory [8]. Thus,

In Eq.(4) Ψ(x,y, z,t) describes the electromagnetic field and it is complex. Using the Fourier integral Ψ(x,y, z,t) can be expressed as [8]

where ω is the angular frequency. Using Eq.(5) in Eq.(4) we get

where ${\nabla }^2=\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right)$, c is the speed of light in free space and k ₀=ω/c is the wave number in free space. Then, identifying ε _ω(x,y, z)µ ₀ c ² as the square of the position-dependent refractive index, n _ω(x,y, z), of the medium, we can re-write Eq.(6) as

which is called the homogeneous Helmholtz equation [8].

In order to incorporate scattering we express Ψ _ω(x,y, z) in Eq.(7) as a perturbed plane wave [8]:

where e ^jkz represents the unscattered plane wave and Ψ̃_n (x,y,z) represents the complex envelope [9]. That is, we have considered the paraxial condition where the rays are not exactly parallel to each other; or in other words, a field with perturbed wave fronts. Using Eq.(8) in Eq.(7), we obtain

Using the identity [10]

in Eq.(9) and simplifying we get [8]

where ${\nabla }_{\mathrm{xy}}^2=\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)$. With the paraxial approximation we consider the envelope, Ψ̃_ω (x,y, z), to be “beam-like” so that its second derivative in the z direction is much smaller in magnitude than its second derivative in the x and y directions. Therefore, the $\frac{{\partial }^2}{\partial z^2}{\tilde{\psi }}_{\omega }(x,y,z)$ term in Eq.(11) can be dropped. Then Eq.(11) reduces to

Let

Here, I represents the irradiance (units:W.m ^-2.Hz ^-1) and ϕ represents the phase. Using Eq.(13) in Eq.(12) and separating the imaginary part we get the following relationship [8, 11]:

Equation (14) is called the transport-of-intensity equation (TIE). It shows how the intensity and the phase are related, and this forms the basis of the phase construction. In the next subsection we show how to retrieve the phase information from the intensity profile, by solving Eq.(14).

2.3. Construction of phase information from the irradiance profile

To construct the phase we re-write the TIE in Eq.(14) as

Equation (15) can be solved for ϕ(x,y, z) numerically using a suitable technique such as the full multigrid algorithm [12, 13], a Green-function method [11] or a fast-Fourier-transform-based method [14, 15].

Out of these techniques for solving the TIE, we adopted the full multigrid algorithm [12, 13, 16]. Equation (15) is a linear, elliptic partial differential equation of the second order and has a unique solution if I (x,y,z)>0 over a simply-connected planar region [17]. We adopt the full multigrid algorithm which solves the TIE exactly [12, 13, 16].

The full multigrid algorithm we used is briefly described below, as explained in [16]. Refer to [16] for further details.

Equation (15) can be expressed as

where $\Gamma =\left(I(x,y,z){\nabla }_{\mathrm{xy}}^2+\frac{\partial \partial I(x,y,z)}{\partial x}\frac{\partial }{\partial x}+\frac{\partial I(x,y,z)}{\partial y}\frac{\partial }{\partial y}\right),u=\varphi (x,y,z)\mathrm{and}f=-k\frac{\partial I(x,y,z)}{\partial z}.$.

In multigrid methods, we discretize the original equation on a uniform grid. We can discretize Eq.(16) as follows:

where i=1,…,M, j=1,…,M for M×M grid points. Also, I _i, _j=I(x _i,y _j, z), ϕ _i, _j=ϕ (x _i,y _j,z), Δ=x _i+1-x _i=y _j+1-y _j and $f_{i,j}=-k\frac{\partial I(x_i,y_j,z)}{\partial z}$. By solving the PTE on two closely separated planes z=z and z=z+δz, we obtain two intensity profiles. Thus, we can use the following approximations in Eq.(17).

and

Since in Eq.(18) and Eq.(19) I represents the irradiance, I and I _PTE are related by

where θ is the zenith angle used in Eq.(1) and d ω̃ is an infinitesimal solid angle [18, 19]. This conversion of the intensity is required because in the PTE we use the ray model of optics, but in Maxwell’s equations we use the wave model; and these two models deals with two different definitions of intensity, radiance and irradiance, respectively.

Since the intensity and its partial derivatives with respect to x, y and z can be approximately calculated from the two intensity profiles, as shown in Eq.(17), Eq.(18) and Eq.(19), the only unknown in Eq.(17) is ϕ _i, _j. Hence, we can use the full multigrid algorithm to solve Eq.(17) for ϕ _i, _j and thus the phase can be retrieved on each grid point.

We can discretize Eq.(16) on a uniform grid with mesh size h as

If ũ_h denotes an approximate solution to Eq.(21), then the error in ũ_h is

and the residual or the defect is

Since Γ_h is a linear operator, the error satisfies

In order to find the next approximate solution, we need to make an approximation to Γ_h in order to find v _h. Classical iteration methods, such as Jacobi or Gauss-Seidel can be used to do this.

The next approximation is generated by

Next, we form an appropriate approximation Γ_H of Γ_h on a coarser grid with mesh size H. Then the residual equation, Eq.(24), is approximated by

Since Γ_H has smaller dimension, Eq.(26) is easier to solve than Eq.(24). In the full multigrid algorithm, the first approximation is obtained by interpolating from a coarse-grid solution and at the coarsest level we start with the exact solution [16]. Using the full multigrid algorithm as detailed above we can solve Eq.(16) for u. Thus, the phase at each grid point, ϕ _i, _j is retrieved.

Before proceeding with an example that makes use of the TIE and PTE in the context of modeling wave-propagation specific effects for electromagnetic pulses traversing foreign structures implanted in biological tissue, we briefly discuss the validity and accuracy of both the TIE and the PTE in this context.

Regarding the validity and accuracy of deterministic phase retrieval using the transport-ofintensity equation, Eq.(14), we make the following three remarks: (a) The TIE has been widely employed for quantitative phase retrieval using monochromatic and polychromatic electromagnetic fields in both the visible-light and X-ray region, given a series of defocused intensity images. The TIE has also been successfully solved for the phase using matter waves such as electrons and neutrons. For a review of this work, see reference [20]. (b) Since the TIE is the continuity equation associated with the paraxial equation [11], an exact solution to which is furnished by the Fresnel diffraction integral [8], its regime of validity is restricted to paraxial beam-like fields. Interestingly, it may be used with both coherent and partially-coherent fields, a point which has been studied from both a theoretical [21] and an experimental [22] perspective. (c) Errors, in the phase retrieved using a TIE analysis, are primarily due to two sources: the finite-difference approximation to the right-hand-side of the TIE (Eq.(14)) that is given in Eq.(19), together with the presence of noise in the detected images. While the latter factor is irrelevant in the context of the analysis in the present paper, as shall become clear in the following sections, errors in the retrieved phase due to the former effect need to be considered. For an analysis of both factors, see reference [23]. The upshot of this analysis is that the error in the TIE-retrieved phase, due to a non-infinitesimal spacing δz (cf. Eq.(19)), leads to a blurring of the retrieved phase which becomes negligibly small if δz tends to zero from above. Reference [23] develops an expression for the optimal δz in the presence of a given level of noise, this optimal defocus distance being proportional to the cube root of the standard deviation of the noise. Typical reconstruction accuracies from experiments involving TIE-based phase retrieval are on the order of 1–5% [8, 22].

In this paper, we use the transport theory to model light propagation through tissue. Most of the recent advances in describing the transfer of laser energy in tissue are based on transport theory [24]. This theory is preferred in tissue optics instead of analytic approaches using the Maxwell equations because of inhomogeneities in biological tissue [24]. Several different models for numerically solving the photon transport equation have been reported in the literature and comparison of these models with the Monte-Carlo method [25] or measured data have been reported [24, 26, 27]. Thus, the application of the PTE for light propagation in biological tissue is established to be valid and accurate. The LRKF method we use for solving the PTE in this paper uses a combination of the discrete ordinates method [28] and functional expansion methods [29]. The validity and applicability of these methods in solving the PTE had been previously established [28, 29].

3. Example: wave propagation through a slit in a metal screen implanted in tissue

In this section, we apply the previously constructed theory to an object buried in tissue.

Figure 1 shows a composite object composed of a layer of biological tissue and a layer of a metal screen with a slit; and Figure 2 shows the end elevation of Figure 1. A short laser pulse is incident on the tissue layer as shown. In general, due to the index mismatch at the interface, radiation gets reflected. Even though the proposed method can easily handle such reflections at interfaces, due to the increased mathematical complexity in formulation which masks the main points of our algorithm, we limit our analysis to an index-matched surrounding at the left boundary of the tissue layer. However, the reflections at the tissue-metal screen interface and metal screen-tissue interface will be taken into consideration with a detailed formulation.

This example focuses on obtaining the magnitude and the phase of the field at z=z _A and then converting these to the corresponding electric and magnetic fields, so that the field due to the slit in the metal screen can be modeled. Then, at the exit of the metal screen, the electric and magnetic fields can be converted back to the intensity profile so that the tissue layer beyond this plane can be modeled by solving the PTE.

 

Fig. 1. (Color online) Metal screen implanted in biological tissue.

Download Full Size | PDF

 

Fig. 2. End elevation of the tissue-metal screen model.

Download Full Size | PDF

For modeling light propagation through biological tissue, (i.e. up to the tissue-metal screen interface), we use the PTE given by Eq.(1). In this paper, without loss of generality, we consider that there is no source contained inside the medium which results in F (z,u,ϕ,t)=0 in Eq. (1).

The Laguerre Runge-Kutta-Fehlberg (LRKF) method [3] can be used to solve the PTE from z=0 to z=z _A-. Thus, we can obtain the radiance profile at the plane just before the tissue-metal screen interface (i.e. at z=zA-). However, in order to model the propagation of the laser pulse beyond this plane, Maxwell’s equations should be used. Maxwell’s equations require the phase of the field in addition to the magnitude. Thus, the phase information of the field at z=z _A- should be retrieved in order to model the light propagation through the slit in the metal screen.

In order to apply the phase retrieval technique detailed above, first the radiance profile, I _PTE obtained by solving Eq.(1) should be converted to an irradiance profile I, using the relationship in Eq.(20). Thus, the irradiance at z=z _A-, I _A-, and at z=z _A--_{δ z}, I _A--_{δ z}, can be obtained by solving the PTE for radiance and integrating over the hemisphere. Then, we use the approximations in Eq.(18) and Eq.(19) in order to solve the TIE. That is,

and

The full multigrid algorithm [12, 13] is then used to solve the TIE, given by Eq.(14), for the phase, ϕ(x,y, z). Thus, the phase at the tissue-metal screen interface is retrieved using the intensity values at two infinitesimally separated planes.

Once the phase is retrieved, if the incident electric field is known, the field at the tissue-metal screen interface can be obtained. If the incident polarization vector is E ₀, as shown in Figure 1, the electric field at z=z ⁺ _A can thus be written as

Then, the corresponding magnetic field at z=z ⁺ _A can be obtained from

Thus, we have obtained the incident electric and magnetic fields at the interface. In order to calculate the field distribution just after the metal screen with the slit, we have adopted the technique introduced by Neerhoff and Mur [4, 5]. We consider a time-dependent incident profile as opposed to the time-independent profile used by Neerhoff and Mur. However, since the time variation is very slow, their technique [4, 5] can be applied for our case as outlined below. We consider TM polarization and the magnetic field is assumed to be approximately time harmonic and constant in y direction as shown in figure 3. Thus

where e _y is the unit vector in y direction. Since the time variation of U (x, z,t) is very slow, it approximately satisfies the Helmholtz equation. Hence,

where j=1,2,3 and k _j is the wave number in region j. The field in region 1 can be decomposed into three components:

where U ⁱ represents the incident field, U ^r represents the field that would be reflected if there were no slit in the screen and U ^d represents the diffracted field in region 1 due to the presence of the slit [5]. Each term on the right hand side of Eq.(33) approximately satisfies the Helmholtz equation. Also it can be shown [4, 5],

 

Fig. 3. Propagation of an incident wave through a slit in a thick metal screen.

Download Full Size | PDF

With the above set of equations and standard boundary conditions for a perfectly conducting screen, there exists a unique solution for the diffraction problem [5]. Thus, the field in region 3, close to the metal screen can be obtained using the 2-dimensional Green’s theorem as discussed in [4] and [5].

Once the electric field (E _{d 2}) just after the metal screen is obtained, these can be combined to obtain the intensity (i.e. the irradiance) using the relationship

where v and ε are the propagation speed and the permittivity in the medium, respectively. Once the irradiance profile on the plane z=d ₂ is obtained, it should be converted back to a radiance profile so that the PTE can be used to model the light propagation beyond this plane. Figure 4 shows a strategy that can be used for mapping the irradiance profile to the radiance profile, as needed for solving the PTE. In Fig. 4, the axes (x,y, z) represent the global coordinate system used in solving the PTE; also shown is the ray-centred spherical coordinate system used for describing the irradiance-to-radiance mapping forward hemisphere.

 

Fig. 4. An illustration of the strategy used for mapping radiance to irradiance.

Download Full Size | PDF

Based on the work of Ramamoorthi et al. [19], this strategy uses a hemisphere positioned centrally at the ray propagation direction and uses the relationship between radiance and irradiance given by Eq.(20), and spherical harmonic representation to achieve this task. Once the irradiance profile on the plane z=d ₂ is converted back to a radiance profile at each point, in the forward hemisphere, on the interface, the LRKF method can be used to model the light propagation through the remaining layers of tissue.

4. Numerical results and discussion

We have simulated the proposed technique using Matlab. In the simulations, all the units were normalized. The PTE, Eq.(1), is linear in intensity, I _PTE, and thus representing I _PTE as I _PTE/I ₀ does not change the equation. Therefore, we use an arbitrary scale for I _PTE throughout this paper. The time units are normalized by T _s, spatial units by Z _s and scattering and absorption coefficients by 1/Z _s.

The input pulse was taken to be a Gaussian pulse (see Fig.5) described mathematically by

where T is the factor determining the width of the input pulse while t ₀ determines the time at which pulse attains its peak value.

In this paper, we have set T _s=T. We choose this normalization factor due to the fact that the Laguerre approximation of the Gaussian pulse is very accurate for pulses with T=1 or greater. Therefore, with this scaling it is possible to obtain very accurate results even for very narrow pulses, which are used in many biomedical applications. For pulses with other shapes, it is recommended that a least square fit is used to obtain a Gaussian approximation, subsequently setting T _s to be the width of that Gaussian pulse. We have set Z _s=v×T̄. Here, T̄ can be chosen to suit the particular application. However, these scaling factors should be chosen carefully so that the matrices that are used in the LRKF method remain well-conditioned. For the simulations presented in this paper, without loss of generality, we have chosen T/T̄=1 so that Z _s=v×T.

 

Fig. 5. The incident radiance profile on the tissue layer (with arbitrary units).

Download Full Size | PDF

Figure 5 shows the irradiance profile, in a particular direction, incident on the centre of the tissue layer from the left hand side on figure 1. The proposed technique does not depend on the type of the input source. Therefore, in order to minimize the additional mathematical complexity which might mask the main idea of the proposed technique, without loss of generality, we have assumed an input source with the same radiance profile in all directions in the forward hemisphere, as depicted in Fig. 5. However, the proposed technique can be applied to other kinds of input sources; for example, one may use the afore-mentioned technique proposed by Ramamoorthi et al. [19], to construct an input radiance profile with non-uniform profile.

 

Fig. 6. (Color online) The irradiance profile on a plane just before the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

 

Fig. 7. (Color online) The electric field distribution on a plane just before the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

Figure 6 shows the irradiance profile at z=2, on a plane just before the tissue-metal screen interface, obtained by solving the PTE using the afore-mentionedLRKF method. Here, we have modeled the tissue layer with a normalized scattering coefficient of 0.3, normalized absorption coefficient of 0.5 and the Henyey-Greenstein phase function [18] with an asymmetry factor of 0.7. The normalized velocity was taken to be 1 while the refractive index of the tissue layer was assumed to be 1.37. Figure 6 shows how the irradiance profile on the (x,y) grid at z=2 varies with time.

 

Fig. 8. (Color online) The magnetic field distribution on a plane just before the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

 

Fig. 9. (Color online) The magnetic field distribution on a plane just after the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

Using the same technique, the irradiance profile at z=1.95 was obtained, and these two profiles were used to retrieve the phase of the field at z=2. For phase retrieval, we first translated the code given in [16] for the full multigrid algorithm to Matlab scripting, and then modified it to solve the TIE, which involved slight modifications to a few subroutines. Then, the irradiance values and the phase values were combined according to Eq.(29) to construct the electric field at z=2, which is shown in figure 7.

 

Fig. 10. (Color online) The electric field component in the x-direction on a plane just after the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

 

Fig. 11. (Color online) The electric field component in the z-direction on a plane just after the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

Equation (30) was used to calculate the magnetic field distribution on a plane just before the metal screen, using the electric field distribution. This result is shown in figure 8. Then, the field distribution on a plane just after the screen was obtained using the technique introduced by Neerhoff and Mur [4], as discussed in the previous section. The magnetic and electric fields thus obtained are shown in figures 9, 10 and 11. The irradiance profile constructed according to Eq.(36) is shown in figure 12.

The proposed technique can be extended to model more complicated structures, such as photonic crystals, implanted in biological tissue.

 

Fig. 12. (Color online) The irradiance profile on a plane just after the tissue-metal screen interface (with arbitrary units).

Download Full Size | PDF

5. Conclusion

This paper introduced a novel strategy by which Maxwell’s equations and the photon transport equation can be seamlessly integrated to analyze electromagnetic radiation in tissue-like media. The need for such analysis stems from the perceived biomedical applications in sensing and characterizing tissue properties. By using the present technique it is possible to analyze diffraction effects within the frame work of radiative transfer theory.

The proposed technique can be used to assist the development of biomedical instruments that can be used for noninvasive diagnosis of diseases. Moreover, it enables us to calculate the light energy distribution within tissue structures with surgically implanted foreign structures. Estimating the light energy distribution in tissue is essential for dosimetry for either photodynamic or for thermal coagulation therapy [30]. Laser treatment is very popular at present and techniques developed for analyzing reflected or transmitted light from tissue can be used as a diagnostic tool or for evaluation of the progress of laser treatment [30].

The maximum safe exposure of laser light for the skin is 0.1 Joule/cm ² per pulse or 1.0 Watt/cm² for continuous exposure [31]. However, structures such as photonic crystals can be implanted to obtain enhanced signals by properly engineering the photon density of states. The proposed technique can be used to model such foreign structures implanted in tissue. Since our main focus is on how to map the photon transport equation to Maxwell equations, for our simulations we have used a simple foreign structure, a metal screen with a slit, in order to introduce the proposed technique without the additional mathematical complexity of modeling more complicated structures.

In modeling wave propagation through biological tissue with foreign structures implanted, the PTE models wave propagation through the tissue layer. At the interface, the phase is retrieved from the irradiance profile and thus the electromagnetic field is determined. The wave propagation through the foreign structure is modeled using Maxwell’s equations. Then, the electromagnetic field is converted back to an irradiance profile at the exit of the foreign structure.