## Abstract

In order for diffuse optical tomography to realize its potential of obtaining quantitative images of spatially varying optical properties within random media, several potential experimental systematic errors must be overcome. One of these errors is the calibration of the emitter strength and detector efficiency/gain. While in principle these parameters can be determined accurately prior to an imaging experiment, slight fluctuations will cause significant image artifacts. For this reason, it is necessary to consider including their calibration as part of the inverse problem for image reconstruction. In this paper, we show that this can be done successfully in a linear reconstruction model with simulated continuous-wave data. The technique is general for frequency and time domain data.

© Optical Society of America

## 1. Introduction

During the past couple of years, research in diffuse optical tomography has rapidly progressed in moving the imaging modality from the development realm of computer simulations and phantom experiments to the application realm of imaging animal and human subjects [1–5]. This transition is forcing us to look more closely at experimental details that degrade image quality, such as the treatment of boundary conditions for solving the forward problem, source and detector (optodes) coupling to the tissue, optode positional uncertainties, and intrinsic tissue heterogeneity to name a few. There has been much work on validating and accelerating forward models and developing regularized inversion algorithms for improving quantitative image measures [6], but little work has discussed the importance and solution of experimental details.

The reconstruction of spatially varying optical properties is an ill-posed non-linear inverse problem. A linear approximation is obtained if it is assumed that the variations can be described as small perturbations, *δµ _{a}* and

*δµ*’, around a known background

_{s}*µ*and

_{a,0}*µ*’. Thence the measured fluence Φ is given by a perturbation Φ

_{s,0}_{pert}on a background fluence Φ

_{o}. In this approach therefore, we first experimentally measure the photon fluence Φ and then must calculate the fluence perturbation Φ

_{pert}caused by spatial variation, (

*δµ*,

_{a}*δµ*’) in the optical properties. Accurate determination of Φ

_{s}_{pert}requires calibration of the source and detector amplitude factors (

*s*and

*d*respectively) and a reasonable estimate of Φ

_{o}. A few approaches have been discussed for obtaining a reasonable estimate of Φ

_{o}including: 1) solving for the background optical properties using all of the measurements by first assuming zero spatial perturbations, and 2) including this estimate into the non-linear inverse problem [7–9]. In an infinite domain problem, and assuming homogeneous background, O’Leary

*et al*. have proposed using a difference of equally distant measurements for which Φ

_{o}is equal and thus cancels [10]. Thus, the difference of Φ

_{pert}between the two measurements is obtained accurately (assuming proper

*s*and

*d*calibration) without reference to Φ

_{o}.

Errors in the calibration of *s* and *d* are likely to be compensated in the image reconstruction through highly localized absorption and scattering perturbations appearing adjacent to the optodes. This type of image noise looks like high frequency spikes appearing preferentially near the optode locations, and has been observed by a number of groups [11]. A number of schemes for addressing this problem have been proposed, including median filtering [12] and the use of a regularization scheme that penalizes the reconstruction closer to the boundary with an exponentially varying weight [11, 13]. Pogue *et al*. used the same variable regularization method to suppress variation near the optodes to help improve image quality of objects far from the optodes [14]. However, this approach has the undesired effect of reducing image quality for objects near the optodes, i.e. the reconstruction produces a bias in the reconstructed object. Statistically speaking it is tantamount to assuming strong prior information in conflict to the evidence available from the data.

In this paper, we describe how calibration of *s* and *d* can be included in the image reconstruction algorithm. When considering the log of the measured fluence as the data the perturbation to the photon fluence is linearly dependent on the logarithm of *s* and *d*, and nonlinearly dependent on the optical properties. When in addition we assume linear dependence on perturbation in the optical properties (the Rytov approximation) we arrive at a fully linear problem. We present simulation results for CW data that reveal significant image improvement despite uncertainties greater than 50% and biases greater than a factor of 10. This method can also be adapted for frequency and time domain data.

## 2. Theory

#### 2.1 Forward Problem

A rigorous theory for the migration of photons through tissue has been developed based on the radiative transport equation. This approach recognizes that near-infrared photons in highly scattering media essentially undergo a random walk because the scattering probability is much greater than the absorption probability, and therefore their propagation through tissue can be described by a diffusion approximation to the transport equation. The photon diffusion equation is [15–17],

