We study the modeling and simulation of steady-state measurements of light scattered by a turbid medium taken at the boundary. In particular, we implement the recently introduced corrected diffusion approximation in two spatial dimensions to model these boundary measurements. This implementation uses expansions in plane wave solutions to compute boundary conditions and the additive boundary layer correction, and a finite element method to solve the diffusion equation. We show that this corrected diffusion approximation models boundary measurements substantially better than the standard diffusion approximation in comparison to numerical solutions of the radiative transport equation.
© 2012 Optical Society of America
Non-invasive boundary measurements of light scattered by tissues are important for biomedical applications [1, 2]. By extracting information from these measurements regarding the optical properties of tissues (e.g. absorption and scattering), one may gain valuable insight into tissue health. For example, in diffuse optical tomography (DOT) and fluorescence diffuse optical tomography (fDOT), one seeks to reconstruct the optical properties of tissues from measurements of light at the boundary of the domain. The applications include, for example, detection and classification of breast cancer, monitoring of infant brain tissue oxygenation level and functional brain activation studies, for reviews see e.g. [3–5]. In quantitative photoacoustic tomography (QPAT) one seeks to estimate concentration of chromophores, such as hemoglobin and melanin, inside tissues by combining the optical contrast and ultrasound propagation [6–8].
Image reconstruction problems in DOT, fDOT and QPAT are non-linear ill-posed inverse problems. There are no direct methods for the solution of these problems, and thus they are typically stated as minimization problems such as regularized output least squares. The iterative solution of this problem requires repetitive solutions of the forward model. Therefore, it is essential to have a computationally feasible forward model that describes light propagation in tissues accurately.
The theory of radiative transport governs light propagation in tissues [1, 2]. This theory takes into account absorption and scattering due to inhomogeneities in the medium. The major challenge in using radiative transport theory to study light propagation in tissues is that it is mathematically complicated due mostly to the large number of variables in the radiative transport equation (RTE). The large dimensionality of the RTE makes even computational methods challenging.
Tissues typically scatter light strongly and absorb light weakly. For that reason, one often replaces the RTE by the diffusion approximation (DA) [1–4]. In the DA, one assumes that the light becomes nearly isotropic due to strong multiple scattering. The DA is much simpler to solve than the RTE. However, applying the DA to model boundary measurements is problematic. It is well known that the DA is not valid near sources or boundaries. This is because the assumption that light is nearly isotropic is too restrictive to take into account sources and boundary conditions. Nonetheless, the DA has been used to model boundary measurements with some success despite these limitations. Regardless, there still exists a need for more accurate models of boundary measurements. Consequently, the prescription of “correct” boundary conditions [9–13] and source terms [14–16] for the DA has been a long-standing issue.
There have been some works that have taken into account sources and boundaries correctly by combining the solutions of the RTE and the DA to form a hybrid method. Wang and Jacques  used Monte Carlo simulations for the RTE in combination with the DA. Tarvainen et al  developed a coupled method combining solutions of the RTE and DA both within a finite element framework for both space and angle. This coupled method can take into account correctly boundaries, sources as well as low-scattering regions in the interior of the domain. Recently, Gao and Zhao  developed a sophisticated numerical method to solve the RTE. This method employs multigrid methods in both space and angle where the coarsest grid level is consistent with the DA when scattering is strong and absorption is weak.
Recently, Kim  presented an asymptotic analysis of the RTE leading to the so-called corrected diffusion approximation (cDA). The additive correction to the DA is given by a boundary layer solution. This boundary layer solution corrects for the error made by the DA near the boundary. It vanishes rapidly away from the boundary (on the order of one scattering mean free path) so that the standard DA approximates the solution deep within the interior of the domain. The asymptotic analysis used to derive this boundary layer solution was established in the early 1970’s in the neutron transport community [21, 22]. In fact, Pomraning and Ganapol  used this asymptotic analysis to study boundary conditions for the DA in detail.
In this paper, we use the cDA to model boundary measurements. Using a finite element method (FEM) to solve the DA, we are able to consider general spatial domains. By computing the boundary layer solution only for boundary points where we are modeling measurements, we show that the cDA provides a superior approximation to the solution of the RTE requiring only a small amount of more work than solving the DA itself. In particular, we consider a spatial domain Ω ⊂ 2 with boundary ∂Ω. Although the simulations presented here are limited to two-dimensional (2D) case, the theory derived in  and reviewed here is represented in dimensional independent form allowing for the realization of the method both in two and three dimensions (3D). In 3D, the cDA consist of solving 3D diffusion equation and one-dimensional radiative transport equation for each of the measurements locations. Thus, the extension of method to 3D is straightforward.
The remainder of the paper is as follows. In Section 2, we give an overview of the asymptotic analysis of the RTE leading to the cDA. In Section 3 we give details of the calculations needed to compute the boundary condition coefficients for the diffusion equation and the boundary layer solution. In addition, we give details of the numerical method used to solve the diffusion equation. In Section 4 we show results from our computational simulations. Section 5 is the conclusions.
2. Asymptotic analysis of the radiative transport equation
The steady-state radiative transport equation
To solve Eq. (1) in Ω × Sn−1, we must supplement boundary conditions of the formAppendix B, Eq. (B.4)). Boundary condition Eq. (3) prescribes the radiance over only the directions pointing into the domain.
Boundary measurements are given by the exitance Γ(rb) defined asAppendix B, Eq. (B.15)). Furthermore, other quantity of interest is the fluence rate ϒ(r), defined as an integral of the radiance over angular directions
We consider the case in which scattering is strong and absorption is weak. To make this assumption explicit, we introduce a small, dimensionless parameter 0 < ε ≪ 1 according toEq. (6) into Eq. (1), we obtain Eq. (7) in the limit as ε → 0+ of the form 20]. In what follows, we summarize the results from that work.
Seeking the interior solution Φ in the formEq. (3). For that reason, we add a boundary layer solution Ψ to correct the interior solution near the boundary. This boundary layer solution decays rapidly away from the boundary on a length scale that is O(ε). Thus, ϕ ∼ Φ0 − εnκŝ · ∇Φ0 deep in the interior of Ω far away from the boundary ∂Ω.
To compute Ψ near a particular boundary point rb ∈ ∂Ω, consider a coordinate system (ρ,z) where ρ is a vector parallel to the tangent plane at rb and z is the coordinate along −n̂. Let z = εζ, μ = ŝ · (−n̂) and ŝ⊥ = ŝ + μn̂. Then, each of the terms in the boundary layer solution Ψ(ρ,ζ, ŝ) = Ψ0(ρ,ζ, ŝ) + εΨ1(ρ,ζ,ŝ) + O(ε2) satisfies boundary value problems for the one-dimensional radiative transport equations of the form
To ensure asymptotic matching condition Eq. (15) is satisfied, one must set
For the special case in which the boundary source is axisymmetric about n̂ so that ϕ0 = ϕ0(rb,μ), the right-hand side of Eq. (14a) vanishes identically and Eq. (16) reduces to the familiar Robin boundary conditionEquation (10) together with the boundary condition Eq. (17) is known as the diffusion approximation to the RTE. The standard diffusion approximation, which is derived from the spherical harmonics expansion of the RTE, is shown in Appendix A.
Upon solution of boundary value problem Eq. (21), the corrected diffusion approximation evaluated at rb ∈ ∂Ω is then given byEquation (22) gives the asymptotic solution up to O(ε2).
3. Numerical implementations
In this section, we give a method used to compute the boundary condition coefficients given by Eqs. (18) – (20) and the solution of the boundary value problem Eq. (21) given the diffusion approximation. Furthermore, we describe briefly a finite element method to solve Eq. (10) subject to boundary condition Eq. (17).
3.1. Computing boundary condition coefficients and boundary layer solutions
To compute the coefficients of the boundary conditions for the diffusion approximation and the boundary layer solution for Ω ⊂ 2, we need to study the canonical half space problem:24, 25]
Since boundary condition Eq. (23b) prescribes an even function of θ at ζ = 0 and because the scattering phase function given in Eq. (24) is rotationally invariant, the solution is an even function of θ:Eq. (23) as
3.1.1. Plane wave solutions
We use plane wave solutions to solve boundary value problem Eq. (27). Plane wave solutions are special solutions of Eq. (27a) of the form Ψ = eλζV(μ). Substituting this ansatz into Eq. (27a), we obtain the eigenvalue problem26, 27].
To calculate plane wave solutions numerically, we use the discrete ordinate method. In particular, we use the Gauss-Chebyshev quadrature ruleEq. (27a) with the Gauss-Chebyshev quadrature rule and evaluating that result at μi, we obtain Equation (32) is a discrete eigenvalue problem suitable for numerical computations. Notice that in Eq. (32) we have added a small regularization parameter 0 < δ ≪ 1 which is equivalent to adding a small amount of absorption. This ensures that the eigenvalues we calculate numerically are real and distinct. The numerical error incurred by introducing δ is exponentially small. For our numerical calculations, we typically have chosen δ = 10−8.
Suppose we solve numerically Eq. (32). We will obtain N eigenvalues λn and eigenvectors Vn(μi). For each pair [λn,Vn(μi)] satisfying Eq. (32), the pair [−λn,Vn(−μi)] satisfies Eq. (32) also. As a result, we order and index the eigenvalues according to
3.1.2. Boundary condition coefficients
In Eqs. (18) – (20), the coefficients a, b and f are defined in terms of a projection operator 𝒫. This operator is given as an expansion in plane wave solutions derived in . Here, we state the result. Let pi for i = N/2 + 1,···,N be defined as
3.1.3. Boundary layer solution
Now that we have computed the boundary condition coefficients in boundary condition Eq. (17), let us suppose that we have solved Eq. (10) subject to boundary condition Eq. (17). We show how we solve this boundary value problem numerically in Section 3.2. Then, the boundary layer solution Ψ satisfying boundary value problem Eq. (21) can be computed as an expansion in plane wave solutions. This expansion is derived in . Here, we give the result for Ψ(0,μi) for −1 ≤ μ ≤ 1 which is what we need to correct the DA at the boundary. Let Hij be defined as
3.2. Computing the diffusion approximation
3.3. Summary of the algorithm
To summarize, the procedure for a numerical solution of the cDA is given as
4. Numerical results
The performance of the cDA was tested with 2D simulations. Simulation domain Ω was a circle with radius of 20 mm centered at the origin. The source ϕ0(r,ŝ) was located at (x,y) = (−20,0) mm with cosine shape giving the largest value in inward direction and value zero in direction of tangent to the boundary
Three types of test cases were considered: a homogeneous medium with matched refractive indices inside and outside the domain for three different values of scattering coefficient μs, the same cases with mismatched refractive indices, and a heterogeneous medium.
The results of the cDA were compared with the results of the standard DA (see Appendix A) and the RTE. The cDA was solved as explained above in Section 3.3. The standard DA was solved with the FEM similarly as in . The FE-approximation of the RTE was implemented similarly as in [32,33] in the case of matched refractive indices. The implementation is explained in more detail in Appendix B.
The FE-mesh for the spatial discretization of the domain contained approximately 4687 nodal points and 9196 triangular elements for the homogeneous test cases and 6114 nodal points and 11 988 elements for the heterogeneous test case. For the angular discretization of the RTE 64 equally spaced angular directions were used. Parameter ε was chosen as the ratio of the total mean free path l = (μs + μa)−1 to the size of the domain L 
4.1. Matched refractive indices
For the first case, we consider a medium with matched refractive indices at the boundary. The refractive indices inside and outside the medium were nin = 1 and nout = 1, respectively. The scattering coefficient was given three value μs = 50, 5 and 0.5 mm−1. The absorption coefficient and the anisotropy factor were constants, μa = 0.01 mm−1 and g = 0.8, respectively.
Radiances at the boundary point (x,y) = (0,−20) for μs = 5 mm−1 computed using the cDA, the DA and the RTE are shown in left image of Fig. 1. To compare the performance of the cDA against the DA, two different quantities were investigated. First, the fluence rate inside the domain was computed using Eq. (5), and secondly the exitance at the boundary was computed using Eq. (4). The exitances computed using the cDA, the DA and the RTE are shown in Fig. 2. In addition, the relative errors of fluence rates computed using the cDA (top row) and the DA (bottom row) against the RTE for different values of μs are shown in Fig. 3. Furthermore, the relative errors of exitances are shown in Fig. 4. To give a quantitative estimate of the errors, the means of the relative errors of fluence rates ΔΦ0 and the exitances ΔΓ were computed. The results are given in Table 1. We also recorded the computation times of the different models. These are given in Table 2.
As it can be seen from Fig. 1, the radiance computed using the cDA agree relatively well with the RTE and satisfy the zero boundary condition in inward direction. Note that the RTE may give small negative or positive values in inward direction due to numerical reasons. The DA gives negative radiance in inward direction which are unphysical. The relative error of fluence rate is approximately 2 % for the cDA and between 4–6 % for the standard DA when scattering coefficient is μs = 50 and μs = 5 mm−1 as it can be seen from Fig. 3 and Table 1. For μs = 0.5 mm−1 both the cDA and the standard DA give large errors since the assumptions of the DA are not valid anymore. Figure 4 shows that as scattering becomes large compared to absorption, the relative error of exitance decreases for the cDA just as the asymptotic theory predicts.
The computation times in Table 2 show that solving the cDA is feasible when compared with the standard DA. In addition, solving both the cDA and the DA are much faster than solving the RTE. Thus, the cDA models light propagation more accurately than the standard DA and requires only a small amount of more work.
4.2. Mismatched refractive indices
For the second case, we consider a medium with mismatched refractive indices at the boundary. The refractive indices inside and outside the medium were nin = 1.33 and nout = 1, respectively. Other optical parameters were the same as before. Radiances at the boundary point (x,y) = (0,−20) for μs = 5 mm−1 computed using the cDA, the DA and the RTE are shown in right image of Fig. 1. Again, the quantities of interest were the fluence rate inside the domain and the exitance at the boundary. The exitances computed using different models are shown in Fig. 5. In addition, the relative errors of fluence rates and exitances are shown in Figs. 6 and 7, respectively. As before, the means of the relative errors of fluence rates and exitances were computed and the results are given in Table 3.
As it can be seen Fig. 1, the approximation to the radiance given by the cDA agrees better with the RTE than the approximation given by the DA. Furthermore, the results in Fig. 6 and Table 3 show the relative error of fluence rate is between 2–3 % for the cDA when μs = 50 mm−1 and μs = 5 mm−1. For the DA the relative error is larger giving the largest errors at the boundary. For the case μs = 0.5 mm−1 neither of the approximations are valid. Figure 7 shows that the relative error of exitance decreases for the cDA when the ratio of absorption and scattering decreases due to the asymptotic theory behind the model. In contrast, the relative error of the standard DA may not have a global error bound over the whole domain. In that case, the relative error might be large close to the boundary while giving satisfactory results inside the domain. In addition, decrease in ratio of absorption and scattering does not ensure decrease in relative error for the DA. Thus, the cDA models light propagation more accurately than the standard DA when compared with the RTE.
4.3. Heterogeneous medium
For the third case, we consider a medium with heterogeneous optical properties. Simulated optical parameters of the medium are shown in Fig. 8. The scattering and absorption coefficients of the background were (μs,μa) = (10,0.01) mm−1. Furthermore, optical parameters of scattering and absorbing inclusions were (μs,μa) = (20,0.01) mm−1 and (μs,μa) = (10,0.02) mm−1, respectively. The anisotropy factor in the whole domain was g = 0.8. The refractive indices inside and outside the medium were nin = 1.33 and nout = 1, respectively. Asymptotic parameter ε = 0.005 was computed from Eq. (50) using the optical parameters of the background. The exitances are shown in Fig. 9. The relative errors of fluence rates computed using the cDA and the DA are shown in Fig. 10. Furthermore, the relative errors of exitances are shown in Fig. 11.
The mean of the relative error of fluence rate is 1.69 % for the cDA over whole domain. In contrast, the mean of the relative error of fluence rate is 5.07 % for the DA. In addition, the DA gives larger errors close to the boundary than inside the domain as it can be seen from Fig. 10. Furthermore, the relative error of exitance is smaller for the cDA (mean 5.38 %) than that of for the DA (mean 18.73 %). Thus, the cDA gives more accurate results over whole domain than the standard DA due to a global error bound given by the asymptotic theory.
In this work, recently introduced corrected diffusion approximation was numerically implemented. In the cDA, an additive correction term is computed for the DA at the boundary based on asymptotic analysis of the RTE. The procedure for computing the cDA requires only small modifications to the existing solvers for the DA. In particular, one only needs to modify the spatial coefficients of the DA and the Robin boundary condition. In addition, an additive correction term, which satisfies a one-dimensional radiative transport equation, is readily solved using the plane wave solutions.
The performance of the cDA was tested with 2D simulations. The results were compared with the results of the DA and the RTE. The results show that the cDA models boundary measurements of scattered light more accurately than the standard DA with only a small increase in computation time.
Appendix A. The standard diffusion approximation
The most typical approach to derive the standard DA from the RTE is expand the radiance, the source term, and the phase function into series using the spherical harmonics and truncate the series [1,3]. The standard DA can be regarded as a special case of the first order approximation. In the standard DA framework, the approximation that is used for the radiance is of the form [1, 3, 35]36]. The parameters A takes into account a mismatch in refractive indices inside and outside the domain and it can be derived from the Fresnel’s law [36, 37]
Appendix B. FE-approximation of the RTE with reflection boundary condition
Let us consider the steady-state RTE Eq. (1). Let ŝi and ŝr be the directions of incident and reflected light, respectively, at the boundary with a mismatch in refractive indices. Furthermore, let ŝt be the direction refracted light. We write the boundary condition for the radiance pointing inward into the domain (i.e. for the directions ŝr · n̂ < 0) in terms of the incident radiance pointing outward (ŝi · n̂ > 0). Mapping from direction ŝi to direction ŝr is given by H : ŝi → ŝrEq. (3) can be written as follows Eq. (B.9) are
After the RTE is solved using Eq. (B.9), the direction ŝt of the refracted radiance at the boundary exiting the domain can be computed from Snell’s law which is in vector form K : (ŝi, n̂, nin, nout) → ŝt
Thus, the final form of the radiance after applying the Snell’s law is of the form
The authors would like to thank professor Simon Arridge from useful discussions. The work of O. Lehtikangas and T. Tarvainen has been supported by the Academy of Finland (projects 136220, 140984, and 213476, 250215 Finnish Centre of Excellence in Inverse Problems Research), Finnish Doctoral Programme in Computational Sciences (FICS), Magnus Ehrnrooth foundation and by the strategic funding of the University of Eastern Finland. A. D. Kim acknowledges support from the National Science Foundation (grant DMS-0806039).
References and links
1. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE, New York, 1997).
2. L.V. Wang and H. Wu, Biomedical Optics: Principles and Imaging (Wiley, Hoboken, NJ, 2007).
3. S.R. Arridge, “Optical tomography in medical imaging,” Inv. Prob. 15, R41–R93 (1999). [CrossRef]
4. S.R. Arridge and J.C. Schotland, “Optical tomography: forward and inverse problems,” Inv. Prob. 25123010 (2009). [CrossRef]
6. B.T. Cox, S.R. Arridge, K.P. Köstli, and P.C. Beard, “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt. 45, 1866–1875 (2006). [CrossRef] [PubMed]
7. G. Bal and K. Ren, “Multi-source quantitative photoacoustic tomography in a diffusive regime,” Inv. Prob. 27, 075003 (2011). [CrossRef]
8. A.Q. Bauer, R.E. Nothdurft, T.N. Erpelding, L.V. Wang, and J.P. Culver, “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt. 16, 096016 (2011). [CrossRef] [PubMed]
9. R. Aronson, “Extrapolation distance for diffusion of light,” Proc. SPIE 1888, 297–304 (1993). [CrossRef]
10. R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
11. R. Aronson, “Boundary conditions for the diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]
12. R. Aronson and N. Corngold, “Photon diffusion coefficient in an absorbing medium,” J. Opt. Soc. Am. A 16, 1066–1071 (1999). [CrossRef]
13. J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A 16, 1947–1957 (1999). [CrossRef]
15. X. Intes, B. Le Jeune, F. Pellen, Y. Guern, J. Cariou, and J. Lotrian, “Localization of the virtual point source used in the diffusion approximation to model a collimated beam source,” Waves Random Media 9, 489–499 (1999). [CrossRef]
16. T. Spott and L. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. 39, 6453–6465 (2000). [CrossRef]
17. L.-H. Wang and S. L. Jacques, “Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media,” J. Opt. Soc. Am. A 10, 1746–1752 (1993). [CrossRef]
18. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. P. Kaipio, “Finite element model for the coupled radiative transfer equation and diffusion approximation,” Int. J. Num. Math. Eng. 65, 383–405 (2006). [CrossRef]
19. H. Gao and H. Zhao, “A fast forward solver of radiative transfer equation,” Transp. Theory Stat. Phys 38, 149–192 (2009). [CrossRef]
20. A.D. Kim, “Correcting the diffusion approximation at the boundary,” J. Opt. Soc. Am. A 28, 1007–1015 (2011). [CrossRef]
21. E.W. Larsen and J. B. Keller, “Asymptotic solution of neutron transport problems for small mean free paths,” J. Math. Phys. 15, 75–81 (1974). [CrossRef]
22. G.J. Habetler and B. J. Matkowsky, “Uniform asymptotic expansions in transport theory with small mean free paths, and the diffusion approximation,” J. Math. Phys. 16, 846–854 (1975). [CrossRef]
23. G.C. Pomraning and B. D. Ganapol, “Asymptotically consistent reflection boundary conditions for diffusion theory,” Ann. Nucl. Energy 22, 787–817 (1995). [CrossRef]
24. L.G. Henyey and J.L. Greenstein, “Diffuse radiation in the Galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]
25. J. Heino, S.R. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (2003). [CrossRef]
26. A.D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20, 92–98 (2003). [CrossRef]
27. A.D. Kim, “Transport theory for light propagation in biological tissue,” J. Opt. Soc. Am. A 21, 820–827 (2004). [CrossRef]
31. V. Kolehmainen, S.R. Arridge, W.R.B Lionheart, M. Vauhkonen, and J.P. Kaipio, “Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data,” Inv. Prob. 15, 1375–1391 (1999). [CrossRef]
33. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, S.R. Arridge, and J.P. Kaipio, “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol , 50, 4913–4930 (2005). [CrossRef] [PubMed]
34. A.D. Kim and M. Moscoso, “Diffusion of polarized light,” Multiscale Model. Simul. 9, 1624–1645 (2011). [CrossRef]
35. T. Tarvainen, V. Kolehmainen, S.R. Arridge, and J.P. Kaipio, “Image reconstruction in diffuse optical tomography using the coupled radiative transport-diffusion model,” J. Quant. Spect. Rad. Trans. 1122600–2608 (2011). [CrossRef]
36. M. Schweiger, S.R. Arridge, M. Hiraoka, and D.T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]
38. G. Kanschat, “A robust finite element discretization for radiative transfer problems with scattering,” East West J. Num. Math. 6, 265–272 (1998).
39. S. Richling, E. Meinköhn, N. Kryzhevoi, and G. Kanschat, “Radiative transfer with finite elements,” Astron. Astrophys. 380, 776–788 (2001). [CrossRef]