Frequency-domain methods [1],[2],[3],[4] are among the approaches being explored for use in localizing and characterizing inhomogeneous structures imbedded in turbid media. Most mathematical approaches to analyzing this problem have focussed on solving the Helmholtz equation with finite element methods using either the differential or integral equation. Recently, this problem has been explored in the context of diffraction tomography, with preliminary results indicating that significant speed improvements are possible [5] along with decreasing the number of views needed. [6] Diffraction tomographic theory is well developed for scenarios including arbitrary illumination and detection geometries, [7] near-field effects, [8] inhomogeneities imbedded in dispersive attenuating backgrounds, [9] and limited views. [10] However, the theory was originally developed for objects imbedded in non-attenuating background media with measurements in the far field. Extensions to the theory still rely upon the Fourier-domain insight obtained with the original theory. In this paper, we analyze the forward problem from a diffraction tomographic viewpoint that occurs when using diffuse photon density waves (DPDWs) in a turbid media. We show that significant insight into the forward problem is obtained by deriving the diffraction tomographic problem formulation from first principles because such a derivation explicitly includes effects caused by complex wave numbers. In particular, we show that, in many cases, the Fourier transform of the scattered wave, measured in a plane, is related to the three-dimensional Fourier transform of the object evaluated on a surface which is approximately a plane passing through zero spatial frequency.

For our problem formulation, we allow an arbitrary structure for the wavefront but assume that our measurements are made in a plane. We start our analysis by assuming that the inhomogeneities in the turbid media are small, purely absorptive, perturbations when compared to the background media, enabling us to employ the first Born solution [11] to the wave equation

which is denoted by u_B(r⃗) and is given by

where o(r⃗) is the absorptive inhomogeneity that we are trying to reconstruct, [1] u_o(r⃗) is the solution of the wave equation without the inhomogeneities, g(r⃗) is the Green’s function associated with the homogeneous infinite medium, k is the complex wavenumber of the homogeneous medium, and r⃗ = (x,y,z) is a three-dimensional spatial coordinate. For the case of u_o(r⃗) being a plane wave with a real wave number, it has been shown [7] that the Fourier transform of u_B(r⃗) is related to a curved surface of the Fourier transform of o(r⃗) which intersects the zero spatial frequency location and which, in two dimensions, is a semicircle. We desire to see what the surface geometry is for turbid media. We will do this by substituting in the angular spectrum form of the Green’s function into Eq.(2) and analyzing the result.

The Green’s function in Eq.(2) can be expressed as [12]

where the second equality is the angular spectrum representation of the Green’s function and the square root is determined by requiring the real part to be greater than zero. An interesting property of Eq.(3) is that the Green’s function has an analytic Fourier transform even though the Green’s function is singular at r⃗ = (0,0,0). This is a consequence of the nonzero imaginary component of k, which keeps the denominator from becoming zero because the spatial frequencies are real. Substituting Eq.(3) into Eq.(2) gives us

where

and where Re() denotes the real part and Im() denotes the imaginary part of a complex number. We desire to rewrite Eq.(4) so that the integration over r⃗′ can be viewed as a Fourier transform of all the functions which depend upon r⃗′. To this end, we reexpress Eq.(4) as

where we have explicitly written the integration over r⃗′ in terms of x’, y’, and z’, and we have interchanged the order of integration, integrating over the spatial coordinates before the spatial frequency coordinates. We have lumped the exponential term which depends upon γ_r with the object and source terms, and have explicitly displayed our measurement geometry by indicating that the scattered field is measured in the z=z_o plane. We also assume an observation geometry where all the sources and inhomogeneities are at z values less than z_o. This enables us to replace |z_o - z′| with its argument. We have included this simplification for the terms multiplying γ_i but have retained the absolute value function for the terms multiplying γ_r in order to more closely match Fourier transform pairs as listed in standard tables. [13]

It can be seen in Eq.(6) that the inner integral term is the three-dimensional Fourier transform of the object function multiplied by both the illumination function and an exponential term which results from the complex nature of k. Denoting the Fourier transform of this product by O_uγ, we can rewrite Eq.(6) as

Our interest is in the two-dimensional Fourier-domain properties of u_B(x,y,z_o). Specifically, we desire to know the relationship between the two-dimensional Fourier transform of u_B(x,y,z_o), U_B(ω_x,ω_y,z_o), and the three-dimensional Fourier transform O_uγ(ω_x,ω_y,ω_z). Therefore, we Fourier transform Eq.(7) and interchange the order of integration, resulting in

The inner integrals are separable and evaluate to 4π²δ(ω_x-α_x)δ(ω_y-α_y). Thus Eq.(8) becomes

which is the main result of the paper. To reconstruct the object spectrum from the measured data as given in Eq.(9), it is necessary to either use a filtered backpropagation algorithm modified from the standard diffraction tomographic approach [7] or use Fourier interpolation after deconvolving the illumination and attenuation functions (see Eq.(6)). We are developing both approaches at this time.

From Eqs.(5) and (9) it can be seen that O_uγ (ω_x ,ω_y ,-γ_i ) is a portion of O_uγ(ω_x,ω_y,ω_z) evaluated in a two-dimensional region. We must analyze the γ_i term to characterize this region further. From Eq.(5) we have

