Optical image reconstruction based on the third-order diffusion equations

Huabei Jiang

doi:10.1364/OE.4.000241

1. Introduction

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 [15]. 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]:

\nabla \cdot D (r) \nabla Φ^{(1)} (r) - μ_{a} (r) Φ^{(1)} (r) - \nabla \cdot D (r) \nabla Φ^{(2)} (r) + 6 \nabla \cdot D (r) {\nabla_{1} Φ}^{(3)} (r) + 6 \nabla \cdot D (r) {\nabla_{2} Φ}^{(4)} (r) = - S (r)

- \nabla \cdot D (r) \nabla Φ^{(1)} (r) + \frac{25}{7} \nabla \cdot D (r) {\nabla Φ}^{(2)} (r) - 5 {μ'}_{t} (r) Φ^{(2)} (r) - \frac{60}{7} \nabla \cdot D (r) {\nabla_{1} Φ}^{(3)} (r) - \frac{60}{7} \nabla \cdot D (r) {\nabla_{2} Φ}^{(4)} (r) = 0

\nabla \cdot D (r) \nabla_{1} Φ^{(1)} (r) - \frac{10}{7} \nabla \cdot D (r) {\nabla_{1} Φ}^{(2)} (r) + \frac{90}{7} \nabla \cdot D (r) {\nabla Φ}^{(3)} (r) - 10 {μ'}_{t} (r) Φ^{(3)} (r) = 0

\frac{1}{2} \nabla \cdot D (r) \nabla_{2} Φ^{(1)} (r) - \frac{5}{7} \nabla \cdot D (r) {\nabla_{2} Φ}^{(2)} (r) + \frac{45}{7} \nabla \cdot D (r) {\nabla Φ}^{(4)} (r) - 5 {μ'}_{t} (r) Φ^{(4)} (r) = 0

where $\nabla = \hat{x} \frac{\partial}{\partial x} + \hat{y} \frac{\partial}{\partial y}, \nabla_{1} = \hat{x} \frac{\partial}{\partial x} - \hat{y} \frac{\partial}{\partial y},$ and $\nabla_{2} = \hat{x} \frac{\partial}{\partial y} + \hat{y} \frac{\partial}{\partial x}$ Φ⁽¹⁾, Φ⁽²⁾, Φ⁽³⁾, and Φ⁽⁴⁾ are the first four components in the spherical harmonic expansion of the photon radiance where the first component, Φ⁽¹⁾ , 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, Φ⁽¹⁾ , i.e., -Dn̂·∇Φ⁽¹⁾ = αΦ⁽¹⁾, 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., Φ⁽²⁾ = Φ⁽³⁾ = and Φ⁽⁴⁾ = 0.

Making use of finite element discretizations, the discretized forms of Eqs. (1)–(4) can be written as

\sum_{j = 1}^{N} {Φ_{j}^{(1)} [〈 - \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla ϕ_{j} \cdot \nabla ϕ_{i} 〉 - 〈 \sum_{q = 1}^{Q} μ_{p} ϕ_{p} ϕ_{j} ϕ_{i} 〉] + Φ_{j}^{(2)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla ϕ_{j} \cdot \nabla ϕ_{i} 〉 - {6 Φ}_{j}^{(3)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{1} ϕ_{j} \cdot \nabla ϕ_{i} 〉 - 6 Φ_{j}^{(4)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{2} ϕ_{j} \cdot \nabla ϕ_{i} 〉}

\sum_{j = 1}^{N} {Φ_{j}^{(1)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla ϕ_{j} \cdot \nabla ϕ_{i} 〉 - \frac{25}{7} Φ_{j}^{(2)} [〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla ϕ_{j} \cdot \nabla ϕ_{i} 〉 + \frac{5}{3} 〈 \sum_{p = 1}^{Q} D_{p}^{- 1} ϕ_{q} ϕ_{j} ϕ_{i} 〉] + \frac{60}{7} Φ_{j}^{(3)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{1} ϕ_{j} \cdot \nabla ϕ_{i} 〉 - \frac{60}{7} Φ_{j}^{(4)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{2} ϕ_{j} \cdot \nabla ϕ_{i} 〉}

\sum_{j = 1}^{N} {Φ_{j}^{(1)} 〈 - \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{1} ϕ_{j} \cdot \nabla ϕ_{i} 〉 + \frac{10}{7} Φ_{j}^{(2)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{1} ϕ_{j} \cdot \nabla ϕ_{i} 〉 - \frac{90}{7} Φ_{j}^{(3)} [〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla ϕ_{j} \cdot \nabla ϕ_{i} 〉 + \frac{10}{3} 〈 \sum_{p = 1}^{P} D_{p}^{- 1} ϕ_{p} ϕ_{j} ϕ_{i} 〉]}

\sum_{j = 1}^{N} {Φ_{j}^{(1)} 〈 - \frac{1}{2} \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{2} ϕ_{j} \cdot \nabla ϕ_{i} 〉 + \frac{5}{7} Φ_{j}^{(2)} 〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla_{2} ϕ_{j} \cdot \nabla ϕ_{i} 〉 - \frac{45}{7} Φ_{j}^{(4)} [〈 \sum_{p = 1}^{P} D_{p} ϕ_{p} \nabla ϕ_{j} \cdot \nabla ϕ_{i} 〉 + \frac{5}{3} 〈 \sum_{p = 1}^{P} D_{p}^{- 1} ϕ_{p} ϕ_{j} ϕ_{i} 〉]}

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:

(ℑ^{T} ℑ + λI) Δχ = ℑ^{T} (Φ^{o} - Φ^{c})

where

ℑ = [\begin{matrix} \frac{\partial Φ_{1}^{(1)}}{\partial D_{1}} & \frac{\partial Φ_{1}^{(1)}}{\partial D_{2}} & \dots & \frac{\partial Φ_{1}^{(1)}}{\partial D_{N}} & \frac{\partial Φ_{1}^{(1)}}{\partial μ_{a, 1}} & \frac{\partial Φ_{1}^{(1)}}{\partial μ_{a, 2}} & \dots & \frac{\partial Φ_{1}^{(1)}}{\partial μ_{a, N}} \\ \frac{\partial Φ_{2}^{(1)}}{\partial D_{1}} & \frac{\partial Φ_{2}^{(1)}}{\partial D_{2}} & \dots & \frac{\partial Φ_{2}^{(1)}}{\partial D_{N}} & \frac{\partial Φ_{2}^{(1)}}{\partial μ_{a, 1}} & \frac{\partial Φ_{2}^{(1)}}{\partial μ_{a, 2}} & \dots & \frac{\partial Φ_{2}^{(1)}}{\partial μ_{a, N}} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{\partial Φ_{M}^{(1)}}{\partial D_{1}} & \frac{\partial Φ_{M}^{(1)}}{\partial D_{2}} & \dots & \frac{\partial Φ_{M}^{(1)}}{\partial D_{N}} & \frac{\partial Φ_{M}^{(1)}}{\partial μ_{a, 1}} & \frac{\partial Φ_{M}^{(1)}}{\partial μ_{a, 2}} & \dots & \frac{\partial Φ_{M}^{(1)}}{\partial μ_{a, N}} \end{matrix}]

and ∆χ = (∆D₁, ∆D₂,…∆D_N,∆μ_a,1, ∆μ_a,2,…∆μ_a,
N)^T is the update vector for the optical property profiles. Φ^o = ( $Φ_{1}^{(1), o}$ , $Φ_{2}^{(1), o}$ ,… $Φ_{M}^{(1), o}$ )^T and Φ^c = ( $Φ_{1}^{(1), c}$ , $Φ_{2}^{(1), c}$ ,… $Φ_{M}^{(1), c}$ )^T , where $Φ_{i}^{(1), o}$ and $Φ_{i}^{(1), c}$ 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, Φ⁽¹⁾ , is used in Eq. (9) since it is the dominant component [16], 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].

Fig. 1. (a) Geometry of the test case under study; (b) reconstructed D image for the first test case; (c) reconstructed μ_a image for the first test case.

Download Full Size | PDF

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 [14]. 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.

Fig. 2. (a) Recovered D image for the second case using the third-order codes; (b) recovered μ_a image for the second case using the third-order codes; (c) recovered D image for the second case using the first-order codes; (b) recovered μ_a image for the second case using the firs-order codes;

Download Full Size | PDF

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.

4. Conclusions

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.

Acknowledgments

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.

2. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modeling and reconstruction,” Phys. Med. Biol. 42, 841–853 (1997). [CrossRef] [PubMed]

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]

5. C. L. Matson, “A diffraction tomographic model of the forward problem using diffuse photon density waves,” Optics Express 1, 6–12 (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]

8. K. D. Paulsen and H. Jiang, “Spatially-varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys. 22, 691–702 (1995). [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]

10. K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total variation minimization,” Appl. Opt. 35, 3447–3458 (1996). [CrossRef] [PubMed]

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]

Optical image reconstruction based on the third-order diffusion equations

Abstract

1. Introduction

2. Reconstruction algorithm

3. Numerical examples of reconstruction

4. Conclusions

Acknowledgments

References and links

Cited By

Figures (2)

Equations (14)

Optics Express