## Abstract

We present a model of the forward problem for diffuse photon density waves in turbid medium using a diffraction tomographic problem formulation. We consider a spatially-varying inhomogeneous structure whose absorption properties satisfy the Born approximation and whose scattering properties are identical to the homogeneous turbid media in which it is imbedded. The two-dimensional Fourier transform of the scattered field, measured in a plane, is shown to be related to the three-dimensional Fourier transform of the object evaluated on a surface which in many cases is approximately a plane.

© Optical Society of America

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]

$$=\frac{1}{8{\pi}^{2}}\int \int \frac{1}{\sqrt{{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}}}\mathrm{exp}\{-\mid z\mid \sqrt{{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}}+\mathrm{ix}{\alpha}_{x}+\mathrm{iy}{\alpha}_{y}\}{\mathrm{d\alpha}}_{x}{\mathrm{d\alpha}}_{y}$$

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

$$\times \sqrt{{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}}+i\left(x-\mathrm{x\prime}\right){\alpha}_{x}+i\left(y-\mathrm{y\prime}\right){\alpha}_{y}\}{\mathrm{d\alpha}}_{x}{\mathrm{d\alpha}}_{y}d\overrightarrow{\mathrm{r\prime}}$$

$$=\frac{1}{8{\pi}^{2}}\int o\left(\overrightarrow{r\prime}\right){u}_{o}\left(\overrightarrow{r}\prime \right)\int \int \frac{1}{\gamma}\mathrm{exp}\{-\mid z-\mathrm{z\prime}\mid \left({\gamma}_{r}+i{\gamma}_{i}\right)$$

$$+i\left(x-\mathrm{x\prime}\right){\alpha}_{x}+i\left(y-\mathrm{y\prime}\right){\alpha}_{y}\}{\mathrm{d\alpha}}_{x}{\mathrm{d\alpha}}_{y}d\overrightarrow{\mathrm{r\prime}}$$

where

$$=\mathrm{Re}\left\{\sqrt{{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}}\right\}+i\phantom{\rule{.2em}{0ex}}\mathrm{Im}\left\{{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}\right\}$$

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

$$\times \int \int \int o(\mathrm{x\prime},\mathrm{y\prime},\mathrm{z\prime}){u}_{o}(\mathrm{x\prime},\mathrm{y\prime},\mathrm{z\prime})\mathrm{exp}\left\{-\mid {z}_{o}-\mathrm{z\prime}\mid {\gamma}_{r}\right\}$$

$$\times \mathrm{exp}\left\{-\mathrm{ix\prime}{\alpha}_{x}-\mathrm{iy}\prime {\alpha}_{y}-\mathrm{iz\prime}\left(-{\gamma}_{i}\right)\right\}\mathrm{dx}\prime \mathrm{dy}\prime \mathrm{dz\prime}{\mathrm{d\alpha}}_{x}{\mathrm{d\alpha}}_{y}$$

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

$$\times \int \int \mathrm{exp}\left\{i\left(x{\alpha}_{x}+y{\alpha}_{y}\right)\right\}\mathrm{exp}\left\{-i\left({\mathrm{x\omega}}_{x}+{\mathrm{y\omega}}_{y}\right)\right\}\mathrm{dxdyd}{\alpha}_{x}d{\alpha}_{y}$$

The inner integrals are separable and evaluate to
4π^{2}δ(ω_{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^{2} 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^{2} 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^{2}}<<Im{k^{2}}, one where
Re{k^{2}}>>Im{k^{2}}, 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^{2}}<<Im{k^{2}} is more flattened than
the curve for real k, and the curve where
Re{k^{2}}>>Im{k^{2}} 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^{10} cm/s. For these
values, k^{2} is given by

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^{2}}>>Im{k^{2}}. We must have temporal
frequencies of many gigahertz before
Re{k^{2}}<<Im{k^{2}}. Recent results [16] for human breast tissue show that typical values of
μ_{a} are around 0.05cm^{-1}, which indicates that
Im{k^{2}} 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^{2} is much smaller than the real part of k^{2},
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]

## References and Links

**
1
. **
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] [PubMed]