Φ(**r**,*t*) is the photon fluence at position **r** and time t. *S*(**r**,*t*) is the source distribution of photons. *D*=*v*/(*3µ _{s}*’+

*α µ*) is the photon diffusion coefficient [18, 19],

_{a}*µ*’ is the reduced scattering coefficient,

_{s}*µ*is the absorption coefficient, and

_{a}*v*is the speed of light in the medium. The dependence of D on

*µ*is expressed through the coefficient

_{a}*α*which is variously considered as 3, zero [18, 19] or some other values [20, 21]. In this paper we use zero, although the point is moot since we only look at absorption variations.

When the optical properties are spatially varying, there are two standard approaches to finding approximate solutions: the first Born approximation [15, 22]

and the Rytov approximation [22]

Each approach decomposes the total fluence Φ into Φ_{o} which only depends on the background optical properties *µ _{ao}* and

*µ*’, and Φ

_{so}_{pert}which is linearly related to spatial variations in the optical properties

*δµ*and

_{a}*δµ*’. Historically, the Born approximation is more common, but for the diffuse optical tomography problem the Rytov approximation tends to be more accurate as it accounts for some non-linear saturation due to increasing perturbation in the absorption coefficient [10]. For the Rytov approximation, assuming absorption variations only,

_{s}**r**
_{s} and **r**
_{d} are the position of the source and detector respectively, *G* is the Greens function of the photon diffusion equation for the background optical properties given the boundary conditions. If the background is homogeneous then *G* can be expressed analytically in some simple geometries[16, 23].

In this paper we consider steady state transmission measurements through a slab, i.e. the sources are on one side and the detectors are on the other. We use the extrapolated zero-boundary condition and the method of images to solve for the Greens function, G, of the diffusion equation [17]. What we measure is related to the Greens function by the source power and coupling into the medium as well as the detector gain factors and coupling from the medium. We express these factors as *s* and *d* respectively and therefore

#### 2.2 Inverse Problem

Given experimental measurements of Φ, we can reconstruct images of the spatially varying optical properties by solving a least-squares problem. For the Rytov approximation we minimize

where the index *i* is summed over measurements with each source and detector pair, and Φ* _{Theory}*(

*x*) is given by eq. (3) and eq. (4) where

*x*is a vector giving

*δµ*at each voxel position.

_{a}For the case of fewer measurements than unknowns the linear problem is underdetermined and is given by the (regularized) Moore-Penrose generalized inverse

where *I* is the identity matrix, *λ* is the Tikhonov regularization parameter, and *y* is the measured data. Each element of the matrix A is ${A}_{i,j}=-\frac{v}{{D}_{o}{\Phi}_{o}({\mathbf{r}}_{s,i},{\mathbf{r}}_{d,i})}{\Phi}_{o}({\mathbf{r}}_{s,i},{\mathbf{r}}_{j})G({\mathbf{r}}_{j},{\mathbf{r}}_{d,i})$ where **r**
_{s,i} and **r**
_{d,i} is the position of the i^{th} source and detector and **r**
_{j} is the position of the j^{th} voxel. Within the Rytov approximation ${y}_{i}=\mathrm{ln}\left[\frac{\Phi ({\mathbf{r}}_{s,i},{\mathbf{r}}_{d,i})}{{\Phi}_{o}({\mathbf{r}}_{s,i},{\mathbf{r}}_{d,i})}\right]$. In the results given here, we chose *λ*=10^{-3} of the maximum value within *AA ^{T}* [24, 25].

#### 2.3 Image and Coupling Coefficient Reconstruction — Rytov Approximation

As discussed in the introduction[DAB1], an accurate estimate of **A** and **y** requires accurate knowledge of the background optical properties *µ _{ao}* and

*µ*’ and of the

_{so}*s*and

*d*amplitudes. Here we consider

*µ*and

_{ao}*µ*’ known, and the

_{so}*s*and

*d*amplitude for each source and detector as unknowns in the inverse problem.

