## Abstract

We introduce two filtering methods for near-field diffuse light diffraction tomography based on the angular spectrum representation. We then combine these filtering techniques with a new method to find the approximate depth of the image heterogeneities. Taken together these ideas improve the fidelity of our projection image reconstructions, provide an interesting three dimensional rendering of the reconstructed volume, and enable us to identify and classify image artifacts that need to be controlled better for tissue applications. The analysis is accomplished using data derived from numerical finite difference simulations with added noise.

© 1999 Optical Society of America

## 1. Introduction

Diffuse photon density waves (DPDWs) [1–7] and their time-domain analogs [8–10] can provide quantitative spectroscopic information about tissue structure, tissue chromophores, and tissue metabolism. Tomographic imaging of deep tissues with diffusive waves offers the possibility for functional imaging based on these spectroscopic contrasts. Thus image reconstruction has been explored intensely with varying degrees of success in tissue phantoms [11,12], the brain [13–18], and the breast [19–24]. One theoretical approach that has received significant recent interest is the method of diffraction tomography [25–28] extended to the near-field [29–43]. This scheme is particularly well suited for planar source-detection geometries that arise for example, in breast imaging. The methodology has attracted attention [29–43] because the diffusive wave modes derived using the angular-spectrum representation provide a very convenient framework for imaging in parallel-plate geometries, and for analytic or semi-analytic investigation of resolution and signal-to-noise issues. The method also offers a very fast and direct way to obtain two-dimensional projection images of the turbid media; these images are essentially optical mammograms[21].

In this contribution we study a set of theoretical issues associated with image reconstruction using near-field diffraction tomography. In order to obtain quality images, for example, one apply spectral filters [41–46] to the data at several levels of the image processing (i.e. filters with respect to spatial frequencies in the reconstructions). Well defined rules do not exist for choosing these filter functions; in fact different filter functions introduce uncontrolled systematic errors into the images [41] . Additionally, the images are dramatically improved when the experimenter has knowledge of the approximate depth of the optical heterogeneity. Object depth determination is generally difficult unless one has other means to obtain this information such as multiple optical views of the same sample or image segmentation from other types of diagnostics.

We investigate both critical issues in this paper, and demonstrate algorithms for optimizing projection image fidelity. Data for these studies is derived from numerical simulations with added noise. Two kinds of mathematical filters are introduced; (1) a phase-only filter which does not have any free parameters and allows accurate localization, and (2) a depth and noise dependent low-pass filter with the cut-off frequency as a free parameter. Then the two filtering methods are combined with a technique to find the approximate depth of the heterogeneities. We obtain two dimensional projection images and we demonstrate a three dimensional rendering of these image projections which appears promising for clinical applications. Our work clarifies the mentioned issues in a quantitative way; it makes apparent the limitations of the technique, identifies artifacts and directions for improvement.

## 2. Theory

#### 2.1 Angular spectrum representation of diffuse photon density waves

The angular-spectrum representation provides a set of modes well suited for the
description of the propagation and scattering of diffuse photon density waves in
parallel-plate geometries (fig.(1)). The representation and its applications have
been described by numerous investigators in optics [25–28,47–48] and in the photon migration community [29–43]. For clarity we review the salient
features of this approach below. We assume infinite space boundary conditions in
this discussion; Green’s functions for semi-infinite and slab
geometries have been reported in the literature [33–35] and do not qualitatively change the
conclusions presented herein. The starting point for our treatment is the formal
expression relating the scattered diffuse light field, Φ_{sc} , to a volume integral over heterogeneities within the entire turbid
medium [33]:

Here *T*(**r**′) is
“tumor” function at position **r**′
which, in the Born Approximation, depends on the background diffuse photon
density wave, Φ_{0}(**r**), as well as the
absorption
(*δμ _{a}*(

**r**′)) and scattering (

*δμ*′

_{s}(

**r**′)) deviations from the medium’s background absorption (

*μ*

_{a0}) and scattering coefficients (

*μ*′

_{s0}) respectively.

