1. Introduction

Many models have been developed to solve the one-dimensional photon transport equation (PTE), which assumes horizontally uniform plane-parallel media. The modeling of the three-dimensional photon transport is usually considered difficult, even for the steady state case, because it involves solving integro-differential equations of four or five variables [1]. However, many diverse methods for solving the multidimensional steady state photon transport problem, which discard the time-dependency complications of the transient problem considered here, have been found in literature [2, 3].

The use of short pulse lasers in engineering and biomedical applications, such as optical tomography and laser ablation has motivated research in transient photon transport through biological media [1]. Transient signals have the possibility to probe specific resonances in biological media, which far exceed the capabilities rendered by steady state methods. Most of the existing transient models had been developed for the one-dimensional PTE; but there are a few models developed for the two-dimensional PTE [4]. Only recently have researchers started to work on models for the three-dimensional transient photon transport in scattering media such as biological tissue. Also, there has been renewed interest in having sensors implanted in tissue for monitoring and characterization [5–10].

Many three-dimensional photon transport problems have been modeled by reducing the 3D propagation to equivalent 1D models. However, such mapping inevitably discard the inhomogeneities and anisotropic characteristics of such media. For example, realistic analysis of light propagation in biological tissue must be carried in a framework which accounts for both weak anisotropy and inhomogeneity. Moreover, if sensors are implanted in such scattering media, 1D models in general has difficulty accounting for both scattering and diffraction effects. Also, specific surface features of a model problem may be more accurately incorporated into threedimensional photon transport. Owing to these reasons, 3D models can provide quantitatively superior realistic results compared with its one or two dimensional counterparts.

Most of the recent developments on solving the three-dimensional transient PTE rely on approximate methods such as the diffusion approximation and spherical harmonics approximation, which are used to simplify the PTE [4]. These approximate methods have failed to predict the correct propagation speed within the medium and also the solution accuracy is not satisfactory except for specific cases that consider optically thick media at asymptotically longer time scales [4]. Tan et al. [4] have proposed an integral equation technique to model the transient radiative transfer in 3D homogeneous as well as non-homogeneous participating media. Guo et al. [1] have formulated a complete transient three-dimensional discrete ordinates method to solve the transient PTE in a nonhomogeneous turbid medium with a rectangular enclosure. Guo et al. [11] have modeled heterogeneous biological tissues, using the discrete ordinates method, incorporating Fresnel specularly reflecting boundary condition. They treated the laser intensity as having a diffuse part and a specular part. The reflectivity at the tissue-air interface was calculated by using the Snell’s law and the Fresnel equations [11]. Chai et al. [12] have proposed a finite-volume method to calculate transient radiative transfer in a three-dimensional enclosure. Sawetprawichkul et al. [13] have implemented the Monte Carlo method for modeling the three-dimensional transient photon transport in a parallel computer system.

Oliveira et al. [14] have proposed a semi-analytical numerical method for solving the one-dimensional transient PTE in slab geometry. They represented the time dependence of the radiance and its time derivative using a truncated series of Laguerre polynomials. This results in a set of time-independent equations which are then solved recursively using a Laplace transform method in the space variable with analytical inversion.

In this paper, we propose a novel solution method for solving the three-dimensional transient PTE, which can also be easily adopted for solving one and two dimensional cases. This is an extension of the model we developed for the one-dimensional transient case in [15]. Current work uses a transformation which eliminates the partial derivative term with respect to time. We then expand the radiance using a Laguerre basis which transforms the transient PTE to the steady state version, which is then solved using a finite volume approach with the discrete ordinates method (DOM) [16]. With the use of a Laguerre basis we can represent any arbitrary pulse shape using only a few Laguerre polynomials, considerably reducing subsequent computational overhead. Also, the causality of the system is implicitly imposed by the Laguerre expansion. The use of the discrete ordinates method has several advantages. It requires only a single formulation to invoke higher-order approximations of discrete ordinate quadrature sets; and also, the DOM is applicable to the complete anisotropic scattering phase function and inhomogeneous media [1].

