## Abstract

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

## 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]:

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}$ Φ^{(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

$$=-\u3008S{\varphi}_{i}\u3009-\oint D\hat{n}\xb7\nabla {\Phi}^{\left(1\right)}{\varphi}_{i}\mathrm{ds}+\oint D\hat{n}\xb7\nabla {\Phi}^{\left(2\right)}{\varphi}_{i}\mathrm{ds}-6\oint D\hat{n}\xb7{\nabla}_{1}{\Phi}^{\left(3\right)}{\varphi}_{i}\mathrm{ds}-6\oint D\hat{n}\xb7{\nabla}_{2}{\Phi}^{\left(4\right)}{\varphi}_{i}\mathrm{ds}$$

$$=\oint D\hat{n}\xb7\nabla {\Phi}^{\left(1\right)}{\varphi}_{i}\mathrm{ds}+\frac{25}{7}\oint D\hat{n}\xb7\nabla {\Phi}^{\left(2\right)}{\varphi}_{i}\mathrm{ds}-\frac{60}{7}\oint D\hat{n}\xb7{\nabla}_{1}{\Phi}^{\left(3\right)}{\varphi}_{i}\mathrm{ds}-\frac{60}{7}\oint D\hat{n}\xb7{\nabla}_{2}{\Phi}^{\left(4\right)}{\varphi}_{i}\mathrm{ds}$$

$$=-\oint D\hat{n}\xb7{\nabla}_{1}{\Phi}^{\left(1\right)}{\varphi}_{i}\mathrm{ds}+\frac{10}{7}\oint D\hat{n}\xb7{\nabla}_{1}{\Phi}^{\left(2\right)}{\varphi}_{i}\mathrm{ds}-\frac{90}{7}\oint D\hat{n}\xb7\nabla {\Phi}^{\left(3\right)}{\varphi}_{i}\mathrm{ds}$$

$$=-\frac{1}{2}\oint D\hat{n}\xb7{\nabla}_{2}{\Phi}^{\left(1\right)}{\varphi}_{i}\mathrm{ds}+\frac{5}{7}\oint D\hat{n}\xb7{\nabla}_{2}{\Phi}^{\left(2\right)}{\varphi}_{i}\mathrm{ds}-\frac{45}{7}\oint D\hat{n}\xb7\nabla {\Phi}^{\left(4\right)}{\varphi}_{i}\mathrm{ds}$$

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:

where

and ∆χ = (∆D_{1},
∆D_{2},…∆D_{N},∆μ_{a,1},
∆μ_{a,2},…∆μ_{a,
N})^{T} is the update vector for the optical property profiles.
Φ^{o} = (${\mathrm{\Phi}}_{1}^{\left(1\right),\mathrm{o}}$
,${\mathrm{\Phi}}_{2}^{\left(1\right),\mathrm{o}}$
,…${\mathrm{\Phi}}_{\mathrm{M}}^{\left(1\right),\mathrm{o}}$)^{T} and
Φ^{c} = (${\mathrm{\Phi}}_{1}^{\left(1\right),\mathrm{c}}$ ,
${\mathrm{\Phi}}_{2}^{\left(1\right),\mathrm{c}}$
,…${\mathrm{\Phi}}_{\mathrm{M}}^{\left(1\right),\mathrm{c}}$)^{T} , where
${\mathrm{\Phi}}_{\mathrm{i}}^{\left(1\right),\mathrm{o}}$ and ${\mathrm{\Phi}}_{\mathrm{i}}^{\left(1\right),\mathrm{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, Φ^{(1)} , 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].

## 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.

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. **41**767–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]