*T*(

**r**′) also depends on the speed of light in the medium,

*υ*, the background light diffusion constant

*D*

_{0}, the background diffuse photon density wavevector ${k}_{0}=\sqrt{\left(-\upsilon {\mu}_{a0}\text{}+i\omega \right)/{D}_{0})}$ , and the modulation frequency,

*ω*. It can be conveniently separated into absorptive and scattering parts [49], i.e.,

Here we drop the gradient terms due to the scattering variations for simplicity.
This approximation is described in detail in [49]. *G*
_{0}(**r**,
**r**′) is the usual Green’s function for the
medium of the form:

Typical experiments that employ the angular-spectrum representation for analysis
use a “parallel-plate” geometry (e.g. fig.(1)) where the detected DPDW,
Φ(**r**) = Φ_{sc}(**r**)+ Φ_{0}(**r**) , is
measured at discrete points within a square grid on the planar output (detector
plane, (*x*,*y*,*z _{d}*))
surface. We take our source (at the origin, (0,0,0)) to be a point emitter on
the planar input surface (source plane,
(

*x*,

_{s}*y*,0)). In this geometry the Green’s function is conveniently written in the following form:

_{s}where the z-axis is defined in the direction normal to the plane surfaces, and
(*p*, *q*) are the 2D spatial (k-space)
frequencies conjugate to the *x*- and
*y*-coordinates. The angular spectrum (Weyl expansion)
representation of the Green’s function is [31,50]

where $m=\sqrt{{k}_{0}^{2}-\left({p}^{2}+{q}^{2}\right)}$ is the complex propagator wavenumber and the imaginary part of
*m* is positive, i.e. ℑ(*m*)
> 0. With these definitions, the 2D Fourier Transform of the
scattered field in the detection plane is simply [33]

where hat (̂) indicates quantities in k-space. This equation is the fundamental equation for diffraction tomography.

#### 2.2 Image reconstruction

Equation (7) relates the tumor function in each slice to the 2D
Fourier Transform of the scattered DPDW in the detection plane. Rather than
inverting the Laplace Transform directly, we discretize the integral into a sum
over slices of thickness ∆*z*, i.e

∆*z* is the step size in the z-direction and
*N* =
*z _{d}*/∆

*z*is the total number of slices. This is the key equation for our image reconstruction using the angular-spectrum approach. There are a number of possible ways to solve this problem in 3D [36,40], but from this point onwards we focus on projection images and a 3D variant thereof.

For two dimensional projection images in the x-y planes along the depth direction (z), we assume that the inhomogeneities are isolated and thin, and we drop the sum. That is, we take the tumor function to be zero everywhere except at the phantom (object) depth. Thus a measurement, Fourier transform, and a simple division in the k-space yields the tumor function, i.e

where the subscript “obj” indicates the phantom coordinates.

The inverse Fourier transform of this quantity enables us to solve for absorption and scattering variations. Using

we can solve eq. (2) and (3) to get
*δμ _{a}*(

**r**) and

*δμ*′

_{s}(

**r**). Note that we need to know the background field at the object depth to accurately to obtain optical properties. Fig.(2)is a flowchart illustrating the algorithm described above.

The relevant clinical situation is shown in fig.(1). The breast is embedded in a slab-like box filled with matching material such as Intralipid. The optical properties of Intralipid are chosen to be close to the average optical properties of the breast. Illumination comes from one plane, and the detectors are scanned in a grid on the opposite plane. The amplitude and phase of the DPSW’s with and without the breast are measured and used to obtain the scattered DPDW. In our experiments [21,33,41] we also use multiple optical wavelengths to obtain spectroscopic information[41]. The background field at different depths is not measured and must be estimated using numerical models or a backpropagation technique [25,51]. This is the main difficulty for obtaining 3D information in clinical situations and we are investigating various methods to overcome this problem [33,34].