This paper is organized as follows. In Section 2, we provide a detailed mathematical description about the proposed model. This section is further broken down into three subsections. The first subsection carries a step by step derivation of the proposed algorithm for the more general three-dimensional case. In the second and third subsections, descriptions on how the model reduces to two and one dimensional geometries are provided. In Section 3, we present results obtained by numerically implementing the proposed algorithm. First, a validation of the model is provided using the one-dimensional case, with special parameters for which the results can be intuitively derived. Then, a comparison of the results obtained for all the three dimensions is provided against the results obtained using the algorithm proposed by Guo et al. in [1]. We conclude this paper in Section 4 after highlighting relative advantages of the proposed method.

2. Formulation

This section introduces the proposed solution technique. In the first sub-section, we present a detailed mathematical description of the method, as applied for the three-dimensional PTE. In the second sub-section, we describe how this method can be adopted to solve the two-dimensional PTE. In the last section, a description on how to adopt the proposed method for the one-dimensional case is presented.

2.1. For three-dimensional PTE

The three dimensional transient photon transport equation (PTE) is given by [1]

where I(x,y, z, s, t) is the light intensity (radiance) with unitsW.m ^-2.sr ^-1.Hz ^-1, (x,y, z) are the cartesian coordinates, s is the local solid angle, µ, ξ and η are direction cosines such that µ=cosθ, ξ=sinθ cosϕ, η=sinθ sinϕ and t is time; σ_t (x,y,z) and σ_s(x,y,z) represent the position dependent attenuation and scattering coefficients, respectively. The speed of light in the medium is v, P(s′; s) is the phase function and F(x,y, z, s, t) is the source term. Since we consider a constant refractive index, the speed of light in the medium is a constant. In this paper we assume that there are no internal sources inside the medium and hence F=0. Figure 1 shows the incident direction, s, the scattered direction, s′, and the scattering angle. Figure 2 shows the relationship between the Cartesian and spherical coordinates, showing the zenith and azimuthal angles.

 

Fig. 1. Illustration of the incident direction, scattered direction and the scattering angle.

Download Full Size | PDF

 

Fig. 2. Relationship between Cartesian and spherical coordinates.

Download Full Size | PDF

Before proceeding to solve the PTE, the following substitution is used to map Eq. (1) to a moving reference frame with the pulse:

Using Eq. (2) in Eq. (1), we are able to eliminate the partial derivative term with respect to time in the original PTE, as shown below.

Equation (3) can be simplified to get,

Then the time dependency is represented using Laguerre polynomials. We use Laguerre functions to expand I(x,y, z, s,τ) in a summation.

Laguerre polynomials are the canonical solutions of the differential equation

which is also known as the Laguerre’s equation [17]. These solutions are a sequence of orthogonal polynomials which can be generated by Rodrigues’ formula [17]

These polynomials satisfy the orthogonality property

where δ_nm is the Kronecker delta [18]. Laguerre polynomials are causal [19] and thus the causality of the system is implicitly imposed. Also, any causal function can be expanded using a Laguerre basis because it forms a complete orthogonal basis in real space [20]. Since such expansions are a linear superposition of causal Laguerre polynomials, the causality property of the original system is implicitly retained.

We expand the radiance, I, as a summation of N Laguerre polynomials;

where ∫^∞ ₀ L_n (τ)L_m (τ)e ^-τ dτ=δ_mn. Using Eq. (5) in Eq. (4) and taking moments (i.e. multiplying by L_n (τ)e^-τ and integrating over [0,∞)), the time dependence can be removed resulting in N uncoupled equations, one for each Laguerre coefficient, B_n.

where n=1, …,N. Thus, the transient photon transport equation has been reduced to a set of uncoupled steady state photon transport equations. Equation (6) can be solved using the discrete ordinates method as detailed in [16], and outlined below.

The solid angle dependence of Eq. (6) can be discretized using the discrete ordinates method [21], resulting in L coupled equations for each Laguerre coefficient.

where w_j are quadrature weights, i=1, …,L and n=1, …,N. We divide the tissue specimen into infinitesimally small control volumes and apply the discrete ordinates method. We integrate Eq. (7) over a control volume:

