1. Introduction

The use of contrast agents in Diffuse Optical Tomography[1] (DOT) for disease diagnostics and for probing tissue functionality follows established clinical imaging modalities such as Magnetic Resonance Imaging (MRI) [2,3,4], Ultrasound [4,5,6] and X-ray Computed Tomography[4,7,8] (CT). Contrast agent administration provides for accurate difference images of the same heterogenous tissue volume under nearly identical experimental conditions. This approach yields superior diagnostic information. In MRI of breast cancer for example, it is the relative signal enhancement compared to the baseline image, and the architectural features of this enhancement that aid in characterizing lesions as benign or malignant [9,10]. While architectural features are unlikely to improve breast cancer diagnosis using DOT in current implementations, due to the low resolution achieved [11,12], the relative signal enhancement after contrast agent administration could significantly increase sensitivity and specificity.

The theory presented in this paper is motivated by clinical DOT measurements of the human breast [13]. In these measurements optical data are obtained before and after intravenous administration of the optical contrast agent Indocyanine Green (ICG). ICG is an absorber and a fluorophore in the Near Infrared (NIR); here we focus on its absorption function. In principle DOT images of the breast before and after ICG administration may be reconstructed and subtracted. However, in practice a more robust approach derives images of the differential changes due to the extrinsic perturbation. In the latter case experimental measurements use the exact same geometry within a few minutes of one another, thereby minimizing positional and movement errors and instrumental drift. Furthermore the use of differential measurements eliminates systematic errors associated with the different medium required to calibrate operational parameters of the instrument or to provide a baseline measurement for independent reconstructions of the pre- and post ICG images. In addition the effect of surface absorbers such as hair or skin color variation is also minimized.

The main analytical difficulties of the differential approach arise because the media are inhomogeneous. Thus the total diffuse light field in the ICG perturbation problem does not separate into a homogeneous background field and a scattered field in a straightforward way. Furthermore the average background properties, particularly the absorption, change as a result of contrast agent administration.

In this paper we derive a perturbation solution that can be used to yield differential images of induced absorption perturbations. This is the primary effect of many molecular optical contrast media that are administered in low concentration. The solution is also applicable for imaging the differential variations in tissue scattering. In obtaining these results we will expose the full consequences of inhomogeneous tissue volumes, and identify the conditions wherein these effects are small. We demonstrate the superiority of the method compared to the typical perturbation approach using simulated data of contrast-enhanced breast. The method is also suitable for imaging differential changes as arise in measurements of functional activity such as brain activation or muscle exercise (the only requirement is that baseline measurements of the same tissue volume are taken during the procedure), or in measurements at multiple wavelengths. In the following analysis the term pre-ICG breast is used to indicate the “baseline” breast before the contrast agent injection and the term post-ICG breast marks the breast following contrast agent administration.

2. Theory

Photon transport in highly scattering low absorbing media such as tissue can be approximated as a diffusional process [14]. For heterogeneous media, the diffusion equation in the frequency domain is [15]:

Here U(r⃗) is the photon fluence [W·cm^-2], ω is the source modulation frequency, c is the speed of light in the medium [cm·s^-1], D(r⃗)≈c/3µ′_s(r⃗) is the medium diffusion coefficient [cm²·s^-1], µ′_s(r⃗) the medium reduced scattering coefficient [cm^-1], µ′_a(r⃗) the medium absorption coefficient [cm^-1] and S(r⃗) the source term [W·cm^-3]. In this work we consider solutions of the heterogeneous diffusion equation in the frequency domain employing the perturbation approach [16,17].

The first order perturbation expansion divides the absorption (µ′_a(r⃗)) and diffusion (D′(r⃗)) coefficients of the pre-ICG breast into spatially varying (δµ′_a(r⃗),δD′(r⃗)) and background components (µ′_a0,D′₀), i.e. µ′_a(r⃗)=µ′_a0+δµ′_a(r⃗) and D′(r⃗)=D′₀+δD′(r⃗). Throughout this paper a single ′ denotes pre-ICG tissue volumes. In the Rytov approximation [16] the total photon density wave measured at position r⃗_d due to a source at position r⃗_s is written as the product of two components, i.e.