#### 2.3 Noise amplification at high spatial frequencies

In eq.(9) we divide two waves with complex wavenumbers
(*m* and *k*
_{0}). The scattered field
generally suffers from noise; the origin of this noise can be experimental,
numerical, or due to finite sampling. The noise in the high spatial frequencies
is preferentially amplified by the deconvolution because the Green’s
function in the denominator of equation (9) decays to zero faster for higher values of
*p* and *q*. The amplification of the high
frequency noise and the associated instabilities are common problems associated
with Fourier convolution approaches to image reconstruction [25–28,44–46]. Fig.(3) demonstrates this problem.

One [37,38,41] way of dealing with this problem is to apply a low-pass spatial filter in order to suppress the high frequency components of the tumor function. The quantitative effects of this cut are not known however [41,45,46]. Furthermore, the spatial resolution of the images depend on the highest k-space frequencies available and therefore filtering modifies image resolution (see eq.(12))[31,41,43,44].

## 3. Filtering and normalization for improving depth information / image quality

We have investigated filtering and normalization techniques which provide different types of information and lead to the improvement of image quality within the context of projection imaging. In particular we have found a robust method which enables the experimenter to estimate the approximate depth of the heterogeneities below the detection plane; the method does not require more data, yet offers a prescription for extending the 2D projection approach to quasi-3D imaging and localization. In this section we briefly outline three techniques to optimize image reconstruction. In the remainder of the paper we demonstrate their utility. For simplicity we focus entirely on absorption heterogeneities.

#### 3.1 Phase only Green’s function (real-m filter)

We have found a very crude filtering technique that gives fairly accurate
positional information in two- and even three-dimensional images. In essence the
method is based on the hypothesis that most of the positional information is
contained in the phase of our signals in k-space. In this procedure we modify
the Green’s function, *G*̂_{0}
(*p*, *q*, *Z _{d}*,

*z*′); in particular we put

*G*̂

_{0}(

*p*,

*q*,

*Z*,

_{d}*z*′) ≈ ie

^{iℜ (m)|zd-z′|}. We set its amplitude to unity, so only phase information in the Green’s function is used for the reconstruction. We also apply a Blackman Window on

*T*̂ (see [43] for its properties) which is a standard windowing function used in the Fourier domain to further stabilise the image. This last step is used for better image quality and is not an essential part of the algorithm. Calculation of the projection image proceeds along similar lines except for these steps (see fig. (2)). Although quantitative information about optical properties is largely lost, the real-m filter has no free parameters (except from the optional windowing function) and provides superior information about the xyz central coordinates of heterogeneities.

#### 3.2 Depth estimation (S_{j})

We have empirically found [52] that the following set of operations applied to the
projection images obtained with the real-m filter enable us to accurately
estimate the heterogeneity depth below the detection surface. Briefly, we first
obtain the tumor function for each slice (centered on some
*z _{j}*) using the real-m filter. We then derive an image
of

*δμ*(

_{a}*x*,

*y*,

*z*) for each slice, determine the center position of the object we are trying to characterize, and record the absorption variation at the object central position. Finally we compute the quantity

_{j}where the sum in the denominator is over all pixels in the image slice.
*δμ*̅_{a} is the mean *δμ*
_{axy} in the *j ^{th}* plane.

*S*is in essence a measure of the contrast of the object in the

_{j}*j*plane. We empirically found that the value of

^{th}*z*for which

_{j}*S*is a maximum closely approximates the actual depth of the object (

_{j}*z*). This procedure enables us to select a slice for the projection image and, thus provides a scheme for extending images to three dimensions (see section

_{obj}*5.1 .1*).

#### 3.3 High frequency filtering based on depth and experimental noise floor (G-filtering)