For DPDWs, the functional form of the k² term is given by [14]

where v is the speed of light in the turbid medium, f_t is the temporal frequency of the DPDW, μ_a is the absorption coefficient and μ′_s is the reduced scattering coefficient. We see from Eq.(11) that the real part of k² is always negative and its imaginary part is always positive. As a result, the expression $\sqrt{{\omega }_x^2+{\omega }_y^2-k^2}$ has a negative imaginary part whose absolute value is always less than its real part and which is dominated by the real part as spatial frequency becomes arbitrarily large. It is instructive to compare the standard diffraction tomographic curve in two-dimensional Fourier space with the curve generated by γ_i. We generate a two-dimensional curve by setting ω_x=0 in the expression for γ_i. The results are shown in Figure 1, where we have plotted three curves: one where Re{k²}<<Im{k²}, one where Re{k²}>>Im{k²}, and one where k is real and positive, as it is in standard diffraction tomography. All three curves have been normalized so that |k| = 1. Notice the familiar band-limited circular shape of the curve for real and positive k. It does not intersect zero spatial frequency because we have not included plane-wave illumination effects for the plots. Notice that the two curves corresponding to DPDWs are not bandlimited. The curve where Re{k²}<<Im{k²} is more flattened than the curve for real k, and the curve where Re{k²}>>Im{k²} is essentially a straight line. From Eq.(11), we see that the third curve results when either the absorption coefficient of the medium is large or when the temporal frequency of the DPDW is small. For many types of human tissue, [15] μ_a is on the order of 0.5 cm^-1, μ′_s is on the order of 10 cm^-1, and v is on the order of 2×10¹⁰ cm/s. For these values, k² is given by

 

Figure 1 Plots of the two-dimensional projection of the surface on which the Fourier transform of the convolved object function is obtained with diffraction tomography. The dotted line is for real values of k, the dashed line is for Re{k²}<<Im{k²}, and the solid line is for Re{k²}>>Im{k²}. The horizontal axis is ω_x and the vertical axis is ω_y.

Download Full Size | PDF

where the temporal frequency ft in Eq.(12) is to be expressed in GHz units. Therefore, for DPDWs which are generated using a few hundreds of megahertz, we have Re{k²}>>Im{k²}. We must have temporal frequencies of many gigahertz before Re{k²}<<Im{k²}. Recent results [16] for human breast tissue show that typical values of μ_a are around 0.05cm^-1, which indicates that Im{k²} becomes significant at temporal frequencies of tens of megahertz instead of hundreds of megahertz for the larger values of μ_a.

Because the surface which is defined by -γ_i does not intersect the zero spatial frequency point, a tomographic reconstruction using multiple views will result in a band-limited high-pass version of the object. In standard diffraction tomography, this problem is overcome by using plane wave illumination. The convolution of the Fourier transform of the plane wave with the object (see Eq.(6)) shifts the surface along the z-axis in Fourier space [11] to include the zero spatial frequency point, permitting a complete low-pass-filtered version of the object’s Fourier transform to be reconstructed using multiple views. Conceptually, this can be done with DPDWs as well. However, DPDWs are typically created using a few point illumination sources, which do not provide this convenient Fourier-domain shift. However, looking again at Eq.(6), we see that the object function is multiplied by an exponential function whose Fourier transform G(ω_x,ω_y,ω_z) is given by [13]

where we have set z_o=0 because its value depends only upon the location of our coordinate system. Because O(ω_x,ω_y,ω_z) is convolved with G(ω_x,ω_y,ω_z), the correlation scale of G(ω_x,ω_y,ω_z) (that is, its smoothness) is imposed upon O(ω_x,ω_y,ω_z). When the imaginary part of k² is much smaller than the real part of k², the correlation scale of O(ω_x,ω_y,ω_z) convolved with G(ω_x,ω_y,ω_z) along the ω_z axis is much larger than γ_i. The implication of this fact is that we have

and so the Fourier transform of the measured scattered field is related to the convolution of the Fourier transforms of the object, the illumination, and the exponential function evaluated at the ω_z=0 plane. This result is similar to results from straight ray tomography. The approximation shown in Eq.(14) becomes more accurate for lower temporal frequencies and higher absorption coefficients. Furthermore, because |γ_i| is always less than γ_r and because the correlation scale of G(ω_x,ω_y,ω_z) along the ω_z axis is set by γ_r, we will always have significant correlations between the spatial frequency components along the curve defined by -γ_i and the spatial frequency components along the ω_z=0 axis. Because of this fact, Eq.(14) should always be qualitatively true. This correlation effect is similar to that exploited by synthetic aperture radar tomographic reconstruction techniques in which a small offset portion of the object’s Fourier transform is used, yet which yield surprisingly good reconstructed images. [17]

In summary, we have applied a diffraction tomographic problem formulation to DPDWs in turbid media. We have shown that the measured scattered field’s Fourier transform is related to a convolved version of the object’s Fourier transform evaluated on a Fourier-domain surface which is, in many cases, well approximated by a plane passing through zero spatial frequency. As a result, we can use Fourier-domain concepts to analyze the quality of image reconstructions using concepts such as limited-view restrictions [10] and signal-to-noise resolution limits. [18]