where V=ΔxΔyΔz is the volume of the control volume, (B_n)_xu and (B_n)_xd are average values of B_n over the faces A_xu and A_xd of the control volume, respectively. Similarly,

Here, Bⁱ_n=B_n(x,y, z, sⁱ). Using Eq. (8), Eq. (9) and Eq. (10) in Eq. (7),

where (Bⁱ_n)p is the Laguerre coefficient of the radiance at the centre of the control volume. The average values of B_n on the faces of the control volume and that at the centre are related by [16],

The weights γ_x, γ_y and γ_z can be determined using the scheme proposed by Lathrop [1, 22]:

$\begin{array}{c}{\left({\gamma }^i\right)}_x=max(0.5,1-\frac{\mid {\xi }^i\mid /\Delta x}{2\left(\mid {\eta }^i\mid /\Delta y+\mid {\mu }^i\mid /\Delta z\right)+{\sigma }_t}),\\ {\left({\gamma }^i\right)}_y=max(0.5,1-\frac{\mid {\eta }^i\mid /\mathrm{\Delta y}}{2\left(\mid {\xi }^i\mid /\Delta x+\mid {\mu }^i\mid /\Delta z\right)+{\sigma }_t}),\\ {\left({\gamma }^i\right)}_z=max(0.5,1-\frac{\mid {\mu }^i\mid /\mathrm{\Delta z}}{2\left(\mid {\xi }^i\mid /\Delta x+\mid {\eta }^i\mid /\mathrm{\Delta y}\right)+{\sigma }_t}),\end{array}$.

Using Eq. (12), Eq. (13) and Eq. (14) in Eq. (11) and eliminating (Bⁱ_n)_xd, (Bⁱ_n)_yd and (Bⁱ_n)_zd,

where

Dividing the numerator and the denominator of the right hand side of Eq. (15) by V, we get

In this work, the spatial differencing scheme suggested by Fiveland [23] is used to determine the step sizes Δx, Δy and Δz. According to this scheme

$\begin{array}{c}{\Delta }_x<\frac{{\mid {\xi }^i\mid }_{min}}{{\sigma }_t\left(1-{\left({\gamma }^i\right)}_x\right)},\\ {\Delta }_y<\frac{{\mid {\eta }^i\mid }_{min}}{{\sigma }_t\left(1-{\left({\gamma }^i\right)}_y\right)},\\ {\Delta }_z<\frac{{\mid {\mu }^i\mid }_{min}}{{\sigma }_t\left(1-{\left({\gamma }^i\right)}_z\right)}.\end{array}$.

The algorithm of discrete ordinates for solving the PTE, which is reduced to the steady state version, given by Eq. (7) starts at corner control volumes. The (Bⁱ_n)_xu, (Bⁱ_n)_yu and (Bⁱ_n)_zu values for each corner control volume are given by the Laguerre coefficients of the boundary condition. Hence, the centre values, (Bⁱ_n)p, can be determined from Eq. (17). Once (Bⁱ_n)_p is calculated (Bⁱ_n)_xd, (Bⁱ_n)_yd and (Bⁱ_n)_zd can be determined using Eq. (12) to Eq. (14). Also, (Sⁱ)_p, which is initially assumed to be equal to zero, is updated using the calculated (Bⁱ_n)_p value from Eq. (16). (Bⁱ_n)_xd, (Bⁱ_n)_yd and (Bⁱ_n)_zd of the current control volume are equal to (Bⁱ_n)_xu, (Bⁱ_n)_yu and (Bⁱ_n)_zu of the adjacent control volumes, respectively. A more detailed description on this method can be found in [16].

The overall algorithm proposed in this paper is summarized in ALGORITHM 1 below:

ALGORITHM 1

1. Use the transformation given by Eq. (2) in the transient PTE, Eq. (1).

2. Expand the radiance, I(x,y, z, s, t), using a Laguerre basis.

3. Expand the time-dependent boundary condition using a Laguerre basis.

4. Divide the specimen geometry into infinitesimally small control volumes (finite volume approach).

5. Apply the discrete ordinates method to Eq. (6) in order to solve for each Laguerre coefficient of the radiance, iteratively until convergence is achieved.

6. Combine the Laguerre coefficients to obtain the radiance values.