After we have obtained this rough picture of the object in the turbid medium, we
repeat the entire procedure using the true Green’s function. From our
discussion in section *2.3* however, it is apparent that we need
to devise a scheme to select the spatial frequencies for reconstruction. We
carry out this procedure in a straightforward way. Note that it is critical for
this procedure that the scattered waves are normalized by the source strength.
For our convenience we used a source strength of one photon per unit volume per
unit time. We first plot the raw data to determine the experimental noise level.
For example by plotting Φ̂_{sc} vs. *p* and *q* (e.g. fig.(3)) we will find at sufficiently high
*p* and *q*, that Φ_{sc} hits a “white” noise floor; the ratio of the signal
amplitude at *p* = *q* = 0 to the signal amplitude
at these high frequencies provides a measure of the experimental
signal-to-noise. Next we insert the estimated object depth
*z*′ = *z _{obj}* into the
expression for

*G*̂

_{0}(

*p*,

*q*,

*z*,

_{d}*z*′) , and we determine at what value of m,

*G*̂

_{0}(

*p*,

*q*,

*z*,

_{d}*z*) drops below a threshold based on our experimental S/N. In practice we set this threshold to be about five times the white noise floor. We set

_{obj}*T*̂ = 0 for all

*m*greater than this threshold (see fig (2) for the occurance of this step in the general algorithm.). The result is a depth dependent pillbox filter applied to the tumor function. The rest of the reconstruction is the same. The combination technique offers the possibility of improvement in the accuracy of optical properties and better images.

At this point it may also be desirable to recompute
*S _{j}* using reconstructions based on the G-filter
technique. Usually the maximum

*S*determined with the G-filter is within ≈1 cm of the maximum

_{j}*S*determined using the real-m filter (see section

_{j}*3.1*). In applying this latter procedure the high frequency cut-off depends on depth (

*z*).

_{j}## 4. Data generation

Our “data” was derived from numerical simulations. In
particular the finite difference method with partial current boundary conditions was
used to solve the diffusion equation for arbitrary interior inhomogeneities. The
program is capable of simulating the forward problem in a box geometry. We used a
point source with a strength of one photon per unit volume per unit time; this
choice of source term eliminated the need to normalize the scattered wave when
computing the noise floor and comparing the threshold to the Green’s
function. The source was modulated at 140Mhz. For the purposes of this paper we
simulate data in a 20×20×20cm box and use the central
9×9×5cm region for “experiments”. The
region is far from all boundaries and is a very good approximation to the infinite
medium. The region created from the simulation is a 65 by 65 by 36 grid in x, y and
z with the coordinates respectively, running from -4.5 to 4.5cm in x & y
and 0 to 5 cm in z directions. The step size is 0.14cm corresponding to a voxel size
of (0.14cm)^{3}. The geometry chosen closely mimics our clinical set-up [34,41]. We use a high number of pixels in the numerical simulation
to avoid truncation errors. Each point of the detection grid is centered in the
squares of the same 65×65 region in the detector plane at
*z _{j}* =

*z*= 5cm.

_{d}Optical properties can be assigned to each pixel independently so that we are able to simulate heterogeneities of arbitrary shape. We use thin rectangular slices (i.e thin in z-direction) for objects. Thin slices insure that our projection approximation is reasonable. Figs. (4) and (8) provide 3D renderings of the input phantoms.

Noise is added to the data using a random number generator with a normal distribution of variance one and mean zero. The approach is tuned to provide gaussian noise with variance 0.5% for the amplitude and variance of 0.05° for the phase. These probably overestimate our experimental noise. Noise is added to heterogeneous and homogeneous measurements. Figs.(4) and (12) show the amplitude of the noiseless and noisy scattered fields obtained in Born approximation.

## 5. Results

#### 5.1 Tissue phantoms with single slice heterogeneity

The simulated heterogeneity is shown in fig.(4). It has dimensions 1.4 × 1.4
× 0.28 cm in x,y and z directions respectively. The optical
properties are *μ _{a}* =0.1 cm

^{-1}and

*μ*′

_{s}= 8 cm

^{-1}. The background optical properties are

*μ*= 0.02 cm