In this augmented model, *y _{i}* written in terms of the unknown source coupling amplitudes

*s*, detector coupling amplitudes

_{k}*d*, and the absorption perturbations

_{l}*δµ*becomes

_{a,j}*k*(*i*) and *l*(*i*) identifies the source and detector respectively used for the *i ^{th}* measurement. Note that the

*s*and

*d*appearing in eq. (8) represent the difference from the assumed values as used in Φ

*and the real values for the measurement of Φ. We can write this in matrix form as*

_{o}**y**=

**Bξ**where

**B**=[

**Ã S D**] and

*N _{v}* is the number of voxels,

*N*is the number of sources, and

_{s}*N*is the number of detectors. Scaling

_{d}*δµ*by

_{a,j}*µ*makes the elements dimensionless and of the same order as ln

_{ao}*s*and ln

*d*.

**Ã**=

*µ*

_{ao}**A**is the rescaling of the matrix given in eq. (4), and

**S**and

**D**have block diagonal form.

**S**has a one in the

*k*column for each measurement corresponding to source k, and

^{th}**D**has a one in the

*l*column for each measurement corresponding to detector

^{th}*l*. For instance, with 4 measurements between 2 sources and 2 detectors,

## 3. Simulation Results

For the simulations, we considered transmission through a 6 cm thick slab with background optical properties of *µm _{so}*’=10 cm

^{-1}and

*µ*=0.05 cm

_{ao}^{-1}(see fig. 1). A 1.6 cm diameter absorbing object with

*µ*’=10 cm

_{so}^{-1}and

*µ*=0.15 cm

_{a}^{-1}was centered at (x,y,z)=(1, -1, 3) cm. Sixteen sources were placed in a 4×4 grid at z=0, spanning x=-3 to 3 cm and y=-3 to 3 cm in 2 cm steps. Sixteen detectors were placed in a 4×4 grid at z=6 cm spanning x=-3 to 3 cm and y=-3 to 3 cm in 2 cm steps. All simulated measurements were made with continuous-wave sources.

The Full Born expansion was used to solve the forward problem as it is known to converge to the exact solution [6]. A 1.6×1.6×1.6 cm cube divided into 7×7×7 voxels was centered over the 1.6 cm diameter absorbing sphere to define the region-of-interest for the forward problem. The source and detector amplitudes were chosen randomly from a normal distribution with a mean of 1. The effect of different standard deviations (0%, 40%, and 80%) was investigated. The normal distribution was biased by discarding negative amplitudes. A separate instance of the source and detector amplitude variation was considered for each standard deviation. There was no additive measurement noise (i.e. shot or detector electronic noise) in the simulated data, just the model error associated with the source and detector amplitudes.

For the inverse problem, voxel centers were distributed from x=-3 to 3 and y=-3 to 3 in 0.5 cm steps and z=0.5 to 5.5 in 0.5 cm steps. There were 1859 voxels each with dimensions of 0.5×0.5×0.5 cm. The vector in eq. (9) was initialized to zeros, which is equivalent to an initial assumption of 1 for the source and detector amplitudes. Two different inverse problems were considered: (1) Rytov without coupling, and (2) Rytov with coupling. In each case, Tikhonov regularization was used to obtain the psuedo-inverse since the system matrix is under-determined and ill-conditioned.

In figure 2 we plot the results obtained with the Rytov approximation without and with reconstruction of the coupling coefficients. The results in figure 2 a-c show the effect of neglecting systematic errors in the coupling coefficients. The color scale is linear and ranges over -0.1(blue)≤*δµ _{a}*≤0.1(red) cm

^{-1}. With no variance (i.e. no systematic errors) the object is localized to the correct voxel in X and Y and is within 1 voxel (0.5 cm) in Z. The amplitude of the object is slightly underestimated because of blurring of the object. We reconstruct

*δµ*=0.098 cm

_{a}^{-1}while we expect

*δµ*=0.010 cm

_{a}^{-1}. Notice the image artefact near the sources and detectors. We do not observe this when the first Born approximation is used for the forward problem, and therefore conclude that it arises from mismatch between the exact forward solution and the Rytov approximation used for the inverse model (i.e. it results from model noise unrelated to the source and detector coupling coefficients).

