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

Turgut Durduran; J. P. Culver; M. J. Holboke; X. D. Li; Leonid Zubkov; B. Chance; Deva N. Pattanayak; A. G. Yodh

doi:10.1364/OE.4.000247

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]:

Φ_{sc} (r) = \int_{V} T (r') G_{0} (r, r') d^{3} r' .

Figure 1. The generic near-field diffraction 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.

Download Full Size | PDF

Here T(r′) is “tumor” function at position r′ which, in the Born Approximation, depends on the background diffuse photon density wave, Φ₀(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 ₀, the background diffuse photon density wavevector $k_{0} = \sqrt{(- υ μ_{a 0} + i ω) / D_{0})}$ , and the modulation frequency, ω. It can be conveniently separated into absorptive and scattering parts [49], i.e.,

T_{abs} (r) = - \frac{υ}{D_{0}} Φ_{0} (r) δ μ_{a} (r),

T_{sc} (r) = \frac{3 D_{0} k_{0}^{2}}{υ} Φ_{0} (r) .

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

G_{0} (r, r') = \frac{\exp (i k_{0} ∣ r - r' ∣)}{4 π ∣ r - r' ∣} .

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)+ Φ₀(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:

G_{0} (r, r') = \int \int dpdq {\hat{G}}_{0} (p, q, z_{d}, z') e^{- i (p (x_{d} - x') + q (y_{d} - y'))},

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]

{\hat{G}}_{0} (p, q, z_{d}, z') = \frac{i}{2 m} e^{i m ∣ z_{d} - z' ∣},

where $m = \sqrt{k_{0}^{2} - (p^{2} + q^{2})}$ 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]

{\hat{Φ}}_{sc} (p, q, z_{d}) = \int d z' {\hat{G}}_{0} (p, q, z_{d}, z') \hat{T} (p, q, z) .

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

{\hat{Φ}}_{sc} (p, q, z_{d}) = \sum_{j = 1}^{N} {\hat{G}}_{0} (p, q, z_{d}, z') \hat{T} (p, q, z_{j}) = \sum_{j = 1}^{N} \frac{i Δ z}{2 m} \hat{T} (p, q, z_{j}) e^{i m (z_{d} - z_{j})} .

∆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

\hat{T} (p, q, z_{obj}) = \frac{Φ_{sc} (p, q, z_{d})}{Δ z {\hat{G}}_{0} (p, q, z_{d}, z_{obj})},

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

T (x', y', z_{obj}) = \int \int dpdq \frac{e^{- i 2 π} (p x' + q y')}{Δ z \hat{G} (p, q, z_{obj}, z)} \int \int dxdy e^{- i 2 π (p x' + q y')} Φ_{sc} (x, y, z_{d}),

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.

Figure 2. Simplified flow diagram of the image reconstruction algorithm. Dotted lines are used for optional steps. Brown and blue indicate real-m filtering, green indicates G-filtering steps.

Download Full Size | PDF

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 ₀). 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.

Figure 3. (a)Amplitude of the scattered wave as a function of p for a fixed 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, |p| = π/∆x. The rising “wings” on the sides are due to noise. The noise effects are amplified by ≈ 10³ relative to the noiseless case at large p.

Download Full Size | PDF

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̂₀ (p, q, Z_d, z′); in particular we put Ĝ₀(p, q, Z_d, z′) ≈ ie^{iℜ
(m)|z_d-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_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

S_{j} = \frac{∣ δ μ_{a} - δ {\bar{μ}}_{a} ∣}{\sqrt{\sum_{xy} {(δ μ_{axy} - δ {\bar{μ}}_{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).

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̂₀(p, q, z_d, z′) , and we determine at what value of m, Ĝ₀(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 five times the white noise floor. We set 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_j determined with the G-filter is within ≈1 cm of the maximum S_j determined using the real-m filter (see section 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)³. 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_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° 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 μ_ao = 0.02 cm^-1 and μ_so = 8 cm^-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.

Figure 4. Single Slice Phantom: Left figure 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 field in the detector plane for the phantom shown in the middle (noiseless) and right (noise added) figures.

Download Full Size | PDF

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_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 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.

Figure 5. (a) and (b) projections at z=2.42, z=3.28 respectively, (c) S_j vs z_j (cm) through the transverse center, peaks at z=2.71, (d) projection at z=2.71. All with real-m filter.

Download Full Size | PDF

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.

Figure 6. (a) projection at z=2.71 at real-m peak, (b) S_j vs z_j (cm) through the transverse center peaks at z=4cm , (c) projection at z=4.0. All obtained with G-giltering.

Download Full Size | PDF

We next computed S_j vs z_j using the G-filter (see fig.(6b)). Notice that the peak value of S_j has shifted from its true value. The G-filter image (at p=q=2.36/cm) based on this new z_j 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.

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

p_{max}^{2} + q_{max}^{2} = {(\frac{2 π}{L})}^{2} .

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.

Figure 7. Estimate of resolution(cm) vs distance from source plane (cm) . The changing depth dependent cut-off frequency results in the increase in resolution with distance from source plane (i.e decrease in resolution with depth from the detector plane).

Download Full Size | PDF

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.

Figure 8. Two slice Phantoms: Two leftmost figures 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 field at the detector plane is shown in the two rightmost figures. The left is the noiseless data and the right shows the data after adding noise.

Download Full Size | PDF

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_j 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.

Figure 9. (a) Circles (crosses) show S_j vs z_j (cm) through the center region of deeper (shallower) object obtained with m-filter. Peak for both objects are exact. Then with G-filter we get projections at (b) z=3.71 and (c) at z=4.41.

Download Full Size | PDF

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 (65×65= 4225 values in each projection) and plot all S_j’s in three dimensions. In fig.(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 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.

Figure 10. Three Dimensional Rendering of the G-filter reconstruction. An isosurface at S_j = 0.042 is shown in two different angles in (a) and (b). In (c) and (d) images of S_j in x-z plane through y=2cm and y=-2cm are shown respectively. Compare the results to that of fig.(8)

Download Full Size | PDF

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).

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

Abstract

1. Introduction