## Abstract

We investigated the diffuse optical image errors resulting from systematic errors in the background scattering and absorption coefficients, Gaussian noise in the measurements, and the depth at which the image is reconstructed when using a 2D linear reconstruction algorithm for a 3D object. The fourth Born perturbation approach was used to generate reflectance measurements and k-space tomography was used for the reconstruction. Our simulations using both single and dual wavelengths show large systematic errors in the absolute reconstructed absorption coefficients and corresponding hemoglobin concentrations, while the errors in the relative oxy- and deoxy- hemoglobin concentrations are acceptable. The greatest difference arises from a systematic error in the depth at which an image is reconstructed. While an absolute reconstruction of the hemoglobin concentrations can deviate by 100% for a depth error of ±1 mm, the error in the relative concentrations is less than 5%. These results demonstrate that while quantitative diffuse optical tomography is difficult, images of the relative concentrations of oxy- and deoxy-hemoglobin are accurate and robust. Other results, not presented, confirm that these findings hold for other linear reconstruction techniques (i.e. SVD and SIRT) as well as for transmission through slab geometries.

© Optical Society of America

## 1. Introduction

Diffuse optical tomography (DOT) uses diffuse light to probe and image spatial variations in macroscopic optical properties including scattering and absorption [1–3]. Considerable effort has been devoted to developing a new, inexpensive, non-invasive medical imaging technique that offers new contrast possibilities. Light transport in turbid media with considerably small absorption relative to scattering, like biological tissue, can be approximated as a diffusion process [1]. Finite element [4], Monte Carlo [5, 6], as well as the Born approximation [2, 7, 8] have been used for the light transport forward calculation. These forward solutions determine the distribution of light inside the tissue. This in turn reveals the influence that localized changes in the optical properties have on the measurements. The transformation from localized image changes to measured perturbations is given by a weight function. To deduce the optical properties inside tissue from measurements on the tissue surfaces, the inverse problem must be solved. If the weight function is assumed independent of the perturbations, then the problem is linearized and a variety of linear inverse algorithms can be used for imaging reconstruction [8, 9]. Nonlinear methods, which iteratively update the weight function based on the reconstructed image, have also been use for solving the inverse problem [10, 11]. Linear algorithms are popular because of their speed, while non-linear algorithms are superior because they reduce modeling errors. The outcomes of all algorithms eventually rely on the knowledge of the averaged background optical properties on which an image of the perturbation is reconstructed. An error in the background optical properties is likely to result in an image reconstruction error. No investigation has been made of the image reconstruction errors resulting from the uncertainty in the optical properties of background media when using a linear algorithm.

In this paper, we present absorption coefficient images of absorbing objects reconstructed from continuous-wave (CW) simulated reflectance data calculated for a single source and a two dimensional detector array with various errors in the background optical parameters as well as system noise. Then we apply two wavelengths to examine how these errors and noise affect the reconstruction of oxy- and deoxy-hemoglobin concentration images. We found that the absolute values of these parameters depend significantly on the systematic errors in the optical properties of the background, while the relative changes are in some instances independent of the background error. The forward algorithm we used for simulating the data is the n^{th} order Born approximation. The inverse solution is based on the k-space reconstruction algorithm in reflectance mode [12]. We have also performed tests using other linear reconstruction techniques (SVD [13] and SIRT [2, 14]) in reflectance mode and transmission mode through a slab. The results (not shown) are consistent with those discussed for the k-space algorithm in reflectance mode. We thus draw the general conclusion that the results presented here have general applicability to all linear algorithms and geometries.

## 2. Forward algorithm and K-space theory

The diffusion of light through a highly scattering medium on a macroscopic level can be treated as a scalar wave, called a diffuse photon density wave (DPDW), whether or not the source of light is intensity modulated [15–17]. In a medium with spatially uniform optical properties, this wave propagates unobstructed. However, the presence of spatial variations in the optical properties results in a perturbation of the DPDW, which is accurately modeled as a scattered DPDW [18]. Therefore, light that diffuses through a heterogeneous medium can be described as a superposition of the incident intensity wave and those waves scattered from local variations in the optical properties.

The first Born approximation of the photon diffusion equation indicates that the wave scattered from spatial variations in the absorption coefficient is given by

*ϕ*_{sc}
(**r**
_{d}) is the scattered photon fluence measured at the detector at position **r**
_{d}, such that what we measure at the detector is *ϕ*(**r**
_{s} ,**r**) = *ϕ*_{inc}
(**r**
_{s}, **r**)+*ϕ*_{sc}
(**r**
_{s}, **r**). *ϕ*_{inc}
(**r**
_{s}, **r**) is the incident fluence at position **r** generated in the medium by a source at position **r**
_{s}. *G*(**r**, **r**
_{d}) is the Green’s function solution of the diffusion equation for the given medium. Both *ϕ*_{inc}
(**r**
_{s}, **r**) and *G*(**r**, **r**
_{d}) are calculated given the spatially uniform, average background optical properties of the medium *μ*_{s}
′, the reduced scattering coefficient, and *μ*_{a}
, the absorption coefficient. *D* = *v*/(3*μ*_{s}
′) is the photon diffusion coefficient [19, 20], *v* is the speed of light in the medium, and *δμ*_{a}
(**r**) is the spatial variation of the absorption coefficient from the average background value ma at position **r**. The second and n^{th} order approximation are given by