_{ao}^{-1}and

*μ*= 8 cm

_{so}^{-1}. The phantom is centered at a depth of 2.65cm from the source plane. It’s transverse location centers on y=x= 1.97cm. The amplitude of the “measured” scattered field at the detector plane is shown in its noiseless and noisy versions in fig.(4) for comparison.

### 5.1 .1 Real-m imaging

We now employ the real-m filter scheme described in section
*3.1* on the single object data shown above. Figs.(5a) and (5b) show two dimensional projection images at depths
z=2.42cm and z=3.24cm respectively. We see that the object is clearly
apparent. In fig. (5c), we plot *S _{j}* vs

*z*; we see that the contrast parameter is peaked at

_{j}*z*= 2.71cm, and we use this value to make a “true” real-m filter reconstruction of the data which is shown in fig.(5d). This result is representative of all of the single-object phantoms that we tested.

_{j}Although most of the information about optical properties is lost, the 3D localization of single objects is quite good. The slices with the greatest noise are either near the detection or source planes. The issue of multiple objects at different depths will be addressed later in this paper.

### 5.1 .2 G-filter imaging

Next we use the combination scheme described in section *3.3*
to reconstruct the single object data. Our experimental noise threshold was
3×10^{-3} for all cases. For
*z _{j}* = 2.71cm, the G-filter was set for p=q=
1.01/cm at maximum, and the resulting image in this plane is shown in fig.(6a) . Notice that the image is a little cleaner
than the real-m image with most of the artifacts from the characteristic
ringing near the corners of the reconstructed volume. The optical properties
(although still not accurate for such a deep object) are more reasonable.

We next computed *S _{j}* vs

*z*using the G-filter (see fig.(6b)). Notice that the peak value of

_{j}*S*has shifted from its true value. The G-filter image (at p=q=2.36/cm) based on this new

_{j}*z*is shown in fig.(6c); we see that the image is sharper and has better optical properties. This is very important for multiple wavelength images as we observe that the ratio of reconstructed optical properties is fairly well approximated by this technique. However, its actual location is systematically shifted away from its true value. At present we do not fully understand the origin of this systematic error in a deep way. We find that this shift is related to the amount of power cut out by the filter. By changing the filter spatial frequency cut-off to a lower frequency (increasing the threshold) we find that the shift is increased but the optical properties at the shifted location however approach their correct values.

_{j}Finally we note that since the maximum spatial frequencies used in the
G-filter approach depend on object depth, we can expect our resolution to
decrease with increasing object-detector separation. Following Pattanayak [31] our approximate experimental resolution,
*L*, depends on the maximum allowable *p* and
*q*, i.e

In fig.(7) we show the change in resolution with depth assuming our experimental S/N threshold. While the exact numerical values depend on approximations in [31], the important qualitative effect to note is that the resolution improves dramatically for objects near the detection plane.

#### 5.2 Tissue phantoms with two slice heterogeneities / three dimensional renderings

We now simulate two objects with same dimensions and optical properties but in a more complicated geometry where in principle the scattered waves from one heterogeneity could effect the other. One object is centered at 3.07cm deep with its transverse center at x=y=1.97cm and the other object is 4.35cm deep with its transverse center at x=y=-1.97cm as shown in fig.(8). The amplitude of the “measured” scattered field at the detector plane is shown in its noiseless and noisy version in fig.(8) for comparison.

We performed essentially the same set of reconstruction procedures as described
in sections **3**. and *5.1* . In fig.(9a) we show that *S _{j}* has
two maxima in different planes depending on the heterogeneity under
consideration using real-m filter. Note that the real-m images were more noisy
by comparison to the single object phantom. We again find that the optical
properties are most accurate where

*S*has its shifted maxima with the G-filter. These images are shown in Fig (9b,c). Absolute values are closer to the real values for shallow objects and the positional shift also depends on the depth.