**
2
. **
H.B.
Jiang
,
K.D.
Paulsen
,
U.L.
Osterberg
, and
M.S.
Patterson
, “
Frequency-domain optical-image reconstruction in turbid media - an experimental study of single-target detectability
,”
Appl. Opt.
**
36
**
,
52
–
63
(
1997
) [CrossRef] [PubMed]

**
3
. **
S.A.
Walker
,
S.
Fantini
, and
E.
Gratton
, “
Image reconstruction using back-projection from frequency-domain optical measurements in highly scattering media
,”
Appl. Opt.
**
36
**
,
170
–
179
(
1997
) [CrossRef] [PubMed]

**
4
. **
S.B.
Colak
,
D.G.
Papioannou
,
G.W.
Hooft
, and
M.B. van der
Mark
, “
Optical image reconstruction with deconvolution in light diffusing media
,” in
*
Photon Propagation in Tissues
*
,
B.
Chance
,
D.T.
Delpy
, and
G.J.
Mueller
, eds.,
Proc. Soc. Photo-Opt. Instrum. Eng.
**
2626
**
,
306
–
315
(
1995
)

**
5
. **
X.D.
Li
,
T.
Durduran
,
A.G.
Yodh
,
B.
Chance
, and
D.N.
Pattanayak
, “
Diffraction tomography for biochemical imaging with diffuse-photon density waves,
”
Opt. Lett.
**
22
**
,
573
–
575
(
1997
) [CrossRef] [PubMed]

**
6
. **
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
–
220
(
1997
) [CrossRef] [PubMed]

**
7
. **
A.J.
Devaney
, “
Reconstructive tomography with diffracting wavefields
,”
Inverse Problems
**
2
**
,
161
–
183
(
1986
) [CrossRef]

**
8
. **
A.
Schatzberg
and
A.J.
Devaney
, “
Super-resolution in diffraction tomography
,”
Inverse Problems
**
8
**
,
149
–
164
(
1992
) [CrossRef]

**
9
. **
A.J.
Devaney
, “
Linearised inverse scattering in attenuating media
,”
Inverse Problems
**
3
**
,
389
–
397
(
1987
) [CrossRef]

**
10
. **
A.J.
Devaney
, “
The limited-view problem in diffraction tomography
,”
Inverse Problems
**
5
**
,
501
–
521
(
1989
) [CrossRef]

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

**
12
. **
A.
Baños
Jr.
,
*
Dipole Radiation in the Presence of a Conducting Half-Space
*
(
Pergamon Press, Oxford
,
1966
)

**
13
. **
K.B.
Howell
, “
Fourier transforms
,” in
*
The Transforms and Applications Handbook
*
,
A.D.
Poularikas
, ed. (
CRC Press, Boca Raton
,
1996
)

**
14
. **
D.A.
Boas
,
M.A.
O’Leary
,
B.
Chance
, and
A.G.
Yodh
, “
Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: Analytic solution and applications
,”
Proc. Natl. Acad. Sci. USA
**
91
**
,
4887
–
4891
(
1994
) [CrossRef] [PubMed]

**
15
. **
W.F.
Cheong
,
S.A.
Prahl
, and
A.J.
Welch
, “
A review of the optical properties of biological tissues
,”
IEEE J. Quantum Electronics
**
26
**
,
2166
–
2185
(
1990
) [CrossRef]

**
16
. **
H.
Heusmann
,
J.
Koelzer
, and
G.
Mitic
, “
Characterization of female breasts in vivo by time-resolved and spectroscopic measurements in near infrared spectroscopy
,”
J.Biomed.Opt.
**
1
**
,
425
–
434
(
1996
) [CrossRef]

**
17
. **
D.C.
Munson
Jr.
and
J.L.C.
Sanz
, “
Image reconstruction from frequency-offset Fourier data
,”
Proc. IEEE
**
72
**
,
661
–
669
(
1984
) [CrossRef]

**
18
. **
C.L.
Matson
,
I.A.
Delarue
,
T.M.
Gray
, and
I.E.
Drunzer
, “
Optimal Fourier spectrum estimation from the bispectrum
,”
Computers Elect. Engng
**
18
**
,
485
–
497
(
1992
) [CrossRef]