where the scattered field ϕ′_sc(r⃗_s,r⃗_d,ω), is produced by heterogeneities (δµ′_a(r⃗),δD′(r⃗) and the incident field U′₀(r⃗_s,r⃗_d,ω), is the field that would have been detected from the same medium if these heterogeneities were not present. These terms are incorporated into the diffusion equation, whose formal solution can be expressed as an integral equation using the appropriate Green’s function for the geometry implemented. For simplicity we present this analysis for an infinite medium. This theory however can be easily extended to other simple geometries such as semi-infinite or slab, using weights derived with the method of image sources [18] and the appropriate extrapolated boundary condition [19, 20, 21].

The first order perturbative solution of the heterogeneous diffusion equation yields

where W′_a,(W′_s) represents the absorption (scattering) weight of the voxel at position r⃗, due to a source at r⃗_s and a detector at r⃗_d . In the infinite medium these weights are:

Following the administration of a contrast agent the background optical properties change. The new, post-ICG, total field can be written in a similar form:

Here ϕ″_sc(r⃗_s,r⃗_d,ω) is the field component scattered from the post-ICG heterogeneities (i.e. δµ″_a(r⃗),δD″(r⃗) with respect to the new background optical properties µ″_a0,D″₀) and U″₀(r⃗_s,r⃗_d,ω) is the incident field obtained from the homogeneous background medium with µ″_a0 D″₀. The first order perturbative solution of the heterogeneous diffusion equation yields

where W_a (W_s) represents the absorption (scattering) weight of voxels at position r⃗, due to a source at r⃗_s and for a detector at r⃗_d :

Combining Eq.2 with Eq.7 we obtain the relative scattered field, i.e

We will show that ϕ _sc(r⃗_s,r⃗_d,ω) can be attributed primarily to perturbations created by the contrast agent injection. U′(r⃗_s,r⃗_d,ω) and U″(r⃗_s,r⃗_d,ω) are the actual measurements on the pre- and post- ICG breast respectively, and U′₀(r⃗_s,r⃗_d,ω), U″₀(r⃗_s,r⃗_d,ω) can be determined from the average optical properties of the pre-and post- ICG breast (see discussion section).

During the administration of an absorption contrast agent the scattering properties of tissue are not expected to change. Therefore we assume D′₀=D″₀=D ₀ and δD′(r⃗)=δD″(r⃗). Substitution of Eq.3 and Eq.8 into Eq.12 yields:

Here we have also assumed W″_s(r⃗_s,r⃗_d,r⃗,µ″_a0,D ₀,ω)≈W′_s(r⃗_s,r⃗_d,r⃗,µ′_a0,D ₀,ω). This is a very good approximation when the average absorption change due to the contrast agent is small or in the transmission geometry (see Appendix and discussion for absorption variations below).

Let ${\delta \mathit{\mu }}_a^{\mathit{\text{ICG}}}$ (r⃗) be the total absorption perturbation due to the ICG injection that includes both position-independent and position-dependent contributions. Then µ″_a(r⃗) can be written

so that

The quantity ${\delta \mathit{\mu }}_a^{\mathit{\text{rel}}}$ (r⃗)=µ′_a0-µ″_a0+${\delta \mathit{\mu }}_a^{\mathit{\text{ICG}}}$ (r⃗) represents the position-dependent absorption heterogeneities induced by the contrast agent. The relative scattered field is computed by substitution of Eq. 15 into Eq.13. It depends on contrast agent induced absorption heterogeneities and on pre-ICG tissue absorption heterogeneities.

The second integral in Eq.16,

describes the influence of the pre-existing (intrinsic) absorption heterogeneity of the breast on the relative scattered field. Since the intrinsic heterogeneity is weighted by the difference W″_a-W′_a the influence of this term can be quite small. Using the analytical forms of the weights (Eq.4 and Eq. 9) we can write out Eq.17 explicitly, i.e.

where

The term e ^{i(k′-k″)R(r⃗)} in Eq.18 is approximately unity and S≈0 when the average absorption increase due to the ICG injection is very small (i.e. k′≈k″). Usually however k′≠k″. For example the recommended ICG dosage for humans (0.25mg/kg) introduces an average µ_a increase within the interval [0.005–0.015] cm^-1 depending on breast vascularization [13]. Fig. 1a and b show the amplitude and phase of the term e ^{i(k′-k″)R(r⃗)} respectively, for different R(r⃗), as a function of the post-ICG breast absorption coefficient for a source detector separation |r⃗_d -r⃗_s |=6cm, using the geometry of Figure 1c. The background µ_a =0.05 cm^-1 and the background µ’_s=10cm^-1.

 

Figure 1. (a) Amplitude of the term e ^{i(k′-k″)R(r⃗)} as a function of the average absorption coefficient of the post-ICG breast, (b) Phase of the term e ^{i(k′-k″)R(r⃗)} as a function of the average absorption coefficient of the post-ICG breast, (c) Test geometry for calculations in (a) and (b).

Download Full Size | PDF

The deviation of e ^{i(k′-k″)R(r⃗)} from unity increases for perturbations farther from the line adjoining source and detector (i.e. as |r⃗_d -r⃗|+|r⃗_s -r⃗| grows larger than |r⃗_d -r⃗_s | when a increases). However, the probability for photons to pass through these “distant” perturbations decreases exponentially via the weight W″_a in the integrand of Eq.18. Hence accumulated contributions of the heterogeneities at large a are small. Figure 2 plots the deviations introduced into ϕ_sc (r⃗_s,r⃗_d,ω) by taking S=0. Figure 2a depicts the ratio of the amplitude detected with no approximation to the amplitude detected assuming S=0. Similarly Fig. 2b depicts the phase shift between the phase detected with no approximation and the phase detected assuming S=0. The error is plotted for a single perturbation at different positions a for the geometry depicted in Fig. 1c. The values assumed in Eq. 16 were δµ′_a(r⃗)=0.05 cm^-1, ${\delta \mathit{\mu }}_a^{\mathit{\text{rel}}}$ (r⃗)=0.05 cm^-1 and the background optical properties µ′_a0=0.05cm^-1, µ″_a0=0.05cm^-1 and µ′_s=10cm^-1.

 

Figure 2. (a) Amplitude deviation and (b) Phase shift introduced in the field measured in Eq. 16 when S is assumed zero. The calculation is done as a function of the distance a., for the geometry shown in Fig. 1c, assuming 200MHz, background µ′_a0=0.05cm^-1, µ″_a0=0.05cm^-1, µ′_s=10 cm^-1, ${\delta \mathit{\mu }}_a^{\mathit{\text{rel}}}$ =0.05 cm^-1 and δµ′_a(r⃗)=0.01 cm^-1

Download Full Size | PDF

The simulation of Fig.2 explicitly shows that the errors introduced because of the approximation S=0 are very small for physiologically relevant optical properties (i.e. relatively small δµ′_a(r⃗)) provide the most probable photon paths. The same behavior is exhibited for the scattering weights as shown in the Appendix. Eq.16 thus becomes

Our conclusions do not change when image sources are invoked to satisfy more complex boundary conditions such as semi-infinite or slab geometries. In these cases e ^{i(k′-k″)R(r⃗)} will appear in all the terms corresponding to image sources. Note however that the assumption that S≈0 is best suited for slab geometry where |r⃗_d-r⃗|+|r⃗_s-r⃗|≈|r⃗_d-r⃗_s | for the most probable photon paths. This condition is not always true for reflectance geometry.

Notice that the differential measurements are insensitive to surface artifacts such as small skin absorbers and hair under a certain source or detector. The term δµ′_a(r⃗) in Eq.18 could be used to approximate surface heterogeneities by taking r⃗ to be close to medium surface, near to the corresponding source or detector. The influence of such terms is virtually zero since in such a geometry R(r⃗)≅0 and subsequently S=0.

For image reconstruction, Eq. 20 is discretized into a sum of volume elements [22] (voxels) and the scattered field is obtained at each modulation frequency w for every source-detector pair w. For n voxels and m=o×p×q measurements, where o is the number of sources, p is the number of detectors and q is the number of frequencies employed, the discretization yields a set of coupled linear equations

Inverting the weights’ matrix determines the spatial map of absorption due to contrast agent injection.

3.Results and Discussion

Although equations 20 and 21 resemble the result of typical perturbation analyses, there are fundamental differences and constraints that must be considered when using them. First the parameter imaged is the synthetic perturbation term ${\delta \mathit{\mu }}_a^{\mathit{\text{rel}}}$ (r⃗)=µ′_a0-µ″_a0+${\delta \mathit{\mu }}_a^{\mathit{\text{ICG}}}$ (r⃗). Secondly, the relative scattered field ϕ_sc depends both on the ratio, U″/U′, of the actual pre-ICG and post-ICG measurements, and the multiplicative term U′₀/U″₀=exp(i(k′-k″)·|r⃗_s-r⃗_d|). This term expresses the change in the incident field due to the average absorption coefficient increase of the post-ICG breast. Its use in Eq.12 leads to significant reconstruction improvements [23]. Note that the term U′₀/U″₀ depends on |r⃗_s-r⃗| and not on R(r⃗). Therefore the arguments that led on the elimination of S from Eq.16 cannot be applied to this term since |r⃗_s-r⃗|≫R(r⃗). The term U′₀/U″₀ can be analytically calculated for simple geometries such as infinite, semi-infinite or slab or calculated numerically for more complicated geometries if we know the average optical properties of the pre- and post- ICG breast.

Using simulated data we have examined the performance of Eq.20 to image the contrast-enhanced breast as compared to:

1) The typical Rytov approximation which assumes ${\varphi "}_{\mathit{sc}}=\ln \left(\frac{U"}{U\text{'}}\right)$, namely

This comparison captures a critical feature of our work. The classical perturbation approach that does not consider the average absorption increase due to the extrinsic contrast.

2) Eq.17 including the term S, namely

This comparison investigates the effect of the approximation S=0 assumed in Eq.20.

3) A difference image produced by subtracting the post-ICG image from the pre-ICG image. Both pre- and post- images were produced using Eq.3, assuming a homogeneous medium as baseline. The optical properties of the homogeneous medium were µ′_a0=0.03cm^-1 and µ′_s0=8cm^-1.