_{j}Finally, by calculating *S _{j}* at all the voxels we can
generate three dimensional images of the entire reconstructed volume. That is,
for each of the 36 projections we calculate the value of

*S*for all transverse pixels (65×65= 4225 values in each projection) and plot all

_{j}*S*’s in three dimensions. In fig.(10a) and (10b) we show an isosurface rendering at value

_{j}*S*= 0.042 at two viewing angles. Two dimensional projections in plane y=2cm and plane y=-2cm are also shown in fig. (10c) and fig.(10d) to provide a better sense of contrast. Comparing this to fig.(8) we see that the shallow object is very well reconstructed both in terms of shape and location. The deeper object has lower resolution and is shifted from its original location. In the corners we see the characteristic ringing effects from Fourier domain reconstructions. These latter artifacts, however, are usually fairly easy to identify. We find that, even when the objects are close to the image edges we are capable of separating objects from ringing artifacts. We are investigating presently whether it is possible to extract any quantitative or qualitative optical property information from these images.

_{j}## 6. Conclusion

We have studied projection images of turbid media based on the angular-spectrum
representation. We demonstrated two filtering schemes. One of them, the real-m
filter, provided object location information with no free parameters; the other
method, the G-filter, introduced a systematic approach for eliminating high spatial
frequencies from the reconstructions. Image localization in 3D was demonstrated
using both filters by maximizing the object signal-to-noise
(*S _{j}*) on a slice-by-slice basis. Despite these successes,
some image artifacts remain and must be addressed in future work. These include:
noise near the detection plane and image boundaries, systematic shifts of the object
location depending on filtering schemes and noise floors, and optical property
accuracy.

## Acknowledgements

We gratefully acknowledge useful conversations with Hanli Liu, D. A. Boas and C. L. Matson. A.G. Y. and D. N. P. acknowledge partial support from US Army RO DAMD17-97-1-7272 and NIH 1-RO1-CA75124-01 grants respectively.

## References and links

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

**2. **B. Chance, in *Photon Migration in Tissues* (Plenum Press, 1989).

**3. **B.J Trombeg, L. O. Svaasand, T. Tsay, and R.C.M. Haskell, “Properties of Photon Density Waves in Multiple-Scattering Media,” Appl. Opt. **32**, 607 (1993). [CrossRef]

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

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

**6. **Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R. L. Barbour, “Frequency-domain optical imaging of absorption and scattering distributions using a born iterative mehod,” J. Opt. Soc. Am. A **14**, 325–342 (1997). [CrossRef]

**7. **M. O’Leary, “Imaging with Diffuse Photon Density Waves,” in PhD Thesis (Dept. Physics & Astronomy, U. of Pennsylvania, May 1996) .

**8. **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 (1989). [CrossRef] [PubMed]

**9. **S. R. Arridge, “Forward and inverse problems in time-resolved infrared imaging,” in *Medical Optical Tomography: Functional Imaging and Monitoring*, ed. G. Muller, B. Chance, Rl. Alfano, S. Arridge, J. Beuthan, E. Gratton, M Kaschke, B. Masters, S. Svanberg, and P. van der Zee, Proc SPIE **IS11**, 35–64 (1993).

**10. **D. A. Benaron, D. C. Ho, S. Spilman, J. P. Van Houten, and D. K. Stevenson, “Tomographic time-of-flight optical imaging device,” Adv. Exp. Med. Biol. **361**, 609–617 (1994). [CrossRef]

**11. **E. Gratton and J. B. Fishkin, “Optical spectroscopy of tissue-like phantoms using photon density waves,” Comments on Cell. and Mol. Biophys. **8(6)**, 309–359 (1995).

**12. **J. B. Fishkin, O. Coquoz, E. R. Anderson, M. Brenner, and B. J. Tromberg, “Frequency-domain photon migration measurements of normal and malignant tissue optical properties in a human subject,” Appl. Opt. **36**, 10 (1997). [CrossRef] [PubMed]

