This paper presents a third-order diffusion equations-based optical image reconstruction algorithm. The algorithm has been implemented using finite element discretizations coupled with a hybrid regularization that combines both Marquardt and Tikhonov schemes. Numerical examples are used to compare between the third- and first-order reconstructions. The results show that the third-order reconstruction codes are more stable than the first-order codes, and are capable of reconstructing void-like regions. From the examples given, it has also been shown that the first-order codes fail to both qualitatively and quantitatively reconstruct the void-like regions.
© Optical Society of America
Current optical image reconstruction algorithms [1–11] are almost all based on the first-order diffusion equation, which is only applicable in cases where the scattering is assumed to dominate the absorption and the optode spacing is much larger than the inverse of the reduced scattering coefficient. While these conditions are satisfied for almost all tissues in the near-infrared region, there are concerns in the use of the first-order diffusion equation in optical image reconstruction. The major concern is encountered in imaging multi-layered brain tissue where a clear, non-scattering layer of cerebrospinal fluid (CSF) lies between the inner skull table and the brain surface. The first-order diffusion equation fails to describe light propagation in this region [12–14], which means that alternative approaches for light modeling must be used. Another concern involved in optical image reconstruction using the first-order diffusion equation is highly absorbing regions such as hematoma.
In this paper we attempt to develop a third-order diffusion equations-based reconstruction algorithm that may resolve the above mentioned concerns. Since the third-order diffusion equations can be derived from the Boltzmann transport equation without the assumptions that are applied to the first-order diffusion equation, they are applicable to the CSF layer in brain tissue. Further, because the third-order diffusion equations are hyperbolic type differential equations, they can provide more stable inverse solutions than the parabolic first-order diffusion equation . Numerical examples have been used to confirm these conclusions using the implemented third-order diffusion equations-based reconstruction algorithm in which finite element discretizations coupled with a synthesized Marquardt and Tikhonov regularization scheme have been used.
2. Reconstruction algorithm
In this paper we are interested only in optical image reconstruction in the steady-state, continuous-wave domain. The algorithm in the frequency- or time-domain can be developed in a similar manner. The third-order diffusion equations derived from the Boltzmann photon transport equation can be stated as[16–18]:
where and Φ(1), Φ(2), Φ(3), and Φ(4) are the first four components in the spherical harmonic expansion of the photon radiance where the first component, Φ(1) , is the average diffused photon density. μa is the absorption coefficient. μ′t = μa + (1 - g)μs, where μs is the scattering coefficient and g is the average cosine of the scattering angle. D = 1 / 3μ′t is the diffusion coefficient. x̂ and ŷ are the unit vectors along x and y axes, respectively. S is the light source term.
We must choose appropriate boundary conditions (BCs) in order to solve Eqs. (1)–(4). In this study we apply Type III BCs to the first component, Φ(1) , i.e., -Dn̂·∇Φ(1) = αΦ(1), where n̂ is the unit normal vector for the boundary surface and α is a coefficient that is related to the internal reflection at the boundary, and employ Type I BCs to the remaining components, i.e., Φ(2) = Φ(3) = and Φ(4) = 0.
Making use of finite element discretizations, the discretized forms of Eqs. (1)–(4) can be written as
where 〈 〉 indicates integration over the problem domain and Φ(1)–(4) , D, and μa have been expanded as the sum of coefficients multiplied by a set of locally spatially-varying Lagrangian basis functions ϕ ϕp, and ϕq. ∮ expresses integration over the boundary surface. N is the node number of a finite element mesh. The expansions used to represent D and μa are P and Q terms long where P ≠ Q ≠ N in general; however, in the study reported here P=Q=N.
In optical image reconstruction, the goal is to form D and μa images from presumably uniform initial estimates of the spatial D and μa distributions. Thus, we need a way of updating D and μa from their starting values. Following the inverse procedures outlined in [8,9], the following matrix equation for updating D and μa is obtained:
and ∆χ = (∆D1, ∆D2,…∆DN,∆μa,1, ∆μa,2,…∆μa, N)T is the update vector for the optical property profiles. Φo = ( , ,…)T and Φc = ( , ,…)T , where and are observed and calculated average diffused photon density for i=1, 2, … M boundary locations. Note that only the first component or the average diffused photon density, Φ(1) , is used in Eq. (9) since it is the dominant component , and other components, Φ(2)–(4) , are set to zeros at the boundary. In Eq. (9), the decomposition of the ill-conditioned matrix ℑ T ℑ is stabilized by a synthesized Marquardt and Tikhonov regularization scheme [8,9].
3. Numerical examples of reconstruction
In this section we use numerical examples to test the reconstruction algorithm described in Section 2. The test geometry, shown in Fig. 1 (a), consists of a circular background region (radius=21.5 mm) with an embedded circular target (radius=6.25 mm) offsetting 5 mm. The examples include two test cases with different optical properties assigned in the embedded target and background regions. For the first case, the optical properties for the target are μ′s =2.0 mm-1, μa=0.012 mm-1; the optical properties for the background are μ′s =1.0 mm-1, μa=0.006 mm-1. For the second case, the optical properties for the target are μ′s =0.01 mm-1, μa=0.005 mm-1; the optical properties for the background are μ′s =1.0 mm-1, μa=0.01 mm-1. The first case is used to just demonstrate the implementation of our third-order reconstruction codes. The purpose of the second case is to test if the third-order codes can reconstruct a void-like region and if it can provide more stable reconstructions than the first-order codes when noisy data is used. The optical properties assigned to the void-like region is similar to that in the CSF layer in brain tissue . Multiple excitation and measurement positions were used to produce the boundary information used in the reconstructions. Specifically, we used 16 excitation positions (equally spaced around the circular circumference) and 16 measurement locations (also equally spaced around the circular circumference, but with a shift relative to the excitation positions) for detection of diffusive light. The radial location of each source was positioned inside of the physical boundary by a distance, d=3D for the point source excitation used in the computational algorithm. The finite element mesh used in this study consisted of 241 nodes and 416 triangle elements. The final images reported are the result of iteration until the initial sum of squared errors between measured and computed photon density values at the measurement site locations is reduced 5 orders of magnitude. Reaching this level of reduction in the initial sum of squared errors typically required 20 iterations at a cost of 2 minutes per iteration for the finite element mesh used herein in a Sun Ultra 30 workstation.
In the examples, the “measured” data were generated using a forward higher-order diffusion model with the exact D and μa in place. Fig. 1 (b, c) shows the D and μa images for the first case reconstructed under conditions of no noise. As can be seen, the images are clearly recovered.
For the second case, 2% noise has been added to the “measured” data. Fig. 2 (a, b) presents the successfully reconstructed D and μa images for the second case. In order to provide a comparison, D and μa images reconstructed using our first-order codes are displayed in Fig. 2 (c, d). A number of observations can be made from Fig. 2. The almost non-scattering, void-like target can be qualitatively recovered for both D and μa images using the third-order codes [Fig. 2 (a, b)], whereas it cannot be correctly recovered using the first-order codes [Fig. 2 (c, d)]. Interestingly, the target location for both D and μa images recovered is incorrectly “swapped” to the left when the firs-order codes were used. From Fig. 2 (a, b), one can see that the third-order codes can produce correct reconstructions of the target location and shape. While the reconstructed target size for D image is correct [Fig. 2 (a)], the recovered target size for μa image is larger than the exact target size. When the first-order codes were used, the recovered target size, location and shape for both D and μa images are totally incorrect. From Fig. 2, it can be seen that both D and μa images can be quantitatively reconstructed using the third-order codes [Fig. 2 (a, b)], whereas the recovered values of both D and μa in the target region are all “swapped” with respect to the exact values when the first-order codes were used [Fig. 2 (c, d)]. However, it is interesting to note that the first-order codes produce better background region reconstruction than the third-order codes. Given the facts that the third-order codes can reconstruct the void-like regions from noisy data and the first-order codes cannot do so, one can also see that the third-order codes are more stable than the first-order codes.
We have demonstrated optical image reconstructions using the higher-order diffusion equations-based reconstruction algorithm. The algorithm presented here coupled with the imaging enhancement schemes we developed [10,11] will be able to improve the overall quality of optical image reconstruction. It is anticipated that this algorithm may find its applications in imaging brain tissues.
This work was supported in part by the National Institutes of Health (R55 CA78334) and the Greenville Hospital System/Clemson University Biomedical Cooperative.
References and links
1. R. L. Barbour, H. L. Graber, Y. Wang, J. H. Chang, and R. Aronson, “A perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data,” in Medical Optical Tomography, G. Miller ed., SPIE Institute for Advanced Optical Technologies (SPIE Optical Engineering Press Vol. IS11, 1993), 87–120.
3. M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusion-photon tomography,” Opt. Lett. 20, 426–428 (1995). [CrossRef]
4. X. D. Li, T Durduran, A. G. Yodh, B. Chance, and D. N. Pattanayak, “Diffraction Tomography for biomedical imaging with diffuse-photon density waves,” Opt. Lett. 22, 573–575 (1997). [CrossRef] [PubMed]
6. S. A. Walker, S. Fantini, and E. Gratton, “Image reconstruction by backprojection from frequency domain optical measurements in highly scattering media,” Appl. Opt. 36, 170–179 (1997). [CrossRef] [PubMed]
7. S. B. Colak, D. G. Papaioannou, G. W. ’t Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens, J. B. M. Melissen, and N. A. A. J. van Asten, “Tomographic image reconstruction from optical projections in light-diffusing media,” Appl. Opt. 36, 180–213 (1997). [CrossRef] [PubMed]
9. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical image reconstruction using frequency-domain data: simulations and experiments,” J. Opt. Soc. Am. A 13, 253–266 (1996). [CrossRef]
11. H. Jiang, K. D. Paulsen, U. L. Osterberg, and M. S. Patterson, “Frequency-domain near-infrared photo diffusion imaging: initial evaluation in multi-target tissue-like phantoms,” Med. Phys. 25, 183–193 (1998). [CrossRef] [PubMed]
12. D. T. Delpy and M. Cope, “Quantification in tissue near-infrared spectroscopy,” Phil. Trans. R. Soc. Lond. B 352, 649–659 (1997). [CrossRef]
13. A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43, 1285–1302 (1998). [CrossRef] [PubMed]
14. M. Firbank, S. Arridge, M. Schweiger, and D. Delpy, “An investigation of light transport through scattering bodies with non- scattering regions,” Phys. Med. and Biol. 41767–783 (1996). [CrossRef]
15. In Mathematics and Physics of Emerging Biomedical Imaging (National Academy Press, Washington, D.C., 1996).
16. H. Jiang and K. D. Paulsen, “A finite element based higher-order diffusion approximation of light propagation in tissues,” Proc. SPIE 2389, 608–614 (1995). [CrossRef]
17. D. A. Boas et al, ”Photon Migration within the P3 approximation,” in Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media, Proc. SPIE 2389, pp. 240–247 (1995). [CrossRef]
18. D. A. Boas, “A fundamental limitation of linearized algorithms for diffuse optical tomography,” Optics Express 1, 404–413 (1997), http://epubs.osa.org/oearchive/source/2831.htm [CrossRef] [PubMed]