Diffuse optical tomography is essentially based on reconstructing an image of *δμ*_{a}
(**r**) from measurements of *ϕ*_{sc}
(**r**
_{d}) at multiple source, **r**
_{s}, and detector, **r**
_{d}, positions. Recognizing that eq. (1) resembles the Sommerfeld - Kirchoff scalar diffraction integral, we can use “diffraction tomography” to reconstruct an image of *δμ*_{a}
(**r**) from measurements of *ϕ*_{sc}
(**r**
_{d}) in a plane [21–23].

We start by noting that for imaging applications where measurements are made in reflection mode, we can usually approximate the highly scattering medium as semi-infinite. In this case, placing the free-space boundary in the xy-plane at z = 0, using a collimated pencil beam incident on the medium at x = y = z = 0, and making measurements in the same plane with the source fixed, eq. (1) can be re-written as

where we have taken the Fourier transform with respect to x, y, x_{d}, and y_{d}, and the object function *A*(*ω*_{x}
, *ω*_{y}
, *z*), which we wish to reconstruct an image of, is given by

An analytic inversion of eq. (1) can now be obtained by assuming that the object function is zero everywhere except in a slice of width *h* at position *z*. Thus, the analytic solution for the spatially varying absorption coefficient becomes

This k-space reconstruction method is described in more detail by others [21–23]. Figure 1 shows reconstructed absorption coefficient images of an absorbing object at various depths when there is no error in the background optical properties and no noise in the measurement. Simulation details are given below in the simulations section below. We see that the image is best focused at the exact depth position of the object. This suggests that the k-space algorithm is capable of pseudo-3D reconstruction, as described in more detail in [12]. This feature, however, is not as clear when there is noise in the measurement and the background

## 3. Spectroscopy for Oxy- and Deoxy-hemoglobin

Measuring blood oxygen saturation is an important application of near-infrared spectroscopy (NIRS) and DOT. The absorption coefficient at a certain wavelength and the concentration of oxy- and deoxy-hemoglobin are related by:

where *ε*_{HbO}
and *ε*_{Hb}
are the extinction coefficients of oxy- and deoxy-hemoglobin respectively and [HbO] and [Hb] are the corresponding concentrations. Below we represent the total hemoglobin concentration as [HbT] = [HbO] + [Hb] and the oxygen saturation as SO2 = [HbO] / [HbT]. The extinction spectra of hemoglobin are published in the literature [24–26] and are shown in figure 2. Measurements of the absorption coefficient at two or more wavelengths are necessary in order to solve for [HbO] and [Hb] as shown in eq. (8a) and eq. (8b).

where Δ[*HbO*] and Δ[*Hb*] refer to the change in oxy- and deoxy- hemoglobin concentration corresponding to the reconstructed absorption coefficient, ${\delta \mu}_{a}^{\mathrm{\lambda}}$, at each wavelength. In this paper, we used wavelengths of 780 nm and 830 nm.

## 4. Simulations

In our simulations, the background is a semi-infinite turbid media with a scattering coefficient of 10 cm^{-1} and an absorption coefficient of 0.05 cm^{-1}. The object has dimensions of 0.5 by 0.5 by 0.2 cm located at a depth of 2 cm. It has the same scattering coefficient as the background but a different absorption coefficient. For our single wavelength simulation, the object has an absorption coefficient of 0.1cm^{-1}. For the two-wavelength simulation, the background has an [HbT] of 100 μmole/liter and SO2 = 80%. For the object, [HbT] = 120 μmol/liter and SO2 = 60%. We used a wavelength of 780 nm and 830 nm in the simulation. From eq. (7), the background absorption coefficient is found to be 0.078 cm^{-1} and 0.092 cm^{-1} at 780 nm and 830 nm respectively and the absorption perturbation at the two wavelengths is 0.024 cm^{-1} and 0.012 cm^{-1} respectively. The measurement geometry consists of an array of 12 by 12 detectors uniformly distributed over an area of 6 by 6 cm. The light source is located at the center of the probe.

The procedure of the simulation is as follows: First, we calculate the diffuse reflectance using the fourth order Born approximation for which the object was divided into 5 by 5 by 3 voxels. Then we add Gaussian noise to this calculation to generate measurement noise. In reconstructing a 2D image using the k-space algorithm we calculate the Green’s function using systematically incorrect scattering and absorption coefficients for the background, as well as a systematically incorrect object depth. We then determine the image errors by comparing the reconstructed image value with that of the known object.