In order to perform the comparisons we obtained two MRI coronal slices of a human breast: one before and one after contrast enhancement. Figure 3a depicts the T₁-weighted MR image. This image depicts structure. White regions correspond primarily to adipose (fatty) tissue while dark regions correspond to parenchymal (glandular) tissue. Figure 3b depicts the signal enhancement of the same T1-weighted image due to injection of the MRI contrast agent Gd-DTPA. The Gd-DTPA enhancement is superimposed in color. An infiltrating ductal carcinoma (shown in yellow) demonstrated the highest signal enhancement. Gd-DTPA and ICG have similar distribution patterns. Here we assume that the Gd-DTPA distribution reflects the ICG distribution.

 

Figure 3. (a) T1-weighted MR coronal slice of a human breast, (b) Gd-DTPA distribution (in color) of the same coronal slice. A ductal carcinoma appears in yellow.

Download Full Size | PDF

We converted the MR images to optical property maps, separating four structures based on the image intensity information (by applying appropriate thresholds). The cancer is assumed having two states: pre- and post- ICG contrast. The structures selected and the corresponding absolute optical properties are shown in TABLE I. The optical properties are chosen to simulate breast optical properties as obtained from our breast clinical measurements [13].

TABLE I. Absolute optical properties of the different structures used for the simulations.