The main steps of the discrete ordinate method using the finite volume approach are listed in ALGORITHM 2 below:

1. Start at a corner control volume. (Bⁱ_n)_xu, (Bⁱ_n)_yu and (Bⁱ_n)_zu for this control volume are known from the boundary condition.

2. (Sⁱ)p is assumed to be zero initially for all control volumes.

3. Calculate (Bⁱ_n)p of the current control volume from Eq. (17).

4. Determine (Bⁱ_n)_xd, (Bⁱ_n)_yd and (Bⁱ_n)_zd using Eq. (12) to Eq. (14).

5. The quantities calculated in step 4 become (Bⁱ_n)_xu, (Bⁱ_n)_yu and (Bⁱ_n)_zu for adjacent control volumes.

6. Repeat steps 3 to 5 until all control volumes are traversed.

7. Move to the next discrete ordinate (of the solid angle, s) in the currently considered quadrant.

8. Repeat steps 2 to 7 until all discrete ordinates corresponding to the current quadrant are covered.

9. Move to the next corner control volume.

10. Repeat steps 1 to 9 until all the corners are considered.

11. Update (Sⁱ)_p for all discrete ordinates using Eq. (16).

12. Repeat steps 1 to 11 until (Sⁱ)p is converged.

13. Repeat steps 1 to 12 for all Laguerre coefficients of I.

2.2. For two-dimensional PTE

The two dimensional transient photon transport equation (PTE), for a medium without any internal sources, is given by

Instead of the transformation given by Eq. (2), we have to use the transformation

in order to make the partial derivative term with respect to time of Eq. (18) vanish. After using Eq. (2) in Eq. (18) to eliminate the time derivative term, the radiance, I, is expanded using a Laguerre basis using Eq. (5) and take moments as done for the three-dimensional case. This results in the two-dimensional steady state PTE:

Then, ALGORITHM 2 can be used Eq. (20) to solve for B_n and using Eq. (5) the radiance, I, can be determined. However, when applying ALGORITHM 2, in the equations used for the three-dimensional case, the terms corresponding to the variable y vanish.

2.3. For one-dimensional PTE

The one dimensional transient photon transport equation (PTE), for a medium without any internal sources, is given by

For this case, instead of the transformation given by Eq. (2), we have to use the transformation

Use of Eq. (22) in Eq. (21) eliminates the time derivative term in Eq. (21). Then, expanding the transformed one-dimensional PTE using a Laguerre basis and taking moments we get,

where n=1, …,N. Then we discretize the solid angle, s, of Eq. (23) using the discrete ordinates method, which results in

where i=1, …,L and n=1, …,N. Thus, there are L coupled equations for each Laguerre coefficient, B_n and Eq. (24) corresponds to the i ^th discrete ordinate of s. We can write this set of equations in matrix form as;