**13. **W. Bank and B. Chance, “An Oxidative Effect in metabolic myopathies - diagnosis by noninvasive tissue oximetry,” Ann. Neurol. **36**, 830 (1994). [CrossRef] [PubMed]

**14. **Y. Hoshi and M. Tamura, “Near-Infrared Optical Detection of Sequential Brain Activation in The Prefrontal cortex during mental tasks,” Neuroimage. **5**, 292 (1997). [CrossRef] [PubMed]

**15. **A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain functions,” Trends. Neurosci. **20**, 435 (1997) . [CrossRef] [PubMed]

**16. **B. Chance, Q. M. Luo, S. Nioka, D. C. Alsop, and J.A. Detre, “Optical investigations of physiology: a study of intrinsic and extrinsic biomedical contrast,” Phil. Trans. Roy. Soc. London B. **352**, 707 (1997). [CrossRef] [PubMed]

**17. **BW Pogue and KD Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of apriori magnetic resonance imaging structural information,” Opt. Lett. **23**,1716–1718 (1998). [CrossRef]

**18. **R. M. Danen, Y. Wang, X. D. Li, W. S. Thayer, and A.G. Yodh, “Regional imager for low-resolution functional imaging of the brain with diffusing near-infrared light,” Photochem. Photobiol. **67**, 33–40 (1998). [CrossRef] [PubMed]

**19. **J.H. Hoogenraad, M.B.van der Mark, S.B. Colak, G.W. ’t Hooft, and E.S.van der Linden, “First Results from the Philips Optical Mammoscope,” Proc.SPIE / BiOS-97 (SanRemo, 1997) .

**20. **S.K. Gayen, M.E. Zevallos, B. B. Das, and R. R. Alfano and “Time-sliced transillumination imaging of normal and cancerous breast tissues,” in *Trends in Opt. And Photonics*, ed. J. G. Fujimoto and M. S. Patterson.

**21. **X.D. Li, J. Culver, D. N. Pattanayak, A. G. Yodh, and B. Chance, “Photon Density Wave Imaging With K-Space Spectrum Analysis: clinical studies - background substraction and boundary effects,” Technical Digest Series - CLEO ’98 , **6**, 88–89 (1998).

**22. **S. Fantini, S. A. Walker, M. A. Franceschini, K. T. Moesta, P. M. Schlag, M. Kaschke, and E. Gratton. “Assessment of the size, position, and optical properties of breast tumors in vivo by non-invasive optical methods” Appl. Opt. , **37**, 1982–1989 (1998). [CrossRef]

**23. **M. A. Franceschini, K. T. Moesta, S. Fantini, G. Gaida, E. Gratton, H. Jess, W. W. Mantulin, M. Seeber, P. M. Schlag, and M. Kaschke. “Frequency-domain instrumentation techniques enhance optical mammography: Initial clinical results” Proc. Natl. Ac ad. Sci. USA , **94**, 6468–6473 (1997). [CrossRef]

**24. **S. Fantini, M. A. Franceschini, G. Gaida, E. Gratton, H. Jess, and W.W. Mantulin, “Frequency-domain optical mammography: edge effect corrections,” Med. Phys. **23**, 149 (1996). [CrossRef] [PubMed]

**25. **E. Wolf, “Principles and Development of Diffraction Tomography” in *Trends in Optics*, ed. A. Consortini (Academic Press, San Diego, 1996). [CrossRef]

**26. **A.J. Devaney, “Diffraction Tomography,” Inv. Meth. In Electromagnetic Imaging, 1107–1135 .

**27. **E. Wolf, “Inverse Diffraction and a New Reciprocity Theorem,” J. Opt. Soc. Am. **58**, 1568 (1968).

**28. **E. Wolf , “Three Dimensional Structure Determination of Semi-Transparent Objects From Holographic Data,” Opt. Commun. **1** , 153–156 (1969). [CrossRef]