## 5. Results and discussions

In figure 3 we plot the systematic image errors as a function of the depth position at which the image is reconstructed. It shows that the image value is very sensitive to errors in the depth: 1 mm errors in the depth will cause between a 50 and 100% error in the reconstructed absorption coefficient. Underestimating the depth results in a smaller systematic error than overestimating the depth as the negative systematic errors **cannot** be greater than -100%. These large systematic errors originate from the exponentially decaying nature of the diffuse photon density wave inside the turbid media. This suggests that the systematic image errors will have a nonlinear relationship with the errors in the background absorption and scattering coefficient.

Figure 4 shows the systematic image errors resulting from the errors in the background absorption and scattering coefficients. We can see that 20% errors in the background absorption results in an ~25% error in the image value while 20% errors in the background scattering coefficient causes an ~15% error in the image value. We also notice that an underestimation of the background optical properties gives a smaller error than the same overestimation for the same reason as underestimating the depth.

Figure 5 shows the uncertainties in the image value caused by measurement noise. In figure 5a, we plot the image uncertainties as a function of the Gaussian noise magnitude, which is defined as the ratio of the standard deviation of the noise with the measured total reflectance. It indicates that 10^{-5}% noise will cause ~5% uncertainty in the image value. In figure 5b, the Gaussian noise is fixed at 10^{-5}% and the object depth is varied. We see that the image uncertainty rises significantly as the depth of the object increases. This means that imaging a deep absorbing object is more prone to uncertainty than imaging a shallow object. Although this result is obvious, it does give us a method for designing the measurement geometry to optimize the depth sensitivity for a given application. Note that for noise levels as small as 10^{-5}% we obtain image uncertainties greater than 10%. Noise levels in practice are expected to be on the order of 0.1%. This extreme sensitivity to noise in our example arises for two reasons: 1) the object contrast is small (<25%) and 2) the measurement geometry is not optimal.

The above results show that the reconstructed image value is sensitive to the errors in the background optical properties, the depth at which the 2D reconstruction is performed, and the measurement noise. These results confirm the difficulty of quantitative diffuse optical tomography. In the following, we demonstrate that it is possible to reduce the image errors by reconstructing the relative absorption changes. A practical example is the reconstruction of perturbations in the concentrations of oxy- and deoxy-hemoglobin.

Figure 6 shows errors in the reconstructed Δ[Hb], Δ[HbO] and Δ[Hb]/Δ[HbO] as a function of the depth at which the images are reconstructed. It indicates that the errors in the reconstructed Δ[Hb] and Δ[HbO] are highly sensitive to the reconstruction depth (errors ranging from -100% to +1000%), but that the errors in the reconstructed ratio Δ[Hb]/Δ[HbO] remain less than 10% given a ±2 mm uncertainty in the depth. This reduction in the systematic error in the reconstructed image results from error cancellation in the ratio, that is the systematic errors in Δ[Hb] and Δ[HbO] are nearly the same. Δ[Hb] and Δ[HbO] errors due to the errors in the background optical properties are summarized in the next section.

## 6. Summary

Table 1 and Table 2 summarize the numerical results of our one-wavelength and two-wavelength simulations respectively. These results demonstrate that quantitative image reconstruction using linear algorithms need accurate information about the background optical properties and especially the depth of the object. Small noise in the measurement can cause large uncertainties in the image values and these uncertainties increase significantly as the object depth increases. The reconstruction of the relative ratio Δ[*Hb*]/Δ[*HbO*] is much less sensitive to the depth of the object than the absolute measures of Δ[*Hb*] and Δ[*HbO*]. The ratio is also less sensitive to errors in the background optical properties. This distinction between absolute and relative imaging will play an important role in the application of diffuse optical tomography to studying hemodynamics. Ultimately, these results motivate the development of 3D reconstruction algorithms to better localize object position and to better determine background optical properties in order to minimize the systematic errors addressed here.

## 7. Acknowledgements

The authors gratefully acknowledge funding from NIH R29NS37053 and from CIMIT, the Center for Innovative Minimally Invasive Therapies. This research was funded in part by the US Army, under Cooperative Agreement No. DAMD17-99-2-9001. This publication does not necessarily reflect theposition or the policy of the Government, and no official endorsement should be inferred. We thank Tom Gaudette for his assistance in preparing the manuscript.

## References and links

**1. **A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today **48**, 34–40 (1995). [CrossRef]

**2. **M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett. **20**, 426–428 (1995). [CrossRef]

**3. **“Trends in Optics and Photonics Series,” R. Alfano ed., Advances in Optical Imaging and Photon Migration,Orlando, FLA, (OSA,1996).