View Table

Scattering has been assumed constant for all structures for simplicity. The resulting absorption maps are shown in figure 4. The medium surrounding the breast was arbitrarily simulated as a highly absorbing diffuse medium (µ_a =0.30cm^-1 µ_s ’=8cm^-1). The average absorption of the pre- and post- ICG breast were found to be µ′_a0=0.0473 cm^-1 and µ″_a0=0.0589 cm^-1 so that average absorption increase due to the ICG is µ″_a0-µ′_a0=0.0116 cm^-1.

 

Figure 4. Absorption maps used for the simulation. (a) pre-ICG breast (b) post-ICG breast.

Download Full Size | PDF

The maps of figure 3 served as an input to a finite-differences implementation of the frequency-domain diffusion equation. The simulation assumed 7 sources and 21 detectors as shown in Figure 5. The frequency employed was 200MHz. No noise was added in the forward data.

 

Figure 5. Sources and detector arrangement used for the simulation. The region reconstructed is outlined with a green double line.

Download Full Size | PDF

For reconstruction purposes, the region of interest (indicated in Fig. 5 as a green double line) was segmented into 35×25 voxels. The inversion was performed using the Algebraic Reconstruction Technique (ART) [16], with relaxation parameter λ=0.1. Convergence was assumed when an additional 100 iterations did not change the result more than 0.1%. The simulated image and the reconstructed ${\delta \mathit{\mu }}_a^{\mathit{\text{rel}}}$ (r⃗) for the three cases examined are shown in Figure 6.

 