As the standard deviation in the coupling coefficients is increased, it is possible to still make out the object, but it has reduced contrast relative to surrounding voxels and its amplitude is smaller than that in voxels closest to the sources and detectors. Note that the artefact amplitude strongly exceeds the color scale and has been truncated to preserve the sensitivity of the color scale to reveal the object of interest. The model noise caused by incorrect calibration of the source and detector amplitudes (i.e. systematic errors in *s* and *d*) is clearly degrading image quality by increasing artefact.

Next, consider Fig 2 d-f which shows the reconstruction results that simultaneously determines the source and detector coupling coefficients. The color scale is linear and ranges over -0.05(blue) ≤*µ _{a}*≤0.05(red) cm

^{-1}. When there is no variance in the coupling coefficients, we see that the location of the object is determined properly in X and Y and is within 0.5 cm in Z. The resolution (blurring or point spread function) is poorer relative to the case without simultaneous reconstruction of the coupling coefficients. This is because we now have more unknowns in the reconstruction being considered. Furthermore, notice that the artefact observed previously in fig. 2a due to error between the exact and Rytov approximation to the forward problem, has been reduced. We therefore conclude that reconstruction for the coupling coefficients has slightly compensated for model mismatch in the Rytov approximation. This is a qualitative comparison between the reconstruction without and with simultaneous determination of the coupling coefficients. A more quantitative comparison would closely examine the singular value spectrum and chose the regularization parameter based on the signal-to-noise ratio in the measurement.

Examining the images reconstructed with increasing variance in the coupling coefficients, relative to the same case in fig. 2 a-c, we see that the image quality is significantly improved. The images, in fact, are as good as the case with no variance. It is hard to notice differences in the images given this color scale, but they do exist. It is important to note that reconstruction does not accurately determine the source and detector coefficients, but does accurately determine the source detector coefficient product (e.g. *s _{k} d_{l}*) for each measurement. For the examples shown here we see that the products are reconstructed with an accuracy of a few percent despite uncertainties exceeding 80%. Finally, when we include additive shot or electronic measurement noise in our simulations, we then observe a greater sensitivity to the measurement noise as the source and detector variance increases (results not shown).

## 4. Conclusions

We have proposed an augmented model for the measurement of light in turbid media that assumes the diffusion model for propagation together with coupling amplitudes at source and detector optodes. These coupling coefficients represent additional unknowns in the reconstruction problem. By taking the log of measured fluence as the data, these coupling coefficients are linearly related to the measured data, and by assuming the Rytov approximation for small perturbations a linear problem for image reconstruction and coefficient recovery was formulated. The method can also be adapted for the Born Approximation but leads to non-linear dependence on the coupling coefficients through the product of *s* and *d*.

We illustrated the method using simulated CW data and absorption recovery only, for small perturbations. However the method is easily generalized to time and frequency domain data and to the recovery of absorption and scattering. In addition we can formulate the problem within a non-linear image reconstruction scheme as well. These extensions will be reported in a later paper. Preliminary experimental results have also been successful, and will be reported.

## Acknowledgements

We acknowledge the Wellcome Trust for a collaborative travel grant that facilitated this research. DAB acknowledges financial support from Advanced Research Technologies, NIH R29-NS38842, NIH P41-RR14075 and from 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 the position or the policy of the Government, and no official endorsement should be inferred.

## References and links

**1. **S. B. Colak, M. B. van der Mark, G. W. Hooft, J. H. Hoogenraad, E. S. van der Linden, and F. A. Kuijpers,“Clinical Optical Tomography and NIR Spectroscopy for Breast Cancer Detection,” IEEE Journal of Selected Topics in Quantum Electronics **5**, 1143–1158 (1999). [CrossRef]

**2. **S. Fantini, M. A. Franceschini, E. Gratton, D. Hueber, W. Rosenfeld, D. Maulik, P. G. Stubblefield, and M. R. Stankoivic,“Non-invasive optical imaging of the piglet brain in real time,” Opt. Express **4**, 308–314 (1999). [CrossRef] [PubMed]