**4. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med.Phys. **20**, 299–309 (1993). [CrossRef] [PubMed]

**5. **R. Graaff, M. H. Koelink, F. F. M. de Mul, W. G. Zijlstra, A. C. M. Dassel, and J. G. Aarnoudse, “Condensed Monte Carlo simulations for the description of light transport,” Appl. Opt. **32**, 426–434 (1993). [CrossRef] [PubMed]

**6. **L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,”Comput. Methods Programs Biomed. **47**, 131–146 (1995). [CrossRef] [PubMed]

**7. **S. R. Arridge, J. P. Kaltenbach, R. L. Barbour, and G. Muller ed., Bellingham, Wa, (Proc. SPIE,1993).

**8. **S. R. Arridge and J. C. Hebden, “Optical Imaging in Medicine: II. Modelling and reconstruction,” Phys. Med. Biol. **42**, 841–854 (1997). [CrossRef] [PubMed]

**9. **A. J. Devaney, “Reconstruction tomography with diffractive wave-fields,” Inverse Problems **2**, 161–183 (1986). [CrossRef]

**10. **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. **A. D. Klose and A. H. Hielscher, “A transport-theory based reconstruction algorithm for optical tomography,” B. Chance, R. Alfano, and B. Tromberg ed., SPIE BiOS99, San Jose, CA, (SPIE,1999). [CrossRef]

**12. **X. Cheng and D. A. Boas, “Diffuse optical reflectance tomography with continuous-wave illumination,” Opt. Express **3**, 118–123 (1998)http://epubs.osa.org/oearchive/source/5663.htm. [CrossRef] [PubMed]

**13. **W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *Numerical Recipes in C: The Art of Scientific Computing* (Cambirdge Univ. Press, New York,1988) Ch2 p52.

**14. **A. C. Kak and M. Slaney, *Principles of Computerized Tomographic Imaging* (IEEE Press, New York,1988).

**15. **M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Refraction of diffuse photon density waves,” Phys. Rev. Lett. **69**, 2658–2661 (1992). [CrossRef]

**16. **J. M. Schmitt, A. Knuttel, and J. R. Knutson, “Interference of diffusive light waves,” J.Opt.Soc.Am.A **9**, 1832 (1992). [CrossRef] [PubMed]

**17. **J. B. Fishkin and E. Gratton, “Propagation of photon density waves in strongly scattering media containing an absorbing ‘semi-infinite’ plane bounded by a straight edge,” J. Opt. Soc. Am. A **10**, 127–140 (1993). [CrossRef] [PubMed]

**18. **D. A. Boas, M. A. O’Leary, B. Chance, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneties within turbid media: analytic solution and applications,” Proc. Natl. Acad. Sci. USA **91**, 4887–4891 (1994). [CrossRef] [PubMed]

**19. **K. Furutsu, “Diffusion equation derived from the space-time transport equation in anisotropic random media,” J.Math.Phys. **24**, 765–777 (1997).

**20. **T. Durduran, B. Chance, A. G. Yodh, and D. A. Boas, “Does the photon diffusion coefficient depend on absorption?,” J. Opt. Soc. Am. A **14**, 3358–3365 (1997). [CrossRef]

**21. **C. L. Matson, “A diffraction tomographic model of the forward problem using diffuse photon density waves,” Opt. Express **1**, 6–11 (1997).http://epubs.osa.org/oearchive/source/1884.htm [CrossRef] [PubMed]

**22. **X. D. Li, T. Durduran, A. G. Yodh, B. Chance, and D. N. Pattanayak, “Diffraction tomography for biochemical imaging with diffuse-photon density waves,” Opt. Lett. **22**, 573–575 (1997). [CrossRef] [PubMed]

**23. **X. D. Li, T. Durduran, A. G. Yodh, B. Chance, and D. N. Pattanayak, “Diffraction tomography for biomedical imaging with diffuse photon density waves: errata,” Opt. Lett. **22**, 1198 (1997). [CrossRef] [PubMed]

**24. **M. Cope, “The Development of a Near-Infrared Spectroscopy System and Its Application for Noninvasive Monitoring of Cerebral Blood and Tissue Oxygenation in the Newborn Infant,” University College London (1991).

**25. **S. Wray, M. Cope, and D. T. Delpy, “Characteristics of the near infrared absorption spectra of cytochrome aa3 and hemoglobin for the noninvasive monitoring of cerebral oxygenation.,” Biochim Biophys Acta **933**, 184–192 (1988). [CrossRef] [PubMed]

**26. **S. J. Matcher, C. E. Elwell, C. E. Cooper, M. Cope, and D. T. Delpy, “Performance comparison of several published tissue near-infrared spectroscopy algorithms,” Anal. Biochem. **227**, 54–68 (1995). [CrossRef] [PubMed]