Algorithms for 3D localization and imaging using near-field diffraction tomography with diffuse light

We introduce two ltering methods for neareld di use light di raction tomography based on the angular spectrum representation. We then combine these ltering techniques with a new method to nd the approximate depth of the image heterogeneities. Taken together these ideas improve the delity 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 nite difference simulations with added noise. c 1999 Optical Society of America OCIS codes: (170.3010) Image reconstruction techniques;(170.5280) photon migration; (170.3830) mammography; (170.5270) Photon density waves; (170.3660) Light propagation in tissues References and links 1. A. G. Yodh and B. Chance, \Spectroscopy and imaging with di using light," Physics Today 48, 34-40 (March 1995). 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). 4. M. A. O'Leary, D. A . Boas, B. Chance and A.G. Yodh, \Experimental Images of heterogeneous turbid media by frequency-domain di using-photon tomography," Opt. Lett. 20, 426-428 (1995). 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). #9187 $15.00 US Received March 01, 1999; Revised March 29, 1999 (C) 1999 OSA 12 April 1999 / Vol. 4, No. 8 / OPTICS EXPRESS 247 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). 7. M. O'Leary, \Imaging with Di use 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 Refectance And Transmittance for the Non-Invasive Measurement of Tissue Optical Properties," Appl. Opt. 28, 2331 (1989). 9. S. R. Arridge, \Forward and inverse problems in time-resolved infrared imaging," in Medical Optical Tomography: Functional Imaging and Monitoring, ed. G. Muller and 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-fight optical imaging device," Adv. Exp. Med. Biol. 361, 609-617 (1994). 11. Gratton, E. 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). 13. W. Bank and B. Chance, \An Oxidative E ect in metabolic myopathies diagnosis by noninvasive tissue oximetry," Ann. Neurol. 36, 830 (1994). 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). 15. A. Villringer and B. Chance, \Non-invasive optical spectroscopy and imaging of human brain functions," Trends. Neurosci. 20, 435 (1997) . 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). 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). 18. R. M. Danen, Y. Wang, X. D. Li, W. S. Thayer and A.G.Yodh, \Regional imager for lowresolution functional imaging of the brain with di using near-infrared light," Photochem. Photobiol. 67, 33-40 (1998). 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, 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 e ects," Technical Digest Series CLEO '98, 6, 88-89 (1998). 22. Fantini, S., 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). 23. Franceschini, M. A., 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). 24. S. Fantini, M. A. Franceschini, G. Gaida, E. Gratton, H. Jess and W.W. Mantulin, \Frequencydomain optical mammography: edge e ect corrections," Med. Phys. 23, 149 (1996). 25. E. Wolf, \Principles and Development of Di raction Tomography" in Trends in Optics, ed. A. Consortini (Academic Press, San Diego, 1996). 26. A.J.Devaney, \Di raction Tomography," Inv. Meth. In Electromagnetic Imaging, 1107-1135 . 27. E. Wolf, \Inverse Di raction 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). #9187 $15.00 US Received March 01, 1999; Revised March 29, 1999 (C) 1999 OSA 12 April 1999 / Vol. 4, No. 8 / OPTICS EXPRESS 248 29. BQ Chen, JJ Stamnes, and K Stamnes, \Reconstruction algorithm for di raction tomography of di use photon density waves in a random medium". Pure Appl Opt ,7, 1161-1180 (1998). 30. DL Lasocki, CL Matson and PJ Collins, \Analysis of forward scattering of di use photon-density waves in turbid media: a di raction tomography approach to an analytic solution," Opt. Lett. 23,558-560 (1998). 31. D. N. Pattanayak, \Resolution of Optical Images Formed by Di usion Waves in Highly Scattering Media," GE Tech. Info. Series 91CRD241 (1991). 32. X.D.Li, T. Durduran, A.G. Yodh, B. Chance and D.N. Pattanayak, \Di raction Tomography for Biomedical Imaging With Di use Photon Density Waves," Opt. Lett. 22, 573-575 (1997). 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 Di raction Tomography with Di use Photon Density Waves," to be published (1998). 35. X. Cheng and D. Boas, \Di use Optical Refection Tomography Using Continous Wave Illumination," Opt. Express 3, 118-123 (1998), http://epubs.osa.org/oearchive/source/5663.htm. 36. J. C. Schotland, \Neareld Inverse Scattering: Microscopy to Tomography," SPIE 3597 (1999). 37. C. L. Matson, N. Clark, L. McMackin and J. S. Fender, \Three-dimensional Tumor Localization in Thick Tissue with The Use of Di use Photon-Density Waves," Appl. Opt. 36, 214-219 (1997). 38. C. L. Matson, \A Di raction Tomographic Model Of The Forward Problem Using Di use Photon Density Waves," Opt. Express 1, 6-11 (1997), http://epubs.osa.org/oearchive/source/1884.htm. 39. S. J. Norton and T. Vo-Dinh, \Di raction Tomographic Imaging With Photon Density Waves: an Explicit Solution," J. Opt. Soc. Am. A 15, 2670-2677 (1998). 40. J. C. Schotland, \Continous Wave Di usion Imaging," J. Opt. Soc. Am. A 14, 275-279 (1997). 41. T. Durduran, J. Culver, L. Zubkov, M. Holboke, R. Choe, X. D. Li, B. Chance, D. N. Pattanayak, A. G. Yodh, \Di raction Tomography In Di use Optical Imaging; Filters & Noise," SPIE 3597 (1999). 42. J. Ripoll and M. Nieto-Vesperinas, \Refection and Transmission Coe cients of Di use Photon Density Waves,"in press. 43. J. Ripoll and M. Nieto-Vesperinas, \Spatial Resolution of Di use 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). 45. F. J. Harris, \On The Use of Windows For Harmonic Analysis with the Discrete Fourier Transform," Proc. Of IEEE 66, 51-83 (1978). 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). 48. A. J. Devaney, \Reconstructive Tomography With Di racting Wave elds," Inv. Probl. 2, 161-183 (1986). 49. Essentially we assume that the scattering contrast (δμ0 ) is slowly varying. For a detailed des scription 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). #9187 $15.00 US Received March 01, 1999; Revised March 29, 1999 (C) 1999 OSA 12 April 1999 / Vol. 4, No. 8 / OPTICS EXPRESS 249