**3. **V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance,“Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc Natl Acad Sci U S A **97**, 2767–72 (2000). [CrossRef] [PubMed]

**4. **M. A. Franceschini, V. Toronov, M. Filiaci, E. Gratton, and S. Fanini,“On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Express **6**, 49–57 (2000). [CrossRef] [PubMed]

**5. **B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Wells, K. S. Osterman, and U. L. Osterberg,“Hemoglobin imaging of breast tumors with near-infrared tomography,” Radiology214, (in press).

**6. **S. R. Arridge,“Optical Tomography in medical imaging,” Inverse Problems **15**, R41–R93 (1999). [CrossRef]

**7. **J. C. Hebden, E. W. Schmidt, M. E. Fry, M. Schweiger, E. M. C. Hillman, D. T. Delpy, and S. R. Arridge,“Simultaneous reconstruction of absorption and scattering images by multichannel measurement of purely temporal data,” Opt. Lett. **24**, 334–336 (1999). [CrossRef]

**8. **M. Schweiger and S. R. Arridge,“Optical tomographic reconstruction in a complex head model using a priori region boundary information,” Phys Med Biol **44**, 2703–21 (1999). [CrossRef] [PubMed]

**9. **V. Kolehmainen, M. Vauhkonen, J. P. Kaipio, and S. R. Arridge,“Recovery of Piecewise constant coefficients in optical diffusion tomography,” Opt. Express **7**:468–480 (2000). [CrossRef] [PubMed]

**10. **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]

**11. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “Performance of an iterative reconstruction algorithm for near infrared absorption and scatter imaging,” in *Photon Migration and Imaging in Random Media and Tissues*, B. Chance and R. R. Alfano, SPIE1888, 360–371 (1993).

**12. **M. Schweiger, S. R. Arridge, and D. T. Delpy,“Application of the finite-element method for the forward and inverse models in optical tomography,” Journal of Mathematical Imaging and Vision **3**, 263–283 (1993). [CrossRef]

**13. **S. R. Arridge and M. Schweiger, “Inverse Methods for Optical Tomography” in *Information Processing in Medical Imaging (IPMI’93 Proceedings), Lecture Notes in Computer Science*, (Springer-Verlag, 1993).

**14. **B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Osterberg, and K. D. Paulsen,“Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. **38**, 2950–2961 (1999). [CrossRef]

**15. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Academic Press, Inc., San Diego1978).

**16. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28**, 2331–2336 (1989). [CrossRef] [PubMed]

**17. **R. C. Haskell, L. O. Svaasand, T. Tsay, T. Feng, M. S. McAdams, and B. J. Tromberg,“Boundary conditions for the diffusion equation in radiative transfer,” J.Opt.Soc.of Am.A **11**, 2727–2741 (1994). [CrossRef]

**18. **K. Furutsu and Y. Yamada,“Diffusion approximation for a dissipative random medium and the applications,” Phys.Rev.E **50**, 3634 (1994). [CrossRef]

**19. **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]

**20. **D. J. Durian, “The diffusion coefficient depends on absorption,” Optics Letters **23**, 1502–1504 (1998). [CrossRef]

**21. **R. Aronson and N. Corngold, “Photon diffusion coefficient in an absorbing medium,” J. Opt. Soc. Am. A **16**, 1066–1071 (1999). [CrossRef]

**22. **A. C. Kak and M. Slaney, *Principles of Computerized Tomographic Imaging* (IEEE Press, New York1988).

**23. **S. R. Arridge, M. Cope, and D. T. Delpy,“The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys Med Biol **37**, 1531–60 (1992). [CrossRef] [PubMed]

**24. **H. Dehghani, D. C. Barber, and I. Basarab-Horwath, “Incorporating a priori anatomical information into image reconstruction in electrical impedance tomography,” Physiol Meas **20**, 87–102 (1999). [CrossRef] [PubMed]

**25. **V. Kolehmainen, M. Vauhkonen, P. A. Karjalainen, and J. P. Kaipio,“Assessment of errors in static electrical impedance tomography with adjacent and trigonometric current patterns,” Physiol Meas **18**, 289–303 (1997). [CrossRef] [PubMed]