Figure 6. Reconstruction results imaging the region indicated with the green double line on Fig.5 assuming an infinite slab geometry. a)The reconstructed image using Eq.20. b)The reconstructed image using Eq.22. c)The reconstructed image using Eq.23. d) The result of subtracting an image of the post-ICG breast (reconstructed relative to a homogeneous baseline medium) from an image of the pre-ICG breast (reconstructed relatively to the same baseline medium). The optical properties of the homogeneous medium were µ′_a0=0.03cm^-1 and µ′_s =8cm^-1.

Download Full Size | PDF

The superior performance of Eq.20 compared to the typical perturbative formulation (Eq. 22) can be evaluated by examining Fig.6a and b. Although both methods resolve the cancerous lesion with comparable positional and size accuracy, the typical formulation (Fig.6b) yields several strong artifacts close to the boundaries. These artifacts illustrate that in the presence of distributed absorption, the perturbative method converges preferentially to localized regions of high absorption. This is often true when inverting underdetermined systems. Eq.20 on the other hand removes the “average absorption increase” from the measurement vector. Therefore Fig. 6a images weaker perturbations introduced by the ICG injection, relative to the average absorption increase. Since by construction the perturbation method works especially well for weak perturbations [16], it is expected that the use of Eq. 20 will more accurately image the heterogeneous medium. The same behavior is expected for a Born-type [16] perturbative formulation. We note that Fig.6b images the ${\delta \mathit{\mu }}_a^{\mathit{\text{ICG}}}$ (r⃗) and not the ${\delta \mathit{\mu }}_a^{\mathit{\text{rel}}}$ (r⃗) as in Fig, 6a and 6c. Therefore it is reasonable that the reconstructed value for cancer in Fig. 6b is higher than the value reconstructed in Fig.6a and 6c. The difference in reconstructed values equals approximately the average absorption increase in the post-ICG breast (µ″_a0-µ′_a0=0.0116 cm^-1).

Figure 6c has been produced after correcting the measurement vector with S instead of setting it to zero as in Fig. 6a. Only minor differences exist between the two images as had been predicted in Figure 2. In this simulation the pre-ICG cancer had a contrast of 2:1 to the average pre-ICG background value. This contribution has most likely resulted in the minor differences observed between the two images, especially in the structures close to the boundaries.

Figure 6d is the result of subtracting an image of the post-ICG breast from an image of the pre-ICG breast. Here the ${\delta \mathit{\mu }}_a^{\mathit{\text{ICG}}}$ (r⃗) is imaged. The magnitude of the cancer is slightly overestimated and its size is significantly overestimated. Similarly to Fig. 6b, strong artifacts appear close to the boundary. A distributed absorption is also reconstructed which does not correspond to the ICG distribution and is also an artifact. Compared to the other approaches the subtraction yields the most artifacts.