**29. **BQ Chen, JJ Stamnes, and K Stamnes, “Reconstruction algorithm for diffraction tomography of diffuse photon density waves in a random medium”. Pure Appl Opt ,**7**, 1161–1180 (1998). [CrossRef]

**30. **DL Lasocki, CL Matson, and PJ Collins, “Analysis of forward scattering of diffuse photon-density waves in turbid media: a diffraction tomography approach to an analytic solution,” Opt. Lett. **23**,558–560 (1998). [CrossRef]

**31. **D. N. Pattanayak, “Resolution of Optical Images Formed by Diffusion Waves in Highly Scattering Media,” GE Tech. Info. Series **91CRJ3241** (1991).

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

**33. **X.D. Li ,in PhD Thesis (Dept. Physics & Astronomy, U. of Pennsylvania, May 1998).

**34. **X.D. Li, D. N. Pattanayak, J. P. Culver, T. Durduran, and A. G. Yodh, “Near-Field Diffraction Tomography with Diffuse Photon Density Waves,” to be published (1998).

**35. **X. Cheng and D. Boas, “Diffuse Optical Reflection Tomography Using Continous Wave Illumination,” Opt. Express **3**, 118–123 (1998), http://epubs.osa.org/oearchive/source/5663.htm. [CrossRef] [PubMed]

**36. **J. C. Schotland, “Near-field Inverse Scattering: Microscopy to Tomography,” SPIE3597 (1999).

**37. **C. L. Matson, N. Clark, L. McMackin, and J. S. Fender, “Three-dimensional Tumor Localization in Thick Tissue with The Use of Diffuse Photon-Density Waves,” Appl. Opt. **36**, 214–219 (1997). [CrossRef] [PubMed]

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

**39. **S. J. Norton and T. Vo-Dinh, “Diffraction Tomographic Imaging With Photon Density Waves: an Explicit Solution,” J. Opt. Soc. Am. A **15**, 2670–2677 (1998). [CrossRef]

**40. **J. C. Schotland, “Continous Wave Diffusion Imaging,” J. Opt. Soc. Am. A **14**, 275–279 (1997). [CrossRef]

**41. **T. Durduran, J. Culver, L. Zubkov, M. Holboke, R. Choe, X. D. Li, B. Chance, D. N. Pattanayak, and A. G. Yodh, “Diffraction Tomography In Diffuse Optical Imaging; Filters & Noise,” SPIE3597 (1999).

**42. **J. Ripoll and M. Nieto-Vesperinas, “Reflection and Transmission Coefficients of Diffuse Photon Density Waves,” in press.

**43. **J. Ripoll and M. Nieto-Vesperinas, “Spatial Resolution of Diffuse Photon Density Waves,”to be published in J. Opt. Soc. Am.A (1999).

**44. **C. L. Matson, “Resolution, Linear Filtering , and Positivity,” J. Opt. Soc. Am. A **15**, 33–41 (1998). [CrossRef]

**45. **F. J. Harris, “On The Use of Windows For Harmonic Analysis with the Discrete Fourier Transform,” Proc. Of IEEE **66**, 51–83 (1978). [CrossRef]

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

**47. **A.J. Devaney, “Linearised Inverse Scattering in Attenuating Media,” Inv. Probs. **3**, 389–397 (1987). [CrossRef]

**48. **A. J. Devaney, “Reconstructive Tomography With Diffracting Wavefields,” Inv. Probl. **2**, 161–183 (1986). [CrossRef]

**49. ** Essentially we assume that the scattering contrast (*δμ*′_{s}) is slowly varying. For a detailed description we refer to [33] and [34].

**50. **A.J. Banos, in *Dipole Radiation In the Presence of a Conducting Half-Space* (Pergamon Press, New York, 1966).

**51. **J. W. Goodman, in *Introduction To Fourier Optics*, (McGraw-Hill, San Fransisco , 1968).

**52. ** We are aware of a similar normalization scheme by Hanli Liu and her collaborators ( private communications SPIE Jan 1999).