where B _n=[B_n(z, sⁱ] _L,1, P=[P(s ^j;s ⁱ)] _L,L, and, A is a L by L diagonal matrix with diagonal elements µⁱ. The matrix W is also a L by L diagonal matrix with diagonal elements w_j.

Rearranging, Eq. (25) can be written as

where $\mathbf{Y}={\mathbf{A}}^{-1}\left[\frac{{\sigma }_s}{4\pi }\mathbf{PW}-{\sigma }_t\mathbf{I}\right].$. Hence, the original transient PTE is reduced to a one-variable ordinary differential equation. The boundary condition should also be expanded using a Laguerre basis as described for the three-dimensional case. Thus, Eq. (26), subject to the given boundary condition, can be solved using the 4^th order Runge-Kutta-Fehlberg(RKF) method [24]. A more detailed description about solving the one-dimensional PTE using this technique (the Laguerre Runge-Kutta-Fehlberg method) can be found in our paper [15]. In that paper, we have represented the solid angle using the zenith and azimuthal angles. However, the technique is essentially the same, the only difference being the matrix sizes.

3. Numerical results and discussion

We have simulated the proposed technique using Matlab^™. First, we validate the proposed technique for special cases, for which the results are intuitive and readily computable using the one-dimensional PTE. Then we provide a comparison for all the three cases; one, two and three dimensions, with results obtained from the technique proposed by Guo et al. [1]. In [1] they have provided a complete transient three-dimensional discrete ordinates method to solve the three-dimensional PTE in a rectangular enclosure. In that paper, they have extended the discrete ordinates algorithm described above in ALGORITHM 2, to include the transient term. We adopted their method for the two and one dimensional cases as well for the purpose of comparison. In the figures, we assigned the label “Laguerre DOM” for our method and “Transient DOM” for the method proposed in [1]. In the simulations, all the units were normalized. A brief account on the normalization used is provided below.

 

Fig. 3. Incident Gaussian pulse (with I ₀=1 in Eq. (27)).

Download Full Size | PDF

3.1. Normalization of units

The PTE is linear in intensity, I, and thus arbitrary units for intensity can be used. Also, it is a first order, linear differential equation in time and thus time units can also be selected to enhance the numerical accuracy. It is very clear that input pulse parameters influence the numerical accuracy and efficiency of the algorithm. The input pulse was taken to be a Gaussian pulse 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. Thus, we set the intensity and time units relative to major characteristics of the input pulse; its peak intensity, I ₀ and its temporal width controlling parameter, T, which is related to full-width half maximum value by $2\sqrt{2\ln \left(2\right)}T.$. We use the scale I/I ₀ for radiance throughout this paper. The time units are normalized by a scaling factor T_s, spatial units by Z_s and scattering and absorption coefficients by 1/Z_s.

 

Fig. 4. (Color online) Short laser pulse incident on the biological tissue layer

Download Full Size | PDF

In this paper, we have set T_s=T/1.5. 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.5 or greater. This effect is shown in Fig. 5 and Fig. 6. Figure 5 presents a comparison of direct Laguerre approximation (i.e. without scaling with T_s) with the exact plot of the Gaussian pulse, for four different pulse widths. This shows that Laguerre approximation for pulses with T=1 and T=1.5 are very accurate. However, for the T=1 case, a very small bump appears for 6.5<x<7.5 (see Fig. 5b). Such inaccuracies can be always eliminated by choosing larger number of Laguerre coefficients. However, if the computational efficiency is also of concern, then better scaling can be used to reduce the required number of Laguerre coefficients. Therefore, we have chosen the T=1.5 and t ₀=4 approximation. A pulse with different T and t ₀ values should be scaled and shifted accordingly to obtain T=1.5 and t ₀=4 before the Laguerre fit is carried out. The results obtained by the simulations should be scaled and shifted back (reversing this process). This idea is illustrated in Fig. 6. Four different Gaussian pulses with different T and t ₀ values are scaled and shifted (up or down) to obtain the Gaussian pulse with T=1.5 and t ₀=4 and a Laguerre fit is carried out. Then the scaling and shifting is reversed and the resulting fit is plotted with the exact plot of the corresponding pulse. It is clearly shown in Fig. 6 that this scaled Laguerre fit produces very accurate approximations.

 

Fig. 5. Comparison of Laguerre fitting, without scaling, with the exact plot of the corresponding Gaussian function.

Download Full Size | PDF

It is interesting to note that with this scaling, it is possible to obtain very accurate results even for very narrow pulses, that are used in many biomedical applications. For pulses with other shapes, it is recommended that a least square error fit is used to obtain a Gaussian approximation, subsequently setting T_s to the width of that Gaussian pulse divided by 1.5. With the Laguerre fit it is not possible to approximate any pulse shape to an arbitrary accuracy because of numerical considerations. Therefore, in practice, there is a finite domain in which this approximation is valid. However, as it can be seen in Fig. 6 this domain expands as T is increased. However, since in this algorithm the Laguerre polynomials are propagated with the pulse, the approximation will always be accurate in the domain we are interested in (i.e. the domain in which the Laguerre fit is accurate extends few time units beyond the required observation point). A discussion on the observation window in which the Laguerre approximation is accurate can be found in [15].

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 remain well-conditioned. For the simulations presented in this paper, without loss of generality, we have chosen T/T̄=1.5 so that $Z_s=v\times \frac{T}{1.5}$. In this paper we use a line over the letters to denote quantities with normalized units.

 

Fig. 6. Comparison of Laguerre fitting, after scaling with the factor T_s, with the exact plot of the corresponding Gaussian function.

Download Full Size | PDF

3.2. Computational complexity

For the simulations we have considered an observation volume of -1<x̄<+1,-1<ȳ<+1 and 0<z̄<2. However, we have assumed refractive index-matched boundaries in all the three dimensions. Therefore, in order to minimize effects from internal reflections at the boundaries the simulation volume was taken to be larger than the observation volume in all the three dimensions (-5<x̄<+5,-5<ȳ<+5 and 0<z̄<5). The total number of voxels used in the finite volume computations was 108000. The number of voxels considered in the observation volume was 1728. The number of discrete ordinates taken for the simulations was 80. We have used the discrete ordinates for S ₈ approximation given in Table 16.1 of reference [16]. The two main drawbacks of the discrete ordinates method are false scattering and the ray effect [4, 16]. False scattering is due to the spatial discretization errors. When a single collimated beam is traced through an enclosure by the discrete ordinates method, the beam will gradually widen as it moves farther away from its point of origin; and this unphysical effect, even in the absence of real scattering, is called false scattering [16]. False scattering can be reduced by using finer meshes [16]. The ray effect is due to the errors in angular discretization and can be reduced by increasing the sizes of the meshes [16]. Therefore, if a finer spatial mesh is used to reduce the false scattering, a finer angular quadrature scheme should be used to reduce the ray effect [16]. Thus, if one uses the Laguerre DOM method proposed in this paper with a finer mesh, it should be accompanied by an increased number of discrete ordinates in order to minimize numerical errors.

The Gaussian pulse shape can be represented accurately using about 40 Laguerre polynomials. Increasing the number of polynomials up to 63 enables further refinement in the results. The incident Gaussian pulse in Fig. 3 can be approximated by 20 Laguerre polynomials with a normalized root mean squared deviation (NRMSD) of 8.07% in the observation window. The NRMSD can be reduced to 1.5% by using 40 polynomials. A further reduction up to 0.39% can be achieved by using 63 polynomials. However, increasing this number further does not improve the approximation because machine accuracy hinder further numerical improvements.

The simulation presented in this section was carried out in a computer with 2 Intel CPUs at 2.5 GHz. For the three-dimensional simulation of the proposed method (Laguerre DOM) shown in Fig. 13 to Fig. 16, with 63 Laguerre polynomials, it took 9 minutes and 1 second. For simulating the same problem using Transient DOM the same machine took 34 minutes and 6 seconds. In [1] the authors state that their method may take several minutes to dozens of hours for simulations depending on the specified problem and the grid size. The efficiency of the proposed Laguerre DOM might be even better for more complex problems than the one considered here.

Figure 4 shows a short laser pulse of Gaussian shape incident on a biological tissue layer at (x̄=0,ȳ=0, z̄=0). This is the model used for the simulations presented in this paper. For all the simulations, we have used T̄=1.5 and $\overline{t_0}=4$ for the input pulse given by Eq. (27), the normalized velocity to be v̄=1 and the normalized thickness of the tissue layer, at which the results were obtained, to be z̄=2. The isotropic phase function (P(s′;s)=1) was used in the simulations. However, the proposed method is able to handle any other phase function.

 

Fig. 7. Variation of radiance with time, along the incident angle, at z̄=2 without any scattering or absorption.

Download Full Size | PDF

Figures 7, 8 and 9 are used for validating the proposed technique for the one dimensional case. These figures were obtained for special cases for which the results are intuitive. For Fig. 7, the tissue medium was assumed to be lossless, without any scattering or absorption. That is, σ_s and σ_t in the PTE are set to zero. With this setup, there should be neither a decay in the intensity, nor scattering to other directions. Figure 7 shows the variation of radiance with time, along the direction of incidence (ie. µ=1), at z̄=2. As expected, at z̄=2, we have obtained the same Gaussian pulse, without any loss, but with the corresponding time delay.

 

Fig. 8. (Color online) Variation of radiance with time and µ, at z̄=2 with absorption but no scattering.

Download Full Size | PDF

For Fig. 8, the medium was assumed to have a normalized absorption coefficient, $\overline{{\sigma }_a}$, of 0.2 but without any scattering (hence, the normalized attenuation coefficient, $\overline{{\sigma }_t}=0.2$ as well). The direction of incidence for this simulation was taken to be at µ=0.5. Figure 8 shows the variation of radiance with time and the cosine of the zenith angle, µ, at z=2. As expected, a decayed radiance profile exists along the incident direction, µ=0.5, but since there was no scattering no radiance values are seen along other directions.

 

Fig. 9. (Color online) Variation of radiance with time and µ, at z̄=2 with scattering and absorption.

Download Full Size | PDF

For Fig. 9, the medium was assumed to have a normalized absorption coefficient, $\overline{{\sigma }_a}$, of 0.02 and a normalized scattering coefficient, $\overline{{\sigma }_s}$, of 0.98 (hence, $\overline{{\sigma }_t}=1$). The direction of incidence for this simulation was taken to be at µ=0.17. Figure 9 shows the variation of radiance with time and µ. In this figure, the scattering to other directions is exhibited clearly.

Figures 10 to 18 show comparisons of the proposed method with the method proposed in [1], for one, two and three dimensional cases, respectively. In figures 7 to 9, we have shown the radiance profile so that it is possible to explain the known results when there is no attenuation, when only absorption is present, and when both absorption and scattering is present. However, for the comparisons we have obtained irradiance profiles, which have more physical meaning, specially in practical measurements. For figures 10 to 18, $\overline{{\sigma }_t}=1$, $\overline{{\sigma }_s}=0.98$ and the direction of incidence is µ=1.

 

Fig. 10. (Color online) Comparison: Variation of irradiance with time for 1D PTE.

Download Full Size | PDF

Figure 10 shows the variation of irradiance with time at z̄=2, for the one-dimensional PTE using the Laguerre DOM and the Transient DOM.

 

Fig. 11. (Color online) Variation of irradiance with time and x coordinate at z̄=2 for 2D PTE, using Laguerre DOM.

Download Full Size | PDF

Figures 11 and 12 show the variation of irradiance with time and along x-axis at z̄=2, for the two-dimensional PTE using the two methods. In Fig. 11 one can observe that the time shift of the peak value of the irradiance profile increases as x̄ moves away from x̄=0. The reason for this is that light takes the shortest time to reach (x̄=0, z̄=2) point and the time increases as the point moves away from x̄=0 central axis. However, in the simulation obtained using Transient DOM shown in Fig. 12 this physical phenomenon is not clearly visible.

Figures 13 to 16 show the variation of intensity with time, along x and y axes, at z̄=2, for the three-dimensional PTE using Laguerre DOM. In these figures the value of the normalized irradiance (ranging from 0 to 16×10^-4) is represented by the color scale. Figures 13 to 16 were obtained using the same data set. Figure 13 shows a slice plane at ȳ=0. Thus, in Fig. 13 the variation of irradiance with time along the x-axis on (ȳ=0, z̄=2) plane is shown. As expected, the light reaches (x̄=0,ȳ=0) point first as rays traveling at an angle to the central axis (x̄=0,ȳ=0) take longer to reach z̄=2 plane (because the diagonal distance is more than the central axis distance and the speed is constant). For the same reason the irradiance profile decays first at (x̄=0,ȳ=0) and later at other points. The more we move away from the central axis, the longer it takes the irradiance profile to decay.

 

Fig. 12. (Color online) Variation of irradiance with time and x coordinate at z̄=2 for 2D PTE, using Transient DOM.

Download Full Size | PDF

 

Fig. 13. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM (Slice plane at ȳ=0, z̄=2).

Download Full Size | PDF

Figure 14 shows the variation of irradiance along x and y axes on z̄=2 plane at t̄=6. With isotropic scattering and normal incidence we have obtained a circular spatial distribution on the xy plane. From Fig. 13 it can be seen that at t̄=6 points around (x̄=0,ȳ=0) are receiving the second half of the gaussian pulse (i.e. the maximum irradiance has been reached earlier) and points further away from the central axis are receiving light corresponding to the maximum or the first half of the incident pulse. This phenomenon is reflected in Fig. 14. However, in the simulation result shown in Fig. 13 and Fig. 14, an uncharacteristically low number of photons can be seen along the central axis and the majority of diffuse photons appear off axis resulting in a well defined shape. The reason for this phenomenon is not very clear yet, and we intend to investigate this issue using Monte-Carlo techniques in the future. Figures 15 and 16 show that the irradiance profile at z̄=2 is symmetrical along x and y axes due to isotropic scattering and normal incidence.

 

Fig. 14. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM (Slice plane at z̄=2, t̄=6).

Download Full Size | PDF

 

Fig. 15. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM.

Download Full Size | PDF

Figures 17 and 18 show the variation of intensity with time, along x and y axes, at z̄=2, for the three-dimensional PTE using Transient DOM proposed in reference [1]. The temporal profile in Fig. 17 obtained using this method is different to that in Fig. 13. Since we use a Gaussian shaped input, the irradiance profile should start decaying at (x̄=0,ȳ=0, z̄=2) before at other points. The further away the observation point moves from the central axis, the longer it should take for the decaying to start. This difference in Fig. 13 and Fig. 17 may be due to the strong numerical diffusion and false propagation in Transient DOM of Guo and Kumar [1]. In [1] the authors state that an obvious disadvantage of transient analysis using their method is that the numerical diffusion and false propagation are inevitable. In order to minimize the numerical diffusion, the spatial grid and the time step are required to be as fine as possible. The size of time step affects the transient behavior. Figure 3 and Fig. 4 in [1] show influences of time step and spatial grid on the temporal behavior, with a comparison to Monte Carlo simulations. Figure 3 clearly shows that the Transient DOM cannot capture the abrupt rise of transmittance in the early time period. In Laguerre DOM, however, there is no discretization in timescale. We represent the entire time domain of the input pulse as a truncated series of Laguerre polynomials and then solve the PTE for each Laguerre coefficient separately. Also, using the substitution in Eq. (2) we map the PTE to a moving reference frame with the pulse.

 

Fig. 16. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Laguerre DOM.

Download Full Size | PDF

 

Fig. 17. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Transient DOM (Slice plane at ȳ=0, z̄=2).

Download Full Size | PDF

 

Fig. 18. (Color online) Variation of irradiance with time, x coordinate and y coordinate at z̄=2 for 3D PTE, using Transient DOM (Slice plane at z̄=2, t̄=6).

Download Full Size | PDF

In [1], the authors have used Duhamel’s superposition theorem, which can only be applied to linear systems, to incorporate the time-dependent boundary conditions. This theorem relates the solution of the problem subject to the time-dependent boundary condition to the solution subject to time-independent unit step input. The Laguerre DOM represents the time-dependency of the boundary condition using a Laguerre expansion. Then, a set of time-independent equations corresponding to each Laguerre coefficient is solved separately using DOM. Thus, one possible explanation for the difference in temporal and spatial irradiance profiles might be due to the strong numerical diffusion and false propagation in Transient DOM of Guo and Kumar [1].

4. Conclusion

This paper introduced a novel technique for modeling multidimensional transient photon transport, for applications in bio-sensing and in short pulse propagation through turbid media. The proposed method uses a transformation to eliminate the transient term in the transient PTE. It then expands the radiance using a Laguerre basis to automatically account for causality while providing an efficient basis suitable for calculations. This reduces the original transient PTE to a steady state version. Hence, the discrete ordinates method, using a finite volume approach can be used to solve for the radiance.

Since the time dependence is expanded using a Laguerre basis all the sampling points in the time domain is obtained in a single execution, as opposed to time marching techniques used in existing solution methods. This makes the proposed algorithm much faster when one requires the intensity profile at a particular point or a plane over a time interval. In addition, the use of the Laguerre expansion to represent the time dependency enables modeling the system with any arbitrary input pulse shape, using only a few Laguerre polynomials. Specifically, the Gaussian pulse shape, which is used in many practical applications can be accurately represented using a few Laguerre polynomials, as opposed to the discrete sampling used in most of the existing models. Also, this expansion implicitly imposes the causality of the system.