Abstract: We introduce two ltering methods for near-eld diuse light di raction tomography based on the angular spectrum representation. We then combine these ltering techniques with a new method to nd the approximate depth of the image heterogeneities. Taken together these ideas improve the delity 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 nite difference simulations with added noise.
In this contribution we study a set of theoretical issues associated with image reconstruction using near-eld di raction tomography. In order to obtain quality images, for example, one apply spectral lters [41][42][43][44][45][46] to the data at several levels of the image processing (i.e. lters with respect to spatial frequencies in the reconstructions). Well dened rules do not exist for choosing these lter functions; in fact dierent lter 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 dicult 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 delity. Data for these studies is derived from numerical simulations with added noise. Two kinds of mathematical lters are introduced; (1) a phase-only lter which does not have any free parameters and allows accurate localization, and (2) a depth and noise dependent low-pass lter with the cut-o frequency as a free parameter. Then the two ltering methods are combined with a technique to nd 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 claries the mentioned issues in a quantitative way; it makes apparent the limitations of the technique, identies artifacts and directions for improvement.

Angular spectrum representation of di use photon density waves
The angular-spectrum representation provides a set of modes well suited for the description of the propagation and scattering of di use photon density waves in parallel-plate geometries (g.(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][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. For clarity we review the salient features of this approach below. We assume innite space boundary conditions in this discussion; Green's functions for semi-innite and slab geometries have been reported in the literature [33][34][35] and do not qualitatively change the conclusions presented herein. The starting point for our treatment is the formal expression relating the scattered diuse light eld, sc , to a volume integral over heterogeneities within the entire turbid medium [33]: The generic near-eld di raction tomography experiment. The detector is scanned in a 2D grid on the surface of the plane parallel to and displaced from the plane containing the source. The breast is embedded in the box along with Intralipid in order to match its average optical properties.
Here T (r 0 ) is \tumor" function at position r 0 which, in the Born Approximation, depends on the background diuse photon density wave, 0 (r), as well as the absorption (δµ a (r 0 )) and scattering (δµ 0 (r 0 ) ) deviations from the medium's background absorption s 0 (µ a0 ) and scattering coecients (µ ) respectively. T (r 0 ) also depends on the speed s0 of light in the medium, v, the background light diusion constant D 0 , the background p diuse photon density wavevector k 0 = (−vµ a0 + iω)/D 0 ) , and the modulation frequency, ω. It can be conveniently separated into absorptive and scattering parts [49], 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 0 ) is the usual Green's function for the medium of the form: (4) 4πjr − r 0 j Typical experiments that employ the angular-spectrum representation for analysis use a \parallel-plate" geometry (e.g. g.(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 s , 0) ). In this geometry the Green's function is conveniently written in the following form: where the z-axis is dened 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] ^ p where m = k 2 − (p 2 + q 2 ) is the complex propagator wavenumber and the imaginary 0 part of m is positive, i.e. =(m) > 0. With these denitions, the 2D Fourier Transform of the scattered eld in the detection plane is simply [33] Z sc (p, q, z d ) = dz 0 Ĝ 0 (p, q, z d , z 0 )T(p, q, z) .
where hat (^) indicates quantities in k-space. This equation is the fundamental equation for di raction tomography.

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 δµ 0 (r). Note that we need to know the s background eld at the object depth to accurately to obtain optical properties. Fig.(2)is a fowchart illustrating the algorithm described above.  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 eld at dierent depths is not measured and must be estimated using numerical models or a backpropagation technique [25,51]. This is the main diculty for obtaining 3D information in clinical situations and we are investigating various methods to overcome this problem [33,34].

Noise ampli cation at high spatial frequencies
In eq.(9) we divide two waves with complex wavenumbers (m and k 0 ). The scattered eld generally su ers from noise; the origin of this noise can be experimental, numerical, or due to nite sampling. The noise in the high spatial frequencies is preferentially amplied 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 amplication of the high frequency noise and the associated instabilities are common problems associated with Fourier convolution approaches to image reconstruction [25-28, [44][45][46]. Fig.(3) demonstrates this problem.
(a)Amplitude of the scattered wave as a function of p for a xed q and (b) amplitude of the tumor function plotted in the k-space for noiseless (only numerical noise) and (c) noise-added data from the single object (sec. 5.1) . The maximum frequency, jpj = π/x. The rising \wings" on the sides are due to noise. The noise eects are ampli ed by ˇ 10 3 relative to the noiseless case at large p.
One [37,38,41] way of dealing with this problem is to apply a low-pass spatial lter in order to suppress the high frequency components of the tumor function. The quantitative e ects 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 ltering modies image resolution (see eq. (12)) [31,41,43,44].

Filtering and normalization for improving depth information / image quality
We have investigated ltering and normalization techniques which provide dierent 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 o ers a prescription for extending the 2D projection approach to quasi-3D imaging and localization. In this section we briefy 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. We have found a very crude ltering 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,Ĝ 0 (p, q, z d , z 0 ); in particular we put Ĝ 0 (p, q, z d , z 0 ) ˇ ie i<(m)jz d −z 0 j . 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 g. (2)). Although quantitative information about optical properties is largely lost, the real-m lter has no free parameters (except from the optional windowing function) and provides superior information about the xyz central coordinates of heterogeneities.

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 lter enable us to accurately estimate the heterogeneity depth below the detection surface. Briefy, we rst obtain the tumor function for each slice (centered on some z j ) using the real-m lter. We then derive an image of δµ a (x, y, z j ) 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δµ a − δµ a j S j = q (11) P xy (δµ axy − δµ a ) 2 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 j is in essence a measure of the contrast of the object in the j th plane. We empirically found that the value of z j for which S j is a maximum closely approximates the actual depth of the object (z obj ). 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 5.1 .1).

High frequency ltering based on depth and experimental noise foor (G-ltering)
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 rst plot the raw data to determine the experimental noise level. For example by plotting sc vs. p and q (e.g. g. (3)) we will nd at suciently high p and q, that sc hits a \white" noise foor; the ratio of the signal amplitude at p = q = 0 to the signal amplitude at these high frequencies provides a measure of the experimental 0 signal-to-noise. Next we insert the estimated object depth z = z obj into the expression ^f or G 0 (p, q, z d , z 0 ) , and we determine at what value of m, G 0 (p, q, z d , z obj ) drops below a threshold based on our experimental S/N. In practice we set this threshold to be about ^ ve times the white noise foor. We set T = 0 for all m greater than this threshold (see g (2) for the occurance of this step in the general algorithm.). The result is a depth dependent pillbox lter applied to the tumor function. The rest of the reconstruction is the same. The combination technique o ers 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-lter technique. Usually the maximum S j determined with the G-lter is within ˇ1 cm of the maximum S j determined using the real-m lter (see section 3.1 ). In applying this latter procedure the high frequency cut-o depends on depth (z j ).

Data generation
Our \data" was derived from numerical simulations. In particular the nite dierence method with partial current boundary conditions was used to solve the diusion 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 foor 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 20x20x20cm box and use the central 9x9x5cm region for \experiments". The region is far from all boundaries and is a very good approximation to the innite 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 65x65 region in the detector plane at z j = z d = 5cm.
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 o 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 elds obtained in Born approximation.

Tissue phantoms with single slice heterogeneity
The simulated heterogeneity is shown in g. Single Slice Phantom: Left gure shows a 3D rendering of the phantom. Gray region has background properties. The detector plane is assumed to be at z=5cm and the source is at the origin in z=0cm plane. Amplitude of the scattered eld in the detector plane for the phantom shown in the middle (noiseless) and right (noise added) gures.

.1 Real-m imaging
We now employ the real-m lter 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 g. (5c), we plot S j vs z j ; we see that the contrast parameter is peaked at z j = 2.71cm, and we use this value to make a \true" real-m lter reconstruction of the data which is shown in g.(5d). This result is representative of all of the single-object phantoms that we tested.  ization 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 dierent depths will be addressed later in this paper.

.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 3x10 −3 for all cases. For z j = 2.71cm, the G-lter was set for p=q= 1.01/cm at maximum, and the resulting image in this plane is shown in g.(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 j using the G-lter (see g.(6b)). Notice that the peak value of S j has shifted from its true value. The G-lter image (at p=q=2.36/cm) based on this new z j is shown in g.(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 nd that this shift is related to the amount of power cut out by the lter. By changing the lter spatial frequency cut-o to a lower frequency (increasing the threshold) we nd that the shift is increased but the optical properties at the shifted location however approach their correct values. Finally we note that since the maximum spatial frequencies used in the G-lter 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 2 2π 2 2 p + q = .
(12) max max L In g. (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 e ect to note is that the resolution improves dramatically for objects near the detection plane. Estimate of resolution(cm) vs distance from source plane (cm) . The changing depth dependent cut-o frequency results in the increase in resolution with distance from source plane (i.e decrease in resolution with depth from the detector plane).

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 e ect 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 g. (8). The amplitude of the \measured" scattered eld at the detector plane is shown in its noiseless and noisy version in g.(8) for comparison. Two slice Phantoms: Two leftmost gures show 3D renderings. Gray region has background properties. The detector plane is at z=5cm and the source is at the origin in z=0cm plane. Amplitude of the scattered eld at the detector plane is shown in the two rightmost gures. The left is the noiseless data and the right shows the data after adding noise.
We performed essentially the same set of reconstruction procedures as described in sections 3. and 5.1 . In g.(9a) we show that S j has two maxima in dierent planes depending on the heterogeneity under consideration using real-m lter. Note that the real-m images were more noisy by comparison to the single object phantom. We again -0 nd that the optical properties are most accurate where S j has its shifted maxima with the G-lter. 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.   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 j for all transverse pixels (65x65= 4225 values in each projection) and plot all S j 's in three dimensions. In g.(10a) and (10b) we show an isosurface rendering at value S j = 0.042 at two viewing angles. Two dimensional projections in plane y=2cm and plane y=-2cm are also shown in g. (10c) and g.(10d) to provide a better sense of contrast. Comparing this to g. (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 e ects from Fourier domain reconstructions. These latter artifacts, however, are usually fairly easy to identify. We nd 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. x(cm)

Conclusion
We have studied projection images of turbid media based on the angular-spectrum representation. We demonstrated two ltering schemes. One of them, the real-m lter, provided object location information with no free parameters; the other method, the G-lter, introduced a systematic approach for eliminating high spatial frequencies from the reconstructions. Image localization in 3D was demonstrated using both lters 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 ltering schemes and noise foors, and optical property accuracy.