In these simulations the average optical properties were known by simple integration over the optical property map. In our clinical implementation [13] the average optical properties of the pre-ICG breast are calculated by fitting the experimental time-domain data obtained to the appropriate solution of the diffusion equation for the geometry used. Furthermore an algorithm has been developed [24] that calculates the difference µ′_a0-µ″_a0, (necessary to calculate both U′₀/U″₀ and ${\delta \mathit{\mu }}_a^{\mathit{\text{ICG}}}$ (r⃗)) with an accuracy of the order of 10^-3 cm^-1.

To conclude we have described a formulation of perturbation theory, which is particularly well suited for image reconstructions of differences in the absorption properties of tissues as a result of optical contrast agent administration. Importantly, these results enable the extraction of differential contrast agent absorption even within media that are heterogeneous in the absence of the contrast agent. Our primary result is an intuitive equation, which is valid over a large range of conditions. We have shown explicitly what these corrections are and how these corrections can be included in more careful analyses. The results should be applicable for a broad range of other DOT applications wherein baseline and “stimulated” measurements are available, particularly functional imaging in brain and muscle.

APPENDIX

In this section we show that the ratio of the scattering weight W″_s(r⃗_s,r⃗_d,r⃗,µ″_a0,D ₀,ω) to W′_s(r⃗_s,r⃗_d,r⃗,µ′_a0,D ₀,ω) is sufficiently small so that when subtracting Eq.3 from Eq.8 to produce Eq.13, the scattering terms cancel out. Instead of expressing all the terms analytically, we plot the ratio

