The time-dependent one-dimensional photon transport (radiative transfer) equation is widely used to model light propagation through turbid media with a slab geometry, in a vast number of disciplines. Several numerical and semi-analytical techniques are available to accurately solve this equation. In this work we propose a novel efficient solution technique based on eigen decomposition of the vectorized version of the photon transport equation. Using clever transformations, the four variable integro-differential equation is reduced to a set of first order ordinary differential equations using a combination of a spectral method and the discrete ordinates method. An eigen decomposition approach is then utilized to obtain the closed-form solution of this reduced set of ordinary differential equations.
© 2011 Optical Society of America
The photon transport equation (PTE), also known as the radiative transfer equation (RTE), is used in various different disciplines to model light propagation through turbid media and to characterize the scattering properties of media. For example, biomedical applications such as near-infrared fluorescence tomography , astrophysical and cosmological applications such as problems related to the reionization of the intergalactic medium and absorption line signatures of high redshift structures , meteorological applications such as visibility studies (e.g. modeling haze) and modeling shortwave and longwave fluxes transmitted by fields of broken clouds , all require the solution to the radiative transfer equation.
A number of models for solving the steady state (time-independent) PTE has been developed over the past four decades [4, 5], ignoring the time dependency of the intensity profile. Some of these models include the discrete ordinates method , integral transformation techniques and the FN method . Siewert  and Larsen et al.  have worked on the inverse-source problem where the source term is determined from the known angular distribution of radiation that exits the surface. Pomraning et al.  have coupled the diffusion approximation and the radiative transfer equation for boundary treatment for the steady state case.
With the advent of short laser pulses that are applied in various different disciplines, researchers have started developing techniques to solve the transient (time-dependent) PTE. Tan and Hsu  developed an integral equation formulation to treat the general transient PTE, and, Fleck used the Monte Carlo simulations . Kim, Ishimaru and Moscoso have applied Chebyshev collocation techniques for solving the PTE [12, 13]. Kim  has used Chebyshev expansion for the spatial discretization and Crank-Nicholson method for time marching. Handapangoda et al.  proposed a technique to transform the four-variable transient PTE to a set of coupled first order ordinary differential equations. However, in that work they solved this set of equations using the Runge-Kutta-Fehlberg method, which is a numerical technique. In this paper we propose an improvement to this technique, which gives a closed-form solution to the transformed PTE based on eigen decomposition. This approach improves the computational efficiency significantly as shown later.
This paper is organized as follows. In section 2 we present the eigen decomposition method to solve the reduced vectorized version of the PTE. Section 3 presents some numerical results. A comparison of the output values and the execution time of the algorithm is carried out against a those of the standard solution, which is obtained using the Laguerre Runge-Kutta-Fehlberg (LRKF) method . Section 4 carries some concluding remarks.
2. Eigen decomposition method
Suppose Eq. (1) is subject to the boundary condition I (z = 0, u0, ϕ0, t) = f (t)δ(u − u0)δ(ϕ − ϕ0), where δ(x) is the Dirac delta function  and f (t) is the temporal profile of the incident pulse. Using the transformationEq. (1) to a moving reference frame with the pulse, we get Eq. (2) in Eq. (1) eliminates the time derivative term in Eq. (1).
Gaussian quadrature  is used to approximate an integral of the form , where W (x) is the weight function and f (x) is any arbitrary function, by a summation such that .
Application of the discrete ordinates method to the azimuthal angle of Eq. (3) will result in a set of uncoupled equationsEq. (4) using the discrete ordinates method, which results in
In order to remove the dependency on τ, we can expand I (z, ui, ϕr, τ) and f (τ) using Laguerre polynomials with respect to τ, such that and where Bk and Fk are the coefficients corresponding to the Laguerre polynomial Lk (τ). Expanding the transformed one-dimensional PTE, Eq. (5) and the boundary condition using a Laguerre basis and taking moments we get,Eq. (6) corresponds to the ith discrete ordinate of u and rth discrete ordinate of ϕ. We can write this set of equations in matrix form as;
Rearranging, Eq. (8) can be written asEq. (7) can be decomposed into a set of values on the boundary corresponding to each discrete ordinate, which can be represented in a column matrix as .
An explicit solution to Eq. (9) can be found using eigen decomposition. The square matrix Y can be written as a product of three matrices composed of its eigen values and eigen vectors, as follows.Eq. (10) in Eq. (9) we get Eq. (11) is a diagonal matrix, each row of Eq. (11) can be solved independently resulting in Eq. (13), the solution to Eq. (9) can be found using the transformation in Eq. (12). That is, Bn = VX. Thus, the explicit solution to Eq. (9) is given by Eq. (9) can be obtained as illustrated in . For example, suppose λi = p + qi and λj = p − qi are two complex eigen values and vi = a + bi and vj = a − bi are their corresponding complex eigen vectors. It can be shown that we can then obtain two real valued solutions to Eq. (9), wi = epz(a cos qz – b sin qz) and wj = epz(b cos qz + a sin qz) . The general solution to Eq. (9) then becomes Equation (14) can be written in matrix notation as Eq. (15) (ie. ).20].
3. Results and discussion
Figure 1 shows a comparison of irradiance profiles obtained from the work proposed in this paper (Eigen method) and the standard solution (obtained using the LRKF method ), at the exit plane of a layer of a human skin specimen of 1 mm thickness, illuminated by a Gaussian pulse of I0 intensity (of arbitrary units). The Henyey-Greenstein phase function with an asymmetry factor of 0.7 together with σa = 0.5 mm−1, σs = 0.8 mm−1, v = 0.215 mm/ps and T = 1.5 was used for this simulation. The results produced by the two methods for the same number of discrete ordinates were exactly the same as shown by Fig. 1.
In the proposed algorithm numerical errors occur due to the truncation of the Laguerre series and due to the finite number of discrete ordinates chosen for θ and ϕ. However, as discussed in detail in  and , it is possible to obtain a very accurate Laguerre fit, in the observation window, with only 63 terms for a Gaussian pulse with an arbitrary width by using a suitable scaling factor. Hence, we can consider that the numerical errors are introduced only due to the discretization of the zenith and azimuthal angles. Figure 2 shows a comparison of execution times of the Eigen method proposed in this paper with the standard solution. The rate of increase of the execution time with the number of discrete ordinates is significantly low for the Eigen method compared to the standard solution.
From the simulation results it can be concluded that the proposed algorithm is accurate to the same degree compared to the standard solution which is based on a numerical solution to the reduced vectorized version of the PTE. However, the proposed algorithm is significantly faster in execution time compared to the standard solution method. For example, increasing the number of directions from 25 to 196 increased the execution time of the proposed algorithm by only about 11 seconds (6.4%) where as that of the standard method was increased by about 166 seconds (97.0%).
References and links
1. A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202, 323–345 (2005). [CrossRef]
2. N. Y. Gnedin, “Multi-dimensional cosmological radiative transfer with a variable Eddington tensor formalism,” New Astron. 6, 437–455 (2001). [CrossRef]
3. D. M. O’Brien, “Accelerated quasi Monte Carlo integration of the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. 48, 41–59 (1992). [CrossRef]
4. K. Stamnes, S. C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27, 2502–2509 (1988). [CrossRef] [PubMed]
5. J. L. Deuz, M. Herman, and R. Santer, “Fourier series expansion of the transfer equation in the atmosphere-ocean system,” J. Quant. Spectrosc. Radiat. Transf. 41, 483–494 (1989). [CrossRef]
6. C. E. Siewert and J. R. Thomas, “Radiative transfer calculations in spheres and cylinders,” J. Quant. Spectrosc. Radiat. Transf. 34, 59–64 (1985). [CrossRef]
7. C. E. Siewert, “A radiative-transfer inverse-source problem for a sphere,” J. Quant. Spectrosc. Radiat. Transf. 52, 157–160 (1994). [CrossRef]
8. E. W. Larsen, “The inverse source problem in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 15, 1–5 (1975). [CrossRef]
9. G. C. Pomraning and G. M. Foglesong, “Transport-diffusion interfaces in radiative transfer,” J. Comput. Phys. 32, 420–436 (1979). [CrossRef]
10. Z. M. Tan and P. F. Hsu, “An integral formulation of transient radiative transfer,” ASME J. Heat Transfer 123, 466–475 (2001). [CrossRef]
11. J. A. Fleck, “The calculation of nonlinear radiation transport by a Monte Carlo method,” Methods Comput. Phys. 1, 43–65 (1963).
12. A. D. Kim and A. Ishimaru, “A Chebyshev spectral method for radiative transfer equations applied to electromagnetic wave propagation and scattering in discrete random media,” J. Comput. Phys. 152, 264–280 (1999). [CrossRef]
13. A. D. Kim and M. Moscoso, “Chebyshev spectral methods for radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). [CrossRef]
14. C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, “Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue,” IEEE J. Sel. Top. Quantum Electron. 14, 105–112 (2008). [CrossRef]
15. A. D. Kim and M. Moscoso, “Chebyshev spectral methods fro radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). [CrossRef]
16. A. B. Carlson, Communication Systems: An Introduction to Signals and Noise in Electrical Communication (McGraw-Hill, 1986).
17. S. Chandrasekhar, Radiative Transfer (Dover Publications, 1960).
18. K. Stamnes and R. A. Swanson, “A new look at the discrete ordinate method for radiative transfer calculations in anisotropically scattering atmospheres,” J. Atmos. Sci. 38, 387–399 (1981). [CrossRef]
19. W. H. Press, S. A. Teukolsk, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C++, 2nd ed. (Cambridge University Press, 2003).
20. C. H. Edwards and D. E. Penney, Differential Equations: Computing and Modeling, 3rd ed. (Prentice Hall, 2004).
21. C. C. Handapangoda and M. Premaratne, “Implicitly causality enforced solution of multidimensional transient photon transport equation,” Opt. Express 17, 23423–23442 (2009). [CrossRef]