(A1)$$R=\frac{{W\text{'}}_s({\overrightarrow{r}}_s,{\overrightarrow{r}}_d,\overrightarrow{r},{\mu \text{'}}_{a0},{D\text{'}}_0,\omega )}{{W"}_s({\overrightarrow{r}}_s,{\overrightarrow{r}}_d,\overrightarrow{r},{\mu "}_{a0},D_0,\omega )},$$

for the geometry of Fig1a and for the same absorption variation of the post ICG-breast. Figure A1 shows that the amplitude and phase of R remains very small for the most probable photon paths, (i.e. when the distance a is small). For higher a, R increases but similarly to the arguments that reduced Eq.16 to Eq.20, the contribution of terms for higher a is weighted less. Therefore for the small absorption changes considered here and for small variations of the diffusion coefficient δD′(r⃗)=δD″(r⃗) the scattering terms can be ignored when reconstructing the absorption perturbation only due to contrast agent injection. When introducing boundaries, the assumption of small a compared to the source-detector distance for the most probable photon paths works better for transmittance geometry.

 

Figure A1. (a) Amplitude and (b) Phase of R as a function of the average absorption coefficient of the post-ICG breast, plotted for the geometry of Fig 1c.

Download Full Size | PDF

Acknowledgements

The authors acknowledge useful discussions with Bruce Tromberg, Monica Holboke, Joe Culver and Deva Pattanayak. This work was supported in part by the National Institutes of Health grant no. RO1-CA75124 and in part by the National Institutes of Health grant no. RO1- CA60182.

References and links

1. B.W. Pogue, “Focus issue: Biomedical diffuse optical tomography - Introduction,” Opt Express , 4, 230–230 (1999). http://epubs.osa.org/opticsexpress/tocv4n8.htm [CrossRef]   [PubMed]  

2. F. Kelcz and G. Santyr, “Gadolinium-Enhanced Breast MRI,” Critical Reviews in Diagnostic Imaging 36, 287–338 (1995). [PubMed]  

3. F. Barkhof, J. Valk, O.R. Hommes, and P. Scheltens, “Meningeal Gd-DTPA enhancement in multiple-sclerosis,” Am. J. Neuroradiology 13, 397–400 (1992).

4. C. Tilcock, “Delivery of contrast agents for magnetic resonance imaging, computed tomography, nuclear medicine and ultrasound,” Adv Drug Deliver Rev 37, 33–51 (1999). [CrossRef]  

5. R.P. Kedar, D. Cosgrove, V.R. McCready, J.C. Bamber, and E.R. Carter, “Microbubble contrast agent for color Doppler US: Effect on breast masses,” Radiology 198, 679–686 (1996). [PubMed]  

6. M.L. Melany, E.G. Grant, and S. Farooki, et al. “Effect of US contrast agents on spectral velocities: In vitro evaluation”, Radiology 211, 427–431 (1999). [PubMed]  

7. S.E. Thompson, V. Raptopoulos, R.L. Sheiman, M.M.J McNicholas, and P. Prassopoulos, “Abdominal helical CT: Milk as a low-attenuation oral contrast agent,” Radiology 211, 870–875 (1999). [PubMed]  

8. R. Jain, S. Sawhney, P. Sahni, K. Taneja, and M. Berry, “CT portography by direct intrasplenic contrast injection: a new technique,” Abdominal Imaging 24, 272–277 (1999). [CrossRef]   [PubMed]  

9. L.W. Nunes, M.D. Schnall, S.G. Orel, M.G. Hochman, C.P. Langlotz, C.A. Reynolds, and M.H Torosian, “Correlation of lesion appearance and histologic findings for the nodes of a breast MR imaging interpretation model,” Radiographics 19, 79–92 (1999). [PubMed]  

10. S.G. Orel, M.D. Schnall, V.A. Livolsi, and R.H. Troupin, “Suspicious Breast-Lesions - MR-Imaging With Radiologic-Pathological Correlation,” Radiology 190, 485–493 (1994). [PubMed]  

11. B.W. Pogue, M. Testorf, T. McBride, U. Osterberg, and K. Paulsen, “Instrumentation and design of a frequency-domain diffuse optical tomography imager for breast cancer detection,” Opt. Express 1, 391–403 (1997). http://www.opticsexpress.org/oearchive/source/2827.htm [CrossRef]   [PubMed]  

12. V. Ntziachristos, X.H. Ma, and B. Chance, “Time-correlated single photon counting imager for simultaneous magnetic resonance and near-infrared mammography,” Rev. Sci. Instr. 69, 4221–4233 (1998). [CrossRef]  

13. V. Ntziachristos, A.G. Yodh, M.D. Schnall, and B. Chance, “Comparison between intrinsic and extrinsic contrast for malignancy detection using NIR mammography,” Proc. SPIE 3597, 565–570 (1999). [CrossRef]  

14. A. Ishimaru, Wave propagation and scattering in random media., (New York, Academic Press1978).

15. 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,” Opt. Soc. Am. A 10, 127–140 (1993). [CrossRef]  

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

17. S.R. Arridge, P. van der Zee, M. Cope, and D.T. Delpy, “Reconstruction methods for infra-red absorption imaging,” Proc. SPIE 1431, 204–215 (1991) [CrossRef]  

18. M.S. Patterson, B. Chance, and B.C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” J. Appl. Opt. 28, 2331–2336 (1989). [CrossRef]  

19. T.J. Farrell, M.S. Patterson, and B. Wilson, “A diffusion-theory model of spatially resolved, steady-state diffuse reflectance for the non-invasive determination of tissue optical properties in-vivo,” Med. Phys. 19, 879–888 (1992). [CrossRef]   [PubMed]  

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

21. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]  

22. 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]  

23. V. Ntziachristos, A. Hielscher, A.G. Yodh, and B. Chance, “Performance of perturbation tomography with highly heterogeneous media under the P1 approximation,” in Biomedical Optics : Advances in Optical Imaging, Photon Migration and Tissue Optics, OSA Technical Digest CLEO/Europe AMB3-1 : 211–213 (1999).

24. V. Ntziachristos, X.H. Ma, A.G. Yodh, and B. Chance, “Multichannel photon counting instrument for spatially resolved near infrared spectroscopy,” Rev. Sci. Instr. 70, 193–201 (1999). [CrossRef]