Time reversal optical tomography: locating targets in a highly scattering turbid medium

doi:10.1364/OE.19.021956

« Show journal navigation

Time reversal optical tomography: locating targets in a highly scattering turbid medium

Binlin Wu, W. Cai, M. Alrubaiee, M. Xu, and S. K. Gayen »View Author Affiliations

Optics Express, Vol. 19, Issue 22, pp. 21956-21976 (2011)
http://dx.doi.org/10.1364/OE.19.021956

View Full Text Article

Acrobat PDF (1359 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (61)
Cited By
Figures (6)
Metrics

Abstract

A time reversal optical tomography (TROT) method for near-infrared (NIR) diffuse optical imaging of targets embedded in a highly scattering turbid medium is presented. TROT combines the basic symmetry of time reversal invariance and subspace-based signal processing for retrieval of target location. The efficacy of TROT is tested using simulated data and data obtained from NIR imaging experiments on absorptive and scattering targets embedded in Intralipid-20% suspension in water, as turbid medium. The results demonstrate the potential of TROT for detecting and locating small targets in a turbid medium, such as, breast tumors in early stages of growth.

1. Introduction

Optical imaging of targets embedded in a highly scattering turbid medium, such as, a tumor in a breast, is a challenging problem because light is strongly absorbed and scattered by the medium leading to poor signal-to-noise ratio, as well as, loss of phase coherence and polarization. As a consequence distinct, sharp image of the targets may not be formed directly. Various frequency-domain, time-resolved, and steady-state inverse image reconstruction (IIR) [1

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]

–5

W. Cai, S. K. Gayen, M. Xu, M. Zevallos, M. Alrubaiee, M. Lax, and R. R. Alfano, “Optical tomographic image reconstruction from ultrafast time-sliced transmission measurements,” Appl. Opt. 38(19), 4237–4246 (1999). [CrossRef] [PubMed]

] approaches are being pursued to form tomographic images using diffusively scattered light measured at the sample boundary. IIR is an ill-posed problem and the development of reliable and fast approaches remains a formidable task. Recent IIR algorithms, such as Newton-Raphson-Marquardt algorithms [6

B. A. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron. 9(2), 199–209 (2003). [CrossRef]

] and direct linear inversion of 3-D matrices [7

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

], are time consuming. The iterative methods [7

S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993). [CrossRef] [PubMed]

] may not ensure that the obtained result arrives at a “global minimum” or converges to a “local minimum”. Still the potential for developing non-invasive imaging approaches with diagnostic ability motivates the ongoing diffuse optical tomography (DOT) research using NIR light.

Many applications require rather accurate determination of location of target(s) in three dimensions. For example, a recent study involving 35,319 patients underscores the influence of primary tumor location on breast cancer prognosis [9

N. Kroman, J. Wohlfahrt, H. T. Mouridsen, and M. Melbye, “Influence of tumor location on breast cancer prognosis,” Int. J. Cancer 105(4), 542–545 (2003). [CrossRef] [PubMed]

], and makes it imperative that DOT for breast cancer detection be able to obtain three-dimensional (3-D) location of the tumor. While two-dimensional (2-D) IIR approaches may provide only lateral positions, 3-D IIR approaches attempt to retrieve all three position coordinates of the target(s). Various frequency-domain, time-domain, and steady-state DOT approaches have addressed the target localization problem with different measures of success [1

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]

–8

S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993). [CrossRef] [PubMed]

]. Several groups have paid particular attention to retrieving target location. Kepshire et al. developed a subsurface DOT approach to obtain location information of absorbing and fluorescent targets, but observed the sensitivity to vary nonlinearly with depth [10

D. S. Kepshire, S. C. Davis, H. Dehghani, K. D. Paulsen, and B. W. Pogue, “Subsurface diffuse optical tomography can localize absorber and fluorescent objects but recovered image sensitivity is nonlinear with depth,” Appl. Opt. 46(10), 1669–1678 (2007). [CrossRef] [PubMed]

]. Mohajerani et al. reported a fluorescent tomography method for locating fluorescent targets embedded in a heterogeneous medium using partitioning of the fluorophore distribution into an object subspace and a background subspace [11

P. Mohajerani, A. A. Eftekhar, and A. Adibi, “Object localization in the presence of a strong heterogeneous background in fluorescent tomography,” J. Opt. Soc. Am. A 25(6), 1467–1479 (2008). [CrossRef] [PubMed]

]. Godavarty et al. developed another fluorescent tomography approach that used a hemispherical breast phantom, near-infrared light-induced fluorescence from a contrast agent, and finite element method-based reconstruction algorithms to obtain target information up to a depth of 2 cm from breast phantom surface [12

A. Godavarty, A. B. Thompson, R. Roy, M. Gurfinkel, M. J. Eppstein, C. Zhang, and E. M. Sevick-Muraca, “Diagnostic imaging of breast cancer using fluorescence-enhanced optical tomography: phantom studies,” J. Biomed. Opt. 9(3), 488–496 (2004). [CrossRef] [PubMed]

]. Zhao et al. introduced a layer-based sigmoid adjustment method to improve depth resolution of DOT and achieved positioning error within 3 mm for depths from 10 to 30 mm [13

Q. Zhao, L. Ji, and T. Jiang, “Improving depth resolution of diffuse optical tomography with a layer-based sigmoid adjustment method,” Opt. Express 15(7), 4018–4029 (2007). [CrossRef] [PubMed]

]. Optical tomography using independent component analysis (OPTICA) approach developed by Xu et al. uses multi-source probing and multi-detector signal acquisition scheme and a numerical algorithm based on independent component analysis of information theory to obtain 3-D position information of absorbing, scattering and fluorescent targets embedded in highly scattering turbid media, and “model breast” assembled using ex vivo human breast tissue [14

M. Alrubaiee, M. Xu, S. K. Gayen, M. Brito, and R. R. Alfano, “Three-dimensional optical tomographic imaging of scattering objects in tissue-simulating turbid medium using independent component analysis,” Appl. Phys. Lett. 87(19), 191112 (2005). [CrossRef]

–17

M. Xu, M. Alrubaiee, S. K. Gayen, and R. R. Alfano, “Optical diffuse imaging of an ex vivo model cancerous human breast using independent component analysis,” IEEE J. Sel. Top. Quantum Electron. 14(1), 43–49 (2008). [CrossRef]

]. Co-registration approaches that use another modality, such as, ultrasound, magnetic resonance imaging, and x-ray mammography for locating suspect areas and DOT for obtaining images have also been introduced [18

A. Poellinger, J. C. Martin, S. L. Ponder, T. Freund, B. Hamm, U. Bick, and F. Diekmann, “Near-infrared laser computed tomography of the breast first clinical experience,” Acad. Radiol. 15(12), 1545–1553 (2008). [CrossRef] [PubMed]

–21

A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42(25), 5181–5190 (2003). [CrossRef] [PubMed]

In this article we report on the development of a time reversal optical tomography (TROT) [22

W. Cai, M. Alrubaiee, S. K. Gayen, M. Xu, and R. R. Alfano, “Three-dimensional optical tomography of objects in turbid media using the round-trip matrix,” Proc. SPIE 5693, 4–9 (2005). [CrossRef]

–24

B. Wu, W. Cai, M. Alrubaiee, M. Xu, and S. K. Gayen, “Three dimensional time reversal optical tomography,” Proc. SPIE 7892, 78920G (2011). [CrossRef]

] approach for NIR optical imaging of target(s) in a turbid medium, and present initial results of its efficacy using both simulated and experimental data.

Time reversal (TR) invariance, the basic symmetry that commonly holds in microscopic physics, forms the basis for macroscopic TR imaging. TR imaging using the so-called “time-reversal mirrors” (TRMs) has been used as an experimental tool in acoustics with practical applications in medicine, underwater imaging, and nondestructive testing [25

C. Prada, F. Wu, and M. Fink, “The iterative time reversal mirror: a solution to self-focusing in the pulse echo mode,” J. Acoust. Soc. Am. 90(2), 1119–1129 (1991). [CrossRef]

–28

M. Fink, “Time reversal mirrors,” J. Phys. D Appl. Phys. 26(9), 1333–1350 (1993). [CrossRef]

]. The theoretical and numerical techniques involved in time reversal have been used for applications involving both acoustic waves and electromagnetic waves (radar) [28

M. Fink, “Time reversal mirrors,” J. Phys. D Appl. Phys. 26(9), 1333–1350 (1993). [CrossRef]

–33

G. Lerosey, J. de Rosny, A. Tourin, and M. Fink, “Focusing beyond the diffraction limit with far-field time reversal,” Science 315(5815), 1120–1122 (2007). [CrossRef] [PubMed]

Devaney and associates developed a theoretical framework for a TR imaging method with Multiple Signal Classification (MUSIC) for finding the location of scattering targets whose size is smaller than the wavelength of acoustic waves or electromagnetic waves (radar) used for probing the homogeneous or inhomogeneous background medium in which the targets were embedded [34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,35

H. Lev-Ari and A. J. Devaney, “The time reversal techniques re-interpreted: subspace-based signal processing for multi-static target location,” in Proceedings of the 1st IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM '00), (Cambridge, MA, USA, 2000) pp. 509–513.

]. While their initial focus was on back-propagation geometry that used coincident acoustic or electromagnetic transceiver array for interrogating the targets, they later extended the formalism to transmission geometry where sources and detectors were distinct and separated [36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

]. They also generalized the theory which was based on distorted wave Born approximation (DWBA) to account for multiple scattering between the targets [37

F. K. Gruber, E. A. Marengo, and A. J. Devaney, “Time-reversal imaging with multiple signal classification considering multiple scattering between the targets,” J. Acoust. Soc. Am. 115(6), 3042–3047 (2004). [CrossRef]

]. In its basic form TR-MUSIC found target location from knowledge of the response matrix K, which was constructed from multi-static data collected by the transceiver array [34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,35

]. TR-MUSIC provided higher spatial resolution than the conventional TR imaging, especially in the case where targets were not well resolved [34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,35

,38

C. Prada and J. L. Thomas, “Experimental subwavelength localization of scatterers by decomposition of the time reversal operator interpreted as a covariance matrix,” J. Acoust. Soc. Am. 114(1), 235–243 (2003). [CrossRef] [PubMed]

We are adapting and extending the TR-MUSIC approach to the optical domain, i.e. to diffusive optical imaging for detecting and locating targets embedded in a turbid medium. In this paper, TROT is studied in details using both simulated data and data from transillumination NIR imaging experiments in slab geometry. A TR matrix is obtained by multiplying the response matrix formed using experimental or simulated data to its conjugate matrix. The leading non-zero eigenvalues of the Hermitian TR matrix determine the signal subspace due to presence of the targets. The signal subspace is separated from the noise subspace using an L-curve method [5

,39

P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. 34(4), 561–580 (1992). [CrossRef]

,40

E. A. Marengo, F. K. Gruber, and F. Simonetti, “Time-reversal MUSIC imaging of extended targets,” IEEE Trans. Image Process. 16(8), 1967–1984 (2007). [CrossRef] [PubMed]

]. The vector subspace method, MUSIC, along with Green’s functions calculated from an appropriate forward model for light propagation through the turbid medium is then used to determine the locations of the targets. The MUSIC algorithm judges if the calculated Green’s function vector corresponding to a location in the sample is mapped into the signal subspace or the noise subspace.

Several salient features make TROT attractive and potentially more promising than other IIR methods. First the size of the TR matrix is much smaller than those used in other IIR approaches, which makes solution of the eigenvalue problem easier and faster. Second, to determine locations of targets, TR-MUSIC approach runs the program over all voxels only once, and there is no need to carry out an iterative procedure done by other inverse approaches. Other IIR approaches seek to determine the absorption and scattering parameters at all voxels into which the sample is divided. The process is iterative, computationally intensive, and leads to a solution of the inverse problem that is not unique because the problem is ill-posed, even when there is no noise. In contrast, TROT seeks to determine the locations of the targets first and thereafter retrieve other information, such as, the size and optical properties of the limited number of targets in the medium, which requires significantly less computation time. The focus of this paper is on finding the locations of targets.

Our result using simulated data shows that without the presence of noise TROT determines the locations of the embedded targets accurately with high resolution. TROT exhibits promise to locate targets both in simulations and experiments even when substantial noise is present. Images of small targets obtained by this approach are sharper than that obtained by other IIR approaches.

This paper is organized as follows. In section 2, the formalism of the TROT approach is presented. In section 3, the numerical algorithm of TROT is described. In section 4, the efficacy of the formalism is tested using simulated data. Section 5 presents the results when the formalism is applied to experimental data using Intralipid-20% suspension in water as the highly scattering turbid medium. Section 6 discusses the results.

2. Formalism

2.1 Diffusion approximation, perturbation method and response matrix

The starting point for the TROT formalism is the diffusion approximation [41

A. Ishimaru, “Diffusion of a pulse in densely distributed scatterers,” J. Opt. Soc. Am. 68(8), 1045–1050 (1978). [CrossRef]

–43

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(12), 2331–2336 (1989). [CrossRef] [PubMed]

] of the radiative transfer equation (RTE) [44

S. Chandrasekhar, Radiative Transfer (Clarendon Press, Oxford, 1950).

,45

A. Ishimaru, Wave Propagation and Scattering in Random Media, Volume 1: Single Scattering and Transport Theory (Academic, New York, 1978).

]. The perturbation in the light intensity distribution due to small inhomogeneities (targets) embedded in a homogeneous medium, to the first order Born approximation, can be written as [46

S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25(12), 123010 (2009). [CrossRef]

,47

S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef] [PubMed]

]

\begin{matrix} Δ ϕ (r_{d}, r_{s}) = - \int^{​} G (r_{d}, r) δ μ_{a} (r) c G (r, r_{s}) d^{3} r \\ - \int^{​} δ D (r) c \nabla_{r} G (r_{d}, r) \cdot \nabla_{r} G (r, r_{s}) d^{3} r, \end{matrix}

(1)

where r_s, r_d, and r are the positions of a point-like source of unit power, detector and target, respectively;

G (r, r_{s})

and

G (r_{d}, r)

are the Green’s functions that describe light propagations from the source to the target and from the target to the detector, respectively; δµ_a is the difference in absorption coefficient and δD is the difference in diffusion coefficient between the targets and the background medium; and c is the light speed in the medium.

A multi-source interrogation and multi-detector signal acquisition scheme is used to acquire transillumination data, from which the difference in the light intensity distribution due to the targets,

Δ ϕ = ϕ - ϕ_{0}

, is found, where

ϕ

is the light intensity distribution measured on the sample boundary with targets embedded in the scattering medium and

ϕ_{0}

is ideally the light intensity distribution without the targets, which in practice is approximated by an “average” over all the multi-source measurements. A response matrix K is constructed with

- Δ ϕ

, to describe the transport of light from different sources through the embedded objects to the array of detectors [22

,36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

For small, point-like absorptive targets, the matrix elements can be rewritten in a discrete form as:

K_{i j} = \sum_{m = 1}^{M} G^{d} (r_{i}, X_{m}) τ_{m} G^{s} (X_{m}, r_{j}), i = 1, 2, \dots, N_{d}; j = 1, 2, \dots, N_{s},

(2)

where

τ_{m} = δ μ_{a} (X_{m}) c δ V_{m}

is the optical absorption strength of the m^th target, δV_m is the volume of m^th target, r_i, r_j and X_m are locations of the i^th detector, j^th source and m^th target, respectively. Due to the reciprocity of light propagation in the medium,

G (r, r') = G (r', r)

. Thus,

K_{i j} = \sum_{m = 1}^{M} G^{d} (X_{m}, r_{i}) τ_{m} G^{s} (X_{m}, r_{j}),

(3)

and

K = {K_{i j}} = \sum_{m = 1}^{M} g_{d} (X_{m}) τ_{m} g_{s}^{T} (X_{m}),

(4)

where

g_{s} (r)

and

g_{d} (r)

are Green’s function vectors (GFVs) associated with the source array and detector array, respectively. GFVs are defined as

g_{s} (r) = {[G^{s} (r_{1}, r), G^{s} (r_{2}, r), \dots, G^{s} (r_{N_{s}}, r)]}^{T},

(5a)

g_{d} (r) = {[G^{d} (r_{1}, r), G^{d} (r_{2}, r), \dots, G^{d} (r_{N_{d}}, r)]}^{T},

(5b)

where the superscript T denotes transpose; and N_s, N_d and M are the numbers of sources, detectors and targets, respectively. It is assumed the number of targets is less than the number of sources and detectors,

M < min (N_{d}, N_{s})

. It also holds that

K^{T} = {K_{j i}}

describes light propagation from the positions of detectors through the medium and targets to sources.

For a homogeneous background medium, the rank R of matrix K, is equal to the dimension of the source array vector space

G_{s}

spanned by

g_{s} (r_{m})

, and also equal to the dimension of the detector array vector space

G_{d}

spanned by

g_{d} (r_{m})

, where

G_{s} \subseteq C^{N_{s}}

and

G_{d} \subseteq C^{N_{d}}

. For absorptive targets, R is equal to the number of targets M.

Similar forms of the response matrix and GFVs can be obtained for scattering targets. As the dot product in the second term of Eq. (1) implies, each scattering target is represented by three components coexisting at one location. The elements of the K matrix for L scattering target may be written as

\begin{matrix} K_{i j} = \sum_{l = 1}^{L} τ_{l} \nabla_{r} G^{d} (r_{i}, X_{l}) \cdot \nabla_{r} G^{s} (X_{l}, r_{j}) \\ = \sum_{l = 1}^{L} τ_{l} \sum_{α = {x, y, z}} \partial_{α} G^{d} (r_{i}, X_{l}) \partial_{α} G^{s} (X_{l}, r_{j}), \end{matrix}

(6)

where

τ_{l} = δ D (X_{l}) c δ V_{l}

is the optical scattering strength of the l^th target. The K matrix for scattering targets can be written in a manner similar to that for absorptive targets:

K = \sum_{l = 1}^{L} \sum_{α = {x, y, z}} \partial_{α} g_{d} (X_{l}) τ_{l} \partial_{α} g_{s}^{T} (X_{l}) .

(7)

The Green’s function for a slab geometry is [16

M. Xu, M. Alrubaiee, S. K. Gayen, and R. R. Alfano, “Three-dimensional localization and optical imaging of objects in turbid media with independent component analysis,” Appl. Opt. 44(10), 1889–1897 (2005). [CrossRef] [PubMed]

,47

S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef] [PubMed]

]

G (r, r') = G (r', r) = \frac{1}{4 π D} \sum_{k = - \infty}^{\infty} (\frac{e^{- κ r_{k}^{+}}}{r_{k}^{+}} - \frac{e^{- κ r_{k}^{-}}}{r_{k}^{-}}),

(8a)

r_{k}^{\pm} = {[{(x - x')}^{2} + {(y - y')}^{2} + {(z \mp z' + 2 k d)}^{2}]}^{1 / 2},

(8b)

where

κ = {[(μ_{a} - i ω / c) / D]}^{1 / 2}

in frequency domain with amplitude modulation frequency ω, and

k = 0, \pm 1, \pm 2, \dots

. The extrapolated boundaries of the slab are located at

z = 0

and

z = d = L + 2 z_{e}

, respectively, where L is the physical thickness of the slab and the extrapolation length z_e is determined from the boundary condition of the slab [48

R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11(10), 2727–2741 (1994). [CrossRef] [PubMed]

,49

M. Lax, V. Nayaramamurti, and R. C. Fulton, “Classical diffusion photon transport in a slab,” in Laser Optics of Condensed Matter, J. L. Birman, H. Z. Cummins, and A. A. Kaplyanskii, eds. (Plenum, New York, 1987), pp. 229–237.

Under ideal conditions, when all three scattering components of each of the L scattering targets are well-resolved, the rank of K contributed by L scattering targets is 3L. In practice, four components (one for absorption and three for scattering) are calculated for each target, since the targets may have both scattering and absorptive characteristics, or the exact nature may not be known a priori. The dominant characteristic is used to label the target as absorptive or scattering in nature.

2.2 Point Spread Functions

If light emitted by a source of unit power at target position X propagates in the sample medium, the signal measured by the detector array at the sample boundary is

G^{d} (r_{i}, X)

. The signal is then “time-reversed” and back-propagated with the Green’s function of the background medium. TR operation is phase conjugation in Fourier domain [28

M. Fink, “Time reversal mirrors,” J. Phys. D Appl. Phys. 26(9), 1333–1350 (1993). [CrossRef]

,50

C. Prada, S. Manneville, D. Spoliansky, and M. Fink, “Decomposition of the time reversal operator: detection and selective focusing on two scatterers,” J. Acoust. Soc. Am. 99(4), 2067–2076 (1996). [CrossRef]

]. So the signal evaluated at r is [34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

]

\begin{matrix} H_{d} (r, X) = \sum_{i = 1}^{N_{d}} G^{d} (r, r_{i}) G^{d}^{*} (r_{i}, X) = g_{d}^{T} (r) g_{d}^{*} (X) \\ = g_{d}^{†} (X) g_{d} (r) = 〈 g_{d} (X), g_{d} (r) 〉, \end{matrix}

(9)

where * denotes phase conjugate, † denotes adjoint, and

< \cdot >

denotes inner product.

H_{d} (r, X)

is the detector array point spread function (PSF). A source array PSF can be similarly formed as

H_{s} (r, X) = g_{s} {(X)}^{†} g_{s} (r) = 〈 g_{s} (X), g_{s} (r) 〉 .

(10)

Due to the time reversal assumption,

H_{d} (r, X)

peaks at

r = X

, so it can be considered as an image of the source at X formed by the TR detector array. PSF vanishes when r is far away from X. A similar interpretation can be used for

H_{s} (r, X)

2.3 Time reversal and MUSIC

The TR matrix may be constructed to represent light propagation from sources to detectors and back denoted by T_SDDS, or to represent light propagation from detector positions to source positions and back denoted by T_DSSD, a consequence of the reciprocity of light propagation [29

C. Prada, L. Thomas, and M. Fink, “The iterative time reversal process: analysis of the convergence,” J. Acoust. Soc. Am. 97(1), 62–71 (1995). [CrossRef]

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,38

,50

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

]. For frequency-domain data,

T_{S D D S} = K^{†} K

, and

T_{D S S D} = {(K^{T})}^{†} K^{T} = K^{*} K^{T}

, where response data matrix K is formed using modulated intensities, instead of the field with phase information used in the conventional TR. For CW measurements,

T_{S D D S} = K^{T} K

, and

T_{D S S D} = K K^{T}

(K is real and only includes intensity values).

Since T_SDDS and T_DSSD are Hermitian (

T^{†} = T

), they have complete sets of orthonormal eigenvectors v_j (

j = 1, \dots, N_{s}

) and u_i (

i = 1, \dots, N_{d}

), with a common set of non-negative real eigenvalues. For

M < \min (N_{s}, N_{d})

absorptive targets without the presence of noise, the rank of T_SDDS and T_DSSD is M. The eigenvalues

X_{m}

, when

j = 1, \dots, M

, and

λ_{j} \approx 0

, when

j = M + 1, \dots, N_{s}

for T_SDDS and

j = M + 1, \dots, N_{d}

for T_DSSD. The eigen system

{v_{j}, u_{j}, λ_{j} > 0}

j = 1, \dots, M

, is related to the targets. The TR matrix T_SDDS can be written as [22

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

]

H_{d} (X_{m}, X_{m'}) \approx 0

(11)

Subsequent formalism may be different depending on whether the targets are “well resolved” or “poorly resolved.”

2.3.1 Well-resolved targets

If the

m^{t h}

and

m'^{t h}

targets (

m \neq m'

) are well resolved, defined by the conditions:

H_{s} (X_{m}, X_{m'}) \approx 0

and

H_{d} (X_{m}, X_{m'}) \approx 0

, i.e. the GFVs at

X_{m}

and

X_{m'}

are orthogonal,

< g_{d} (X_{m}), g_{d} (X_{m'}) > = H_{d} (X_{m'}, X_{m}) = {‖ g_{d} (X_{m}) ‖}^{2} δ_{m m'}

. So we have

T_{S D D S} = \sum_{m = 1}^{M} {| τ_{m} |}^{2} {‖ g_{d} (X_{m}) ‖}^{2} g_{s}^{*} (X_{m}) g_{s}^{T} (X_{m}),

(12)

where

‖ \cdot ‖

denotes L2 norm [52

N. Bourbaki, Topological Vector Spaces (Springer, 1987).

]. The eigenvectors of T_SDDS are proportional to the phase conjugate of the GFVs associated with the M targets [22

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

], i.e.

T_{S D D S} g_{s}^{*} (X_{m}) = {| τ_{m} |}^{2} {‖ g_{d} (X_{m}) ‖}^{2} {‖ g_{s} (X_{m}) ‖}^{2} g_{s}^{*} (X_{m}) .

(13)

The eigenvectors are

v_{j} = \frac{g_{s}^{*} (X_{j})}{‖ g_{s} (X_{j}) ‖},

(14)

with eigenvalues

λ_{j} = {| τ_{j} |}^{2} {‖ g_{d} (X_{j}) ‖}^{2} {‖ g_{s} (X_{j}) ‖}^{2}

j = 1, \dots, M

. Thus T_SDDS is a projection operator that projects a vector onto the conjugate of the source array vector space

G_{s}

. The j^th non-zero eigenvalue λ_j is directly related to the optical strength τ_j of the j^th target. Similar equations can be derived for T_DSSD, which is a projection operator for the conjugate of the detector array vector space

G_{d}

. The eigenvectors of T_DSSD are

u_{j} = \frac{g_{d}^{*} (X_{j})}{‖ g_{d} (X_{j}) ‖},

(15)

j = 1, \dots, M

, with the same eigenvalues as T_SDDS.

Therefore, for well-resolved targets, the target locations can be determined by the inner product [22

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

]

\begin{matrix} ψ_{j}^{s} = 〈 v_{j}^{*}, g_{s} (X_{p}) 〉 = v_{j}^{T} g_{s} (X_{p}) = \frac{g_{s}^{†} (X_{j})}{‖ g_{s} (X_{j}) ‖} g_{s} (X_{p}) \\ = \frac{1}{‖ g_{s} (X_{j}) ‖} H_{s} (X_{p}, X_{j}), \end{matrix}

(16a)

\begin{matrix} ψ_{j}^{d} = 〈 u_{j}^{*}, g_{d} (X_{p}) 〉 = u_{j}^{T} g_{d} (X_{p}) = \frac{g_{d}^{†} (X_{j})}{‖ g_{d} (X_{j}) ‖} g_{d} (X_{p}) \\ = \frac{1}{‖ g_{d} (X_{j}) ‖} H_{d} (X_{p}, X_{j}), \end{matrix}

(16b)

where X_p is a test target position, which is the position of any voxel in the sample space.

ψ_{j}^{s}

and

ψ_{j}^{d}

peak when X_p is the position of the j^th target. In the classical TR imaging [25

C. Prada, F. Wu, and M. Fink, “The iterative time reversal mirror: a solution to self-focusing in the pulse echo mode,” J. Acoust. Soc. Am. 90(2), 1119–1129 (1991). [CrossRef]

,29

C. Prada, L. Thomas, and M. Fink, “The iterative time reversal process: analysis of the convergence,” J. Acoust. Soc. Am. 97(1), 62–71 (1995). [CrossRef]

,38

,50

], for ideally resolved targets, each eigenvector of the TR operator can be used to focus on one particular target. Here ψ_j^s and ψ_j^d represent focusing of signals from the source and detector planes on to the target position, respectively. Use of the eigenvectors v_j and u_j,

j = 1, \dots, M

, ensures that the j^th target is sorted out. When T_SDDS and source array vector space (T_DSSD and detector array vector space) are used, we call the scheme SDDS (DSSD). Both source and detector arrays can be considered simultaneously to locate the target by calculating

ψ_{j} = ψ_{j}^{d} ψ_{j} {^{s}}^{†} = \frac{1}{‖ g_{s} (X_{j}) ‖ \cdot ‖ g_{d} (X_{j}) ‖} H_{d} (X_{p}, X_{j}) H_{s}^{*} (X_{j}, X_{p}),

(17)

j = 1, \dots, M

, which is computationally equivalent to a process that light emitted from a virtual source of unit power at a test target position X_p, propagates to the TR source array and back to a true target position X_j; then it is re-emitted and further propagates to the TR detector array and back to the original position X_p. ψ_j peaks when the test target position X_p coincides with the true target position X_j associated with the j^th eigenvector.

2.3.2 Poorly-resolved Targets and MUSIC

When the targets are too close to each other or the sources and/or detectors are significantly sparse, the targets are considered to be poorly resolved and the GFVs at X_m and X_m_' are not orthogonal. In such cases, the eigenvectors v_j and u_j do not correspond one-to-one with the GFVs associated with target positions X_m (

j, m = 1, \dots, M

). The image resolution degrades because of contributions from multiple targets. To solve this problem, the subspace-based method, MUSIC was implemented with TR [34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

]. MUSIC algorithm is based on the idea that although the vectors characterizing the targets are no longer orthogonal with each other, they are all located in the signal subspace, which is orthogonal to the noise subspace.

The orthonormal sets

{v_{j}^{*}}

(

j = 1, \dots, N_{s}

) and

{u_{j}^{*}}

(

j = 1, \dots, N_{d}

) span the spaces

C^{N_{s}}

and

C^{N_{d}}

associated with the source and detector arrays, respectively. While

{v_{j}^{*}}

and

{u_{j}^{*}}

, with

λ_{j} > 0

, form the signal subspaces on the source and detector arrays,

S^{s} = {v_{j}^{*}}

and

S^{d} = {u_{j}^{*}}

(

j = 1, \dots, M

), respectively;

{v_{j}^{*}}

and

{u_{j}^{*}}

, with

λ_{j} \approx 0

, form the noise subspaces,

N^{s} = {v_{j}^{*}}

(

j = M + 1, \dots, N_{s}

) and

N^{d} = {u_{j}^{*}}

(

j = M + 1, \dots, N_{d}

), respectively. Thus

C^{N_{s}} = S^{s} \oplus N^{s}

and

C^{N_{d}} = S^{d} \oplus N^{d}

[36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

]. Since the dimensions of the signal subspaces

S^{s}

and

S^{d}

and of the GFV spaces

G_{s}

and

G_{d}

are all equal to M,

G_{s} \equiv S^{s}

and

G_{d} \equiv S^{d}

[51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

]. The GFVs,

g_{s} (X_{m})

and

g_{d} (X_{m})

m = 1, \dots, M

, are linear combinations of v_j* and u_j*,

j = 1, \dots, M

, respectively. Therefore,

g_{s} (X_{m}) \in S^{s}

and

g_{d} (X_{m}) \in S^{d}

m = 1, \dots, M

, associated with m^th target are orthogonal to

v_{j}^{*} \in N^{s}

(

j = M + 1, \dots, N_{s}

) and

u_{j}^{*} \in N^{d}

(

j = M + 1, \dots, N_{d}

), respectively:

〈 v_{j}^{*}, g_{s} (X_{m}) 〉 = v_{j}^{T} g_{s} (X_{m}) \approx 0, j = M + 1, \dots, N_{s},

(18a)

〈 u_{j}^{*}, g_{d} (X_{m}) 〉 = u_{j}^{T} g_{d} (X_{m}) \approx 0, j = M + 1, \dots, N_{d} .

(18b)

The locations of targets can be determined by calculating the following squared sum of inner products:

Q_{s} (X_{p}) = \sum_{j = M + 1}^{N_{s}} {| v_{j}^{T} g_{s} (X_{p}) |}^{2},

(19a)

Q_{d} (X_{p}) = \sum_{j = M + 1}^{N_{d}} {| u_{j}^{T} g_{d} (X_{p}) |}^{2} .

(19b)

Q_{s} (X_{p})

and

Q_{d} (X_{p})

vanish when the test target position X_p is a true target position. Similar to Eq. (17),

Q = Q_{s} Q_{d}

can be calculated with both source and detector arrays considered simultaneously. An alternative approach to accentuate a target position is to plot a pseudo spectrum defined as

P_{s} (X_{p}) = {‖ g_{s} (X_{p}) ‖}^{2} / | Q_{s} (X_{p}) |

(20a)

associated with the source array, or

P_{d} (X_{p}) = {‖ g_{d} (X_{p}) ‖}^{2} / | Q_{d} (X_{p}) |

(20b)

associated with the detector array, or

P (X_{p}) = P_{s} (X_{p}) P_{d} (X_{p})

(20c)

associated with both source and detector arrays [22

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

], where

{‖ g_{s} (X_{p}) ‖}^{2}

and

{‖ g_{d} (X_{p}) ‖}^{2}

are used for normalization. The poles of the pseudo spectrum correspond to target locations. These MUSIC pseudo spectra can also be used to locate well-resolved targets.

Since the dimension of the signal subspace is generally much smaller than that of the noise subspace, it is preferred that in Eq. (19) and Eq. (20), the signal subspace is used rather than the noise subspace for ease of computation. Using the properties of the projection operators associated with the source and detector arrays [22

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

Q_{s} (X_{p})

and

Q_{d} (X_{p})

can be calculated as

Q_{s} (X_{p}) = {‖ g_{s} (X_{p}) ‖}^{2} - \sum_{j = 1}^{M} {| v_{j}^{T} g_{s} (X_{p}) |}^{2},

(21a)

Q_{d} (X_{p}) = {‖ g_{d} (X_{p}) ‖}^{2} - \sum_{j = 1}^{M} {| u_{j}^{T} g_{d} (X_{p}) |}^{2} .

(21b)

When the targets are embedded in a non-uniform medium, or when there is significant noise present, the noise or false targets contribute significantly to the eigenvalues. The near-zero and non-zero eigenvalues are not as well separated as when there are no noise. In this case, the rank of the TR matrix is larger than the number of targets M. The TR matrix may even be full rank. However, as long as M is less than

\min (N_{s}, N_{d})

and eigenvalues contributed by the noise and false targets are smaller than those contributed by the real targets with a threshold

ϵ

, the target positions can be obtained using a pseudo spectrum [36

S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]

,51

A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

] associated with the source array,

P_{s} (X_{p}) = {‖ g_{s} (X_{p}) ‖}^{2} / | Q_{s} {(X_{p})}_{λ_{j} \leq ϵ} |,

(22)

where

Q_{s} {(X_{p})}_{λ_{j} \leq ϵ} = {‖ g_{s} (X_{p}) ‖}^{2} - \sum_{λ_{j} > ϵ} {| v_{j}^{T} g_{s} (X_{p}) |}^{2}

. Pseudo spectra associated with the detector array or with both source and detector arrays can also be obtained similarly. In practice, the threshold is selected to separate the signal and noise subspaces using a method similar to L-curve regularization [39

P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. 34(4), 561–580 (1992). [CrossRef]

When scattering targets are concerned, the GFVs

\partial_{α} g

(

α = x, y, z

), associated with the test target position X_p will be used to calculate the pseudo spectrum. For a target with both absorption and scattering properties at the wavelength of probing light, one GFV corresponding to absorption constructed as g and three GFVs corresponding to scattering target constructed with

\partial_{α} g

, (

α = x, y, z

), are used to calculate the pseudo spectrum over every voxel. Ideally, for an absorptive and scattering target four pseudo-values will be obtained for every target position. If the dominant value corresponds to the absorptive (any of the scattering) GFV the target will be identified as absorptive (scattering) in nature.

3. Algorithm

Implementation of TROT to locate targets embedded in a highly scattering turbid medium involves the steps outlined below. For simplicity, the sizes of source array and detector array are assumed to be the same, i.e.,

N_{d} = N_{s} = N

(a) A response matrix K with size N × N is constructed using experimental data (or estimated data in simulation). Data consist of the perturbations in the light intensity distribution due to the targets, $Δ ϕ = ϕ - ϕ_{0}$ , where $ϕ$ is the light intensity distribution measured on the sample boundary with targets embedded in the scattering medium and $ϕ_{0}$ is ideally the light intensity distribution without the targets. In practice, $ϕ_{0}$ is approximated by an “average” over all the multi-source measurements, while in simulation it can be estimated without such approximation.
(b) A detector array TR matrix, $T_{D S S D} = K K^{T}$ with size N × N for CW measurements is constructed. All the eigenvalues and the eigenvectors of the T_DSSD matrix are computed. The eigenvectors are orthogonal to each other. It is to be noted that in this procedure we only deal with a matrix of dimension N, not a matrix of dimension of N × N as done in traditional inverse procedures.
(c) The non-zero eigenvalues of T_DSSD belonging to the signal subspace are separated from the near-zero eigenvalues belonging to the noise subspace using the L-curve method [5
W. Cai, S. K. Gayen, M. Xu, M. Zevallos, M. Alrubaiee, M. Lax, and R. R. Alfano, “Optical tomographic image reconstruction from ultrafast time-sliced transmission measurements,” Appl. Opt. 38(19), 4237–4246 (1999). [CrossRef] [PubMed]
,39
P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. 34(4), 561–580 (1992). [CrossRef]
,40
E. A. Marengo, F. K. Gruber, and F. Simonetti, “Time-reversal MUSIC imaging of extended targets,” IEEE Trans. Image Process. 16(8), 1967–1984 (2007). [CrossRef] [PubMed]
].
(d) MUSIC approach [34
A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.
,36
S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]
,51
A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]
] is next used to determine the locations of the targets as follows. (i) The 3-D medium is divided into a certain number of voxels. A detector array GFV, $g_{d} (X_{p})$ , associated with an absorptive test target position X_p at p^th voxel is calculated using Diffusion Approximation of RTE. Other proper forward models could be used as well. In order to check if $g_{d} (X_{p})$ is located in the signal subspace or in the noise subspace, a pseudo spectrum associated with the detector array is computed using Eq. (20b), where M is the dimension of the signal subspace found in step (c). If $g_{d} (X_{p})$ is located in the signal subspace, the corresponding pseudo value $P (X_{p})$ in Eq. (20b) will become a maximum. (ii) Pseudo spectra are also calculated using the other three GFVs, $\partial_{α} g_{d} (X_{p})$ , ( $α = x, y, z$ ) for scattering property. (iii) All pseudo values are put together and sorted in a descending order. Since the leading pseudo values at X_p are associated with targets and specific GFVs, the positions of the embedded targets and their nature (absorptive or scattering) are determined. The pseudo spectrum in the whole sample space can be used to plot pseudo tomographic images.

In this approach, only a single run is needed for calculating the pseudo spectrum and no iterative procedure is involved, which makes it faster and computationally less intensive than the traditional IIR approaches. Similar procedure can be used for application of TROT when

N_{d} \neq N_{s}

. The pseudo spectrum associated with either the source array, or the detector array, or both source and detector arrays, as outlined in Eq. (20) can be used to obtain target positions.

It is instructive to compare the computational complexity of TROT formalism with that of typical iterative methods. For a typical iterative method, an equation b = Wx is solved to find the inhomogeneities (targets), where W is a weight matrix with size N_dN_s × N, N is the number of voxels, b is an N_dN_s × 1 vector describing the perturbation in the detected light intensity due to the presence of inhomogeneities, and x is the perturbation or variation in the optical properties from the background values with dimension of N × 1. The computational complexity is typically O(N_dN_sN²) for a single iteration. The computational complexity of TROT is much smaller than that for even one iteration of an iterative method. For the SDDS scheme, the complexity for TROT is O(N_dN_s²) if N_dN_s > NN_k, and O(N_sNN_k) otherwise, where N_k is the dimension of the signal subspace. In the DSSD scheme, the complexity is O(N_sN_d²) if N_dN_s > NN_k, and O(N_dNN_k) otherwise. TROT does not involve any iteration.

In the following sections, TROT will be tested using simulation and experimental data.

4. Testing TROT with simulated data

To test the efficacy of the TROT approach, we first consider a rather challenging task of detecting and locating six targets embedded in a simulated sample which is a 40-mm thick uniform scattering slab. Its absorption and diffusion coefficients are µ_a = 1/300 mm⁻¹ and D = 1/3 mm, respectively. The incident CW beam was step-scanned in an x-y array of 41 × 41 grid points with a step size of 2 mm on the input plane covering an 80 mm × 80 mm area. Light transmitted from the opposite side (output plane) was recorded at 41 × 41 grid points covering the same area. No random noise was added.

The six (M = 6) point-like absorptive targets, with absorption coefficient difference of ∆µ_a = 0.01 mm⁻¹ from the background, were placed at A (24 mm, 26 mm, 9 mm), B (38 mm, 38 mm, 15 mm), C (38 mm, 38 mm, 21 mm), D (40 mm, 38 mm, 21 mm), E (44 mm, 42 mm, 21 mm) and F (30 mm, 30 mm, 31 mm), respectively. The origin (0 mm, 0 mm, 0 mm) was located at the upper-left corner of the input boundary (source plane) of the sample. The medium was divided into 40 × 40 × 20 voxels, with each voxel of size 2 mm × 2 mm × 2 mm. As can be seen from the assigned coordinates, targets C and D are located at two adjacent voxels, and are close to target E, and these three targets are located in the same z layer. Consequently, targets C and D are expected to be very difficult to resolve, and hard to distinguish from target E. Target B and C have the same lateral position x and y, and different depths. Target A is close to the source plane, while F is close to the detector plane.

Using the Diffusion Approximation of the RTE as the model for light propagation in slab geometry, signals arising from light propagation from the source array to the detector array through medium with and without the targets were calculated. The difference between the two sets, which is the perturbation due to the targets, was used as the “simulated data”. The size of the K matrix is N × N = 1681 × 1681. The TR matrix T = KK^T was constructed. Then, 1681 eigenvalues and 1681 eigenvectors of T were found.

The first seven (7) computed eigenvalues in a descending order of magnitude are listed in the first column of Table 1 . The leading twenty eigenvalues are plotted in Fig. 1(a) on a logarithmic scale. The first six (6) eigenvalues are at least 10 orders-of-magnitude higher than the 7-th and other smaller eigenvalues. Hence, the dimension of the signal subspace and the number of targets are determined to be six. The pseudo spectrum (consisting of 40 × 40 × 20 × 4 elements) was calculated using the M eigenvectors in the signal subspace. The values of elements in the pseudo spectrum were sorted in a descending order. The seven leading pseudo values are listed in Table 1 with the corresponding positions of voxels. The six peaks are found to be associated with the GFVs for absorptive targets. Namely, the corresponding six targets are identified to be absorptive targets.

Table 1 Eigenvalues, pseudo spectrum and the corresponding positions

Leading Eigenvalues	Poles of Pseudo Spectrum	Retrieved Position  (x, y, z) mm	Know Position  (x, y, z) mm
2.6697E-010	1.5911E + 015	(44, 42, 21)	(44, 42, 21)
1.1722E-011	8.6376E + 014	(38, 38, 15)	(38, 38, 15)
4.0081E-013	7.9559E + 014	(38, 38, 21)	(38, 38, 21)
3.6676E-014	7.2328E + 014	(40, 38, 21)	(40, 38, 21)
5.2629E-016	6.3010E + 014	(24, 26, 9)	(24, 26, 9)
6.4837E-017	2.1159E + 014	(30, 30, 31)	(30, 30, 31)
2.8337E-039	2.4353E + 005	(38, 38, 19)
…	…

Fig. 1 (a) A plot of first twenty (20) eigenvalues on logarithmic scale. (b) 2-D slice of the pseudo spectrum on z = 21 mm plane showing the location of the three difficult targets described in the text. Similar 2-D slices were also obtained for z = 9-mm, 15-mm, and 31-mm planes (not shown).

All six large pseudo-values are located at the exact known target locations and their values are approximately 9 orders-of-magnitude larger than those at their neighborhood locations. A 2-D slice of the pseudo spectrum on z = 21 mm plane is shown in Fig. 1(b), showing the locations of the three difficult targets.

With the highly encouraging result from simulation even for a considerably challenging task, we proceeded to test the approach for the realistic situation of detecting and locating targets from experimental data.

5. Testing TROT using Experimental Data

5.1 Experimental materials and methods

Three different experiments with three different samples were carried out to test the efficacy of the TROT approach in detecting and locating targets in a turbid medium. All three samples used a 250 mm × 250 mm × 60 mm transparent plastic container filled with Intralipid-20% suspension in water as the background medium. The concentration of Intralipid-20% was adjusted to provide an estimated [53

H. J. van Staveren, C. J. M. Moes, J. van Marie, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. 30(31), 4507–4514 (1991). [CrossRef] [PubMed]

,54

C. Bordier, C. Andraud, E. Charron, J. Lafait, M. Anastasiadou, and A. Martino, “Illustration of a bimodal system in Intralipid 20% by polarized light scattering: experiments and modelling,” Appl. Phys., A Mater. Sci. Process. 94(2), 347–355 (2009). [CrossRef]

] absorption coefficient µ_a ~0.003 mm⁻¹ at 790 nm, and a transport mean free path l_t ~1 mm, which were similar to the average values of those parameters for human breast tissue, while the cell thickness of 60 mm was comparable to thickness of a typical compressed breast.

In the first experiment, the depth (position along z-axis) of an absorptive target was varied to explore how the accuracy of position estimate depended on depth. The target was a glass sphere of diameter ~9 mm filled with ink dissolved in Intralipid-20% suspension in water (µ_s was adjusted to be the same as that of the background medium, while µ_a ≈0.013 mm⁻¹ was about 3 times higher than that of background medium).

In the second experiment, the separation between two absorptive targets was varied to test how close those could be and yet be resolved as separate objects. Both the targets were similar to the target in the first experiment.

In the third experiment, the depth of a scattering target was varied to explore the efficacy of TROT in locating and characterizing a scattering target. The target was a glass sphere of diameter 10 mm filled with Intralipid-20% suspension in water to provide a transport mean free path l_t of 0.25 mm, and a scattering coefficient µ_s ≈11 mm⁻¹.

A multi-source interrogation and multi-detector signal acquisition scheme, shown in Fig. 2 , was used to acquire data. A 100-mW 790-nm diode laser beam was used to illuminate the samples. A 1024 × 1024 pixels charge coupled device (CCD) camera equipped with a 60-mm focal-length camera lens was used on the opposite side of the sample to detect the transmitted light on the boundaries of the slab samples (detector plane). The pixel size was 24 µm. The multi-source illumination scheme was realized by scanning the sample across the laser beam in a two-dimensional x-y array of grid points using a computer-controlled translation stage. The first and third samples were scanned across the laser beam in an array of 9 × 9 grid points, and the second sample was scanned in an array of 15 × 11 grid points, with a step size of 5 mm in all cases. The scanning and data acquisition processes were controlled by a personal computer (PC). Raw transillumination images of the sample were recorded by the PC for each scan position, and stored for subsequent analysis. A typical image, which is a 2-D intensity distribution, is shown in the right side of Fig. 2.

Fig. 2 A schematic diagram of the experimental arrangement for imaging objects embedded in a turbid medium. (Key: CCD = charge coupled device, PC = personal computer) Inset (below) shows the 2-D array in the input plane that was scanned across the incident laser beam, and inset (right) shows a typical raw image.

While we have scanned the sample and kept the source fixed in the experiments reported here, a more clinically relevant approach would be to scan the source and fix the sample. In the experimental arrangement, the source scanning may be accomplished by: (a) delivering the beam using an optical fiber, and translating the delivery end of the fiber in an x-y array using a computer-controlled translation stage; or (b) raster scanning the laser beam using two orthogonal (x-y) galvanometers. The main change in the processing of data would involve alignment of the images so that laser beam positions are overlapped before averaging to generate the background image.

5.2 Data Processing and Analysis

From each image, a region of interest was cropped out and then every 5 × 5 pixels in the cropped image were binned to one pixel to enhance the signal-to-noise ratio. The background image was generated by taking an average of the original images for all scan positions, which is a reasonable approximation since for most of the scan positions the target(s) is (are) not along the direction of the incident beam. Then the background image was also cropped and binned corresponding to the region of interest for each scan position. Perturbation in the light intensity distribution

Δ ϕ

due to targets in each image was found by subtracting the background image from each individual image. The response matrix was constructed using the light intensity perturbations,

- Δ ϕ

. The TR matrix was generated by multiplying the response matrix by its transpose for our continuous-wave (CW) probing scheme. The eigenvalue equation was solved and the signal subspace was selected and separated from the noise subspace. MUSIC was then used to calculate the pseudo spectrum for all voxels in the 3-D space of the sample. For each voxel, four pseudo values, one for absorption and three for scattering as described in Algorithm step (d), were calculated. The voxel size was 0.77 mm × 0.77 mm × 1 mm. By sorting the pseudo spectrum in a descending order, the target(s) were located.

The voxel size to be used in reconstruction and its relation to the target size is an important consideration. In general, smaller voxels provide reconstruction of higher resolution at the cost of increased computational time. Finer details of an extended target may be obtained using smaller voxels. Decreasing the voxel size indefinitely may not improve resolution because of the diffusive nature of light propagation in the turbid medium. However, the computation time increases dramatically, which has been observed by other researchers [55

S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16(7), 5048–5060 (2008). [CrossRef] [PubMed]

]. The optimal voxel size for a given reconstruction problem will depend on factors, such as, target size, experimental geometry, and noise level.

Since the signal used in image reconstruction is taken to be the difference between the image recorded for a scan position and the background image, estimation of the background image is an important issue. This is a common problem for every diffuse optical imaging modality using perturbation method, and needs further elaboration. We accumulated data in the transillumination slab geometry, and generated the background image by averaging images for all scan positions after proper alignment with respect to the incident source position. This averaging method for generating background image worked well for small targets that we used in our experiments, as the ratio of sample volume to target volume was quite high (~500:1). This volume ratio for breast and a tumor in early stages of growth will also be substantially high for the averaging method to be applicable. Assuming a scenario where the volume ratio is substantially smaller than in above examples, a modified approach would be to select recorded images which were minimally affected by embedded targets for averaging [56

M. Xu, Y. Pu, and W. Wang, “Clean image synthesis and target numerical marching for optical imaging with backscattering light,” Biomed. Opt. Express 2(4), 850–857 (2011). [CrossRef] [PubMed]

]. As long as the targets only occupy a limited volume within the host medium, a clean background image can be generated in this fashion. It should also be noted that while estimation of target optical properties, such as absorption coefficient and scattering coefficient, are sensitively dependent on background image estimation, estimation of target positions are not so sensitive.

Several alternative ways of generating background image are suggested in the literature. One experimental approach is to record image using a phantom that has the same average optical properties as the sample, such as human breast [57

T. Nielsen, B. Brendel, R. Ziegler, M. van Beek, F. Uhlemann, C. Bontus, and T. Koehler, “Linear image reconstruction for a diffuse optical mammography system in a noncompressed geometry using scattering fluid,” Appl. Opt. 48(10), D1–D13 (2009). [CrossRef] [PubMed]

]. Along the same line, image of the healthy contralateral breast taken under the same experimental conditions as that of the suspect breast may be used as background image for breast imaging [58

Y. Ardeshirpour, N. Biswal, A. Aguirre, and Q. Zhu, “Artifact reduction method in ultrasound-guided diffuse optical tomography using exogenous contrast agents,” J. Biomed. Opt. 16(4), 046015 (2011). [CrossRef] [PubMed]

]. Some authors have suggested acquiring data at a wavelength for which the target(s) and the background have identical optical properties for assessing the background, e.g., measurement using 805-nm light for which hemoglobin and oxyhemoglobin have the same absorption coefficient may serve as background to image hemoglobin oxygenation [59

B. W. Pogue, M. S. Patterson, H. Jiang, and K. D. Paulsen, “Initial assessment of a simple system for frequency domain diffuse optical tomography,” Phys. Med. Biol. 40(10), 1709–1729 (1995). [CrossRef] [PubMed]

]. Still another approach is to compute the background using an appropriate forward model [18

]. Any of these approaches may be employed for generating the background image for use with the TROT formalism presented here.

The geometries commonly used in DOT include slab, cylindrical, hemispherical, and semi-infinite; and different source-detector combinations have been used to record images in these geometries. As long as multiple source-detector combinations provide multiple angular views of the sample the TROT formalism can be adapted to obtain target location for these geometries. TR imaging and TR-MUSIC was originally developed for reflection (backscattering) geometry that used coincident transceiver array to detect the return signal [28

M. Fink, “Time reversal mirrors,” J. Phys. D Appl. Phys. 26(9), 1333–1350 (1993). [CrossRef]

–30

A. J. Devaney, E. A. Marengo, and F. K. Gruber, “Time-reversal-based imaging and inverse scattering of multiply scattering point targets,” J. Acoust. Soc. Am. 118(5), 3129–3138 (2005). [CrossRef]

,34

A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.

,35

]. With requisite modification in the experimental arrangement TROT would be suited for use in the reflection geometry.

5.3 Results

5.3.1 Single absorptive target at different depths

In the first experiment, only one target was used, the lateral (x, y) position of the target was kept the same at (25.5 mm, 24.7 mm), while seven different depths (position along z-axis) of 15 mm, 20 mm, 25 mm, 30 mm, 35 mm, 40 mm and 45 mm were used. The eigenvalue spectrum plotted on logarithmic scale for the target at z = 30 mm is shown in Fig. 3 . Similar eigenvalue spectra were obtained for other cases.

Fig. 3 A semi-log plot of eigenvalue spectrum with first 40 leading eigenvalues for the target at z = 30 mm.

Both SDDS and DSSD pseudo spectra were calculated using Eq. (20). The target was identified as an absorptive target. In the DSSD pseudo spectrum, the absorptive pseudo value at the peak position is ~41 times of the scattering pseudo value associated with

\partial_{z} g_{d}

, and even larger than those associated with

\partial_{x} g_{d}

, and

\partial_{y} g_{d}

, as shown in Table 2 . Similarly in the SDDS scheme, the absorptive pseudo value at the peak position is ~33 times of the scattering pseudo value with

\partial_{z} g_{s}

, and much larger than the other two.

Table 2 Pseudo values associated with absorptive and scattering components at the peak position

Scheme with GFV (g)	Absorptive pseudo value	Scattering pseudo value
Scheme with GFV (g)	Absorptive pseudo value	$\partial_{x} g$	$\partial_{y} g$	$\partial_{z} g$
DSSD (g_d)	1305.0	1.0	1.0	31.7
SDDS (g_s)	2729.3	14.0	1.1	81.6

Three-dimensional tomographic images were generated using the whole absorption pseudo spectrum for all voxel positions in the sample. The left pane of Fig. 4(a) shows a tomographic image for the target at z = 30 mm. The spatial profiles in the x, y and z directions, shown in the right pane of Fig. 4(a) were used to assess the target location. Similar images were generated for other depths. The retrieved target positions are compared with known positions in Table 3 .

Fig. 4 Pseudo image of the target (left pane) and corresponding spatial intensity profiles (right pane) when the target is located at z = 30 mm: (a) experimental data; (b) simulation without any added noise; and (c) simulation with 20% Gaussian noise added. The pseudo values are calculated using Eq. (20).

Table 3 Positions of one target located at different depths

Known Positions x, y, z (mm)	DSSD Scheme		SDDS Scheme
Known Positions x, y, z (mm)	Retrieved Positions x, y, z (mm)	Error Δx, Δy, Δz (mm)	Retrieved Positions x, y, z (mm)	Error Δx, Δy, Δz (mm)
25.5, 24.7, 15	24.9, 24.4, 17.5	0.6, 0.3, 2.5	24.9, 25.2, 15.5	0.6, 0.5, 0.5
25.5, 24.7, 20	25.7, 24.4, 21.5	0.2, 0.3, 1.5	24.9, 25.2, 20.5	0.6, 0.5, 0.5
25.5, 24.7, 25	25.7, 24.4, 26.5	0.2, 0.3, 1.5	25.7, 24.4, 24.5	0.2, 0.3, 0.5
25.5, 24.7, 30	25.7, 24.4, 30.5	0.2, 0.3, 0.5	25.7, 25.2, 29.5	0.2, 0.5, 0.5
25.5, 24.7, 35	25.7, 25.2, 33.5	0.2, 0.5, 1.5	24.9, 24.4, 36.5	0.6, 0.3, 1.5
25.5, 24.7, 40	24.9, 25.2, 36.5	0.6, 0.5, 3.5	24.9, 25.9, 40.5	0.6, 1.2, 0.5
25.5, 24.7, 45	24.9, 25.2, 39.5	0.6, 0.5, 5.5	24.9, 26.7, 45.5	0.6, 2.0, 0.5

As is evident from Table 3, when DSSD scheme was used, the TROT-assessed lateral positions (x, y) were within 0.6 mm of the known values, which is an excellent agreement. The accuracy of the z-position was found to be optimal when the target was located in the middle plane of the sample, and deteriorated when the target was closer to the source plane or the detection plane. When using SDDS scheme, the TROT-assessed lateral positions were also within 0.6 mm of the known positions, except for z = 40 mm and 45 mm, where the error in y direction was 1.2 mm and 2 mm, respectively. However, remarkable improvement in the accuracy of the z-position estimation was observed, the error Δz being within 0.5 mm for all cases except for z = 35 mm, where the error was 1.5 mm. We ascribe the superior performance of the scheme using T_SDDS, to the much larger size of the detector array (1024 × 1024) than that of the source array (9 × 9) used in the experimental arrangement.

It should be noted that the choice of either DSSD or SDDS scheme depends on experimental parameters, such as, the number and density of sources and detectors, and does not depend on the characteristics of the background medium. When more detectors than sources are used and inter-detector spacing is small, SDDS would provide better resolution than DSSD, and vice versa. However, due to the diffusive nature of light propagation in the turbid medium, increasing the numbers and decreasing the spacing of the sources/detectors beyond a limit may not always improve the results.

While the target position could be obtained from the experimental data, it was observed that the difference between the smaller eigenvalues in the signal subspace and the larger eigenvalues in the noise subspace were not as pronounced as observed in the simulation in Section 4. To assess the effect of noise and to what extent noise may be present in the experimental data; we generated simulated data mimicking the experimental conditions, and added different noise levels. The lateral positions were (25.5 mm, 24.7 mm) and all seven z-positions (depth) of 15 mm, 20 mm, 25 mm, 30 mm, 35 mm, 40 mm, and 45 mm were tested. Typical pseudo images generated for z = 30 mm without and with 20% Gaussian multiplicative noise to compare with the experimental result are shown in Fig. 4(b) and Fig. 4(c), respectively. Simulated data provided the known position coordinates.

The simulated spatial profiles with zero added noise are much sharper than the profiles obtained from experimental data, or from simulated data with 20% added Gaussian noise. Broadening of spatial profile is an indication of the uncertainty in determination of position coordinates. Results from simulation show that uncertainty in position determination increases with added noise, and that experimental data behave in a way similar to simulated data with added noise.

5.3.2 Resolving two absorptive targets

In the second experiment using two targets the depth (z) and height (y) were kept same (z = 30 mm, y = 26.0 mm), while three different center-to-center separations, Δx of ~12.6 mm, 17.6 mm, and 27.6 mm, between them along the x-axis were considered. A cross-sectional pseudo image of the targets when separated by a center-to-center distance of 27.6 mm, generated using the pseudo spectrum is shown in the left pane of Fig. 5(a) . Figure 5(b) shows a similar image for the separation 12.6 mm (separation between nearest edges ~4 mm). A similar image for the separation 17.6 mm was also obtained (not shown in the figure). The profiles in the x, y and z directions through the right target are shown in the right panes of Fig. 5(a)–Fig. 5(c). These profiles were used to assess locations of the targets, and the separation between the two targets. In all cases, the targets were determined to be absorptive, because peaks occurred in the pseudo spectrum with the GFVs corresponding to absorption property.

Fig. 5 (a) (Experiment): TROT generated cross-section pseudo image when the targets are separated by 27.6 mm is shown in the left pane and pseudo-value profiles through the right target along x, y and z directions are shown in the right pane. (b) (Experiment): TROT generated cross-section pseudo image when the targets are separated by 12.6 mm is shown in the left pane and the corresponding spatial profiles through the right target along x, y and z directions are shown in the right pane. (c) (Simulation): TROT generated cross-section pseudo image when two targets are separated by 12.6 mm is shown in the left pane and the corresponding pseudo-value profiles are plot in the right pane. In simulation 10% Gaussian noise is added for comparison with the experimental results. P is pseudo value calculated using Eq. (20).

The known and retrieved positions from the experiments and separations Δx between the two targets appear in Table 4 . In all the cases, the two targets were resolved, even when their center-to-center separation was 12.6 mm apart, nearest sides separated by only ~4 mm. For all retrieved positions, the maximum error in the lateral positions is 3.0 mm, and the maximum error in the axial positions is 1.5 mm. The errors in the lateral positions increase as the targets get closer. We ascribe this increase in error in the lateral position to the crosstalk between the two targets, the peak due to one target being affected by the other. The shift in the peaks is also affected by noise. When the two targets are very close or significant noise is present, the two peaks merge, so that the two targets are not resolved. This behavior was confirmed in simulations.

Table 4 Positions of two targets separated with different distances

Known Separation [Δx (mm)]	Target #	Known Position [x, y, z (mm)]	Retrieved Position [x, y, z (mm)]	Error (mm)	Retrieved Separation [Δx (mm)]
12.6	1	27.6, 26.0, 30	30.3, 24.4, 28.5	2.7, 1.6, 1.5	6.9
12.6	2	40.2, 26.0, 30	37.2, 25.2, 29.5	3.0, 0.8, 0.5	6.9
17.6	1	25.1, 26.0, 30	26.4, 24.4, 28.5	1.3, 1.6, 1.5	14.6
17.6	2	42.7, 26.0, 30	41.0, 25.2, 29.5	1.7, 0.8, 0.5	14.6
27.6	1	20.1, 26.0, 30	19.5, 25.2, 29.5	0.6, 0.8, 0.5	27.6
27.6	2	47.7, 26.0, 30	47.1, 25.2, 30.5	0.6, 0.8, 0.5	27.6

The results were compared with simulated data using similar conditions. For the more challenging case of two targets located at z = 30 mm and separated by 12.6 mm, exact target locations were found when no noise was added. With 10% noise present, the positions of the two targets were found to be (39.0 mm, 24.8 mm, 29.0 mm) and (30.0 mm, 24.8 mm, 29.0 mm) with target separation 9.0 mm, compared to 12.6 mm (known) and 6.9 mm retrieved from experiment. The pseudo image and the corresponding profiles through the right target are shown in Fig. 5(c). Similar images were also obtained for the left target. The retrieved separation between the two targets in simulation with 10% noise was smaller than the actual separation. But the error was less than the experimental result. However, when 20% noise was added, the two peaks merged (not shown here).

The limits on the size of targets, separation between the targets, and the target-to-background contrast ratio that are needed to detect and locate the targets depend on noise level, experimental parameters (such as, number and concentration of sources and detectors), and ultimately on the diffuse nature of light propagation in the turbid medium. Coordinated experimental work and numerical modeling will be needed to assess those limits.

5.3.3 Single scattering target at different depths

The experiment involving the third sample is the same as the first one except that the target was scattering in nature. The scattering target was a 10-mm diameter glass sphere filled with Intralipid-20% suspension in water, whose concentration was adjusted to provide l_t = ~0.25 mm, µ_s = 11.3 mm⁻¹. The same scanning and data acquisition scheme was used as for the absorptive targets and the following z-positions of the target were used: 15 mm, 20 mm, 25 mm, 30 mm, 35 mm, 40 mm, and 45 mm. DSSD scheme was used to calculate the pseudo spectrum. A cross-section pseudo image and the corresponding spatial profiles are displayed in Fig. 6(a) when z = 30 mm. It is compared to the simulation results with 20% Gaussian noise (Fig. 6(b)). The lateral (x, y) spatial profiles of the pseudo image generated using simulated data are considerably wider, while the axial (z) spatial profile is narrower than those obtained using experimental data, and the peak values from the two cases are of the same order. The retrieved target positions are listed in Table 5 . SDDS scheme was also used and provided with similar results.

Fig. 6 Pseudo image of the target (left pane) and corresponding spatial intensity profiles (right pane) when the target is located at z = 30 mm: (a) experimental data; (b) simulation with 20% Gaussian noise added. P is pseudo value calculated using Eq. (20).

Table 5 Positions of one scattering target located at different depths

Known Positions [x, y, z (mm)]	Retrieved Positions [x, y, z (mm)]	Error [Δx, Δy, Δz (mm)]
25.7, 24.5, 15	24.9, 25.9, 18.5	0.8, 1.4, 3.5
25.7, 24.5, 20	27.2, 26.7, 20.5	1.5, 2.2, 0.5
25.7, 24.5, 25	25.7, 26.7, 23.5	0.0, 2.2, 1.5
25.7, 24.5, 30	24.9, 25.2, 32.5	0.8, 0.7, 2.5
25.7, 24.5, 35	24.9, 25.2, 36.5	0.8, 0.7, 1.5
25.7, 24.5, 40	24.9, 25.9, 41.5	0.8, 1.4, 1.5
25.7, 24.5, 45	24.9, 25.9, 45.5	0.8, 1.4, 0.5

In Fig. 5, and more prominently in Fig. 6, the image resolution seems better for experimental data than simulated data. Since the peak values and bandwidth of lines (the poles) in the pseudo spectrum depend strongly on the noise, this difference in image resolution is presumably due to lower noise level in the experiments than that used in simulations.

A comparison of experimental results for scattering and absorptive targets validate the common notion that it is more challenging to locate and image scattering targets than absorptive targets in a highly scattering medium. Also the lateral (x, y) positions are determined with higher accuracy than the axial (z) position. Overall the TROT-retrieved target positions are in good agreement with the known positions.

6. Discussion

The article presents the development of time reversal imaging approach with subspace classification, MUSIC in the optical domain. The results from experiment and simulation show that TROT is a faster and less computation intensive approach for detecting small targets in highly scattering turbid media and determining their locations in 3-D than other inverse image reconstruction techniques. While the dominant features in the pseudo spectrum are related to the square of the difference between the absorption (scattering) coefficient of the targets and that of the background, the approach does not directly determine these parameters. It is common for IIR approaches to estimate the optical properties of every voxel in the sample and identify target(s) from differences of these properties between the sample and the target(s), which is a considerably computation intensive undertaking. On the contrary, TROT identifies the targets as poles of the pseudo spectrum and focuses on determining their positions, which do not require as much computation time. Other IIR approaches involve iteration, while TROT is non-iterative. In TROT the data dimension is lower compared to other IIR approaches, which enables analysis and utilization of very large data sets. These two features together make TROT faster. Fast image reconstruction algorithms are of particular interest.

It was observed that lateral (x, y) positions are better determined than the depth (z). Also the spatial profile is more spread out along z compared to that along x, y. We ascribe this difference to fewer data along z-direction compared to those along x-y planes. Addition of another set of data with light incident and signal collected perpendicular to the z-direction is expected to further improve resolution in this dimension. Even without that addition, TROT determines the target position well.

While we have used slab geometry and CW illumination, the TROT approach may be used for other geometries (such as, cylindrical, and spherical), different types of illumination (e.g. frequency domain and pulsed) and different models for light propagation through the medium. Application and adaption of the TROT formalism to inhomogeneous media and extended targets may require careful consideration of several factors. In a non-uniform, inhomogeneous medium, structures other than the desired targets may appear as “false targets” and may interfere with identification of “real targets”. However, as long as the contributions to the signal by any false target is smaller than that made by real targets, TROT with MUSIC will be useful in detecting and locating targets, by choosing a proper threshold to separate the signal and noise subspaces. What is even more important, expected wavelength dependence of the target spectroscopic properties could be used to assess the difference between the real and false targets in experiments using multi-wavelength interrogation of the sample.

The TROT formalism presented in this article is particularly suited for point-like targets requiring fewer eigenvectors in the signal subspace to construct a pseudo spectrum. However, for extended finite-size targets, the formalism needs to be modified and much more eigenvectors may be needed to calculate the pseudo spectrum [40

E. A. Marengo, F. K. Gruber, and F. Simonetti, “Time-reversal MUSIC imaging of extended targets,” IEEE Trans. Image Process. 16(8), 1967–1984 (2007). [CrossRef] [PubMed]

,60

S. Hou, K. Solna, and H. Zhao, “Imaging of location and geometry for extended targets using the response matrix,” J. Comput. Phys. 199(1), 317–338 (2004). [CrossRef]

,61

F. K. Gruber and E. Marengo, “Reinterpretation and enhancement of signal-subspace-based imaging methods for extended scatterers,” SIAM J. Imaging Sci. 3(3), 434–461 (2010). [CrossRef]

]. These interesting problems for further study are currently being pursued.

Acknowledgements

The research is supported in part by US Army Medical Research and Materiel Command (USAMRMC) under Contract Number W81XWH-07-1-0454. M.X. acknowledges support by USAMRMC (W81XWH-10-1-0526) and NIH (1R15EB009224).

References and links

1.	A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]
2.	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(5), 426–428 (1995). [CrossRef] [PubMed]
3.	S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42(5), 841–853 (1997). [CrossRef] [PubMed]
4.	S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]
5.	W. Cai, S. K. Gayen, M. Xu, M. Zevallos, M. Alrubaiee, M. Lax, and R. R. Alfano, “Optical tomographic image reconstruction from ultrafast time-sliced transmission measurements,” Appl. Opt. 38(19), 4237–4246 (1999). [CrossRef] [PubMed]
6.	B. A. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron. 9(2), 199–209 (2003). [CrossRef]
7.	Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R. L. Barbour, “Frequency-domain optical imaging of absorption and scattering distributions by a Born iterative method,” J. Opt. Soc. Am. A 14(1), 325–342 (1997). [CrossRef]
8.	S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993). [CrossRef] [PubMed]
9.	N. Kroman, J. Wohlfahrt, H. T. Mouridsen, and M. Melbye, “Influence of tumor location on breast cancer prognosis,” Int. J. Cancer 105(4), 542–545 (2003). [CrossRef] [PubMed]
10.	D. S. Kepshire, S. C. Davis, H. Dehghani, K. D. Paulsen, and B. W. Pogue, “Subsurface diffuse optical tomography can localize absorber and fluorescent objects but recovered image sensitivity is nonlinear with depth,” Appl. Opt. 46(10), 1669–1678 (2007). [CrossRef] [PubMed]
11.	P. Mohajerani, A. A. Eftekhar, and A. Adibi, “Object localization in the presence of a strong heterogeneous background in fluorescent tomography,” J. Opt. Soc. Am. A 25(6), 1467–1479 (2008). [CrossRef] [PubMed]
12.	A. Godavarty, A. B. Thompson, R. Roy, M. Gurfinkel, M. J. Eppstein, C. Zhang, and E. M. Sevick-Muraca, “Diagnostic imaging of breast cancer using fluorescence-enhanced optical tomography: phantom studies,” J. Biomed. Opt. 9(3), 488–496 (2004). [CrossRef] [PubMed]
13.	Q. Zhao, L. Ji, and T. Jiang, “Improving depth resolution of diffuse optical tomography with a layer-based sigmoid adjustment method,” Opt. Express 15(7), 4018–4029 (2007). [CrossRef] [PubMed]
14.	M. Alrubaiee, M. Xu, S. K. Gayen, M. Brito, and R. R. Alfano, “Three-dimensional optical tomographic imaging of scattering objects in tissue-simulating turbid medium using independent component analysis,” Appl. Phys. Lett. 87(19), 191112 (2005). [CrossRef]
15.	M. Alrubaiee, M. Xu, S. K. Gayen, and R. R. Alfano, “Localization and cross section reconstruction of fluorescent targets in ex vivo breast tissue using independent component analysis,” Appl. Phys. Lett. 89(13), 133902 (2006). [CrossRef]
16.	M. Xu, M. Alrubaiee, S. K. Gayen, and R. R. Alfano, “Three-dimensional localization and optical imaging of objects in turbid media with independent component analysis,” Appl. Opt. 44(10), 1889–1897 (2005). [CrossRef] [PubMed]
17.	M. Xu, M. Alrubaiee, S. K. Gayen, and R. R. Alfano, “Optical diffuse imaging of an ex vivo model cancerous human breast using independent component analysis,” IEEE J. Sel. Top. Quantum Electron. 14(1), 43–49 (2008). [CrossRef]
18.	A. Poellinger, J. C. Martin, S. L. Ponder, T. Freund, B. Hamm, U. Bick, and F. Diekmann, “Near-infrared laser computed tomography of the breast first clinical experience,” Acad. Radiol. 15(12), 1545–1553 (2008). [CrossRef] [PubMed]
19.	V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A. 97(6), 2767–2772 (2000). [CrossRef] [PubMed]
20.	Q. Zhu, M. Huang, N. G. Chen, K. Zarfos, B. Jagjivan, M. Kane, P. Hedge, and S. H. Kurtzman, “Ultrasound-guided optical tomographic imaging of malignant and benign breast lesions: initial clinical results of 19 cases,” Neoplasia 5(5), 379–388 (2003). [PubMed]
21.	A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42(25), 5181–5190 (2003). [CrossRef] [PubMed]
22.	W. Cai, M. Alrubaiee, S. K. Gayen, M. Xu, and R. R. Alfano, “Three-dimensional optical tomography of objects in turbid media using the round-trip matrix,” Proc. SPIE 5693, 4–9 (2005). [CrossRef]
23.	B. Wu, M. Alrubaiee, W. Cai, M. Xu, and S. K. Gayen, “Optical imaging of objects in turbid media using principal component analysis and time reversal matrix methods,” in Computational Optical Sensing and Imaging, OSA Technical Digest (CD) (Optical Society of America, 2009), paper JTuC10. http://www.opticsinfobase.org/abstract.cfm?uri=COSI-2009-JTuC10
24.	B. Wu, W. Cai, M. Alrubaiee, M. Xu, and S. K. Gayen, “Three dimensional time reversal optical tomography,” Proc. SPIE 7892, 78920G (2011). [CrossRef]
25.	C. Prada, F. Wu, and M. Fink, “The iterative time reversal mirror: a solution to self-focusing in the pulse echo mode,” J. Acoust. Soc. Am. 90(2), 1119–1129 (1991). [CrossRef]
26.	M. Fink, C. Prada, F. Wu, and D. Cassereau, “Self-focusing in inhomogeneous media with time-reversal acoustic mirrors,” in IEEE Ultrasonics Symposium Proceedings (Montreal, Que., Canada, 1989), vol. 2, pp. 681–686.
27.	M. Fink, “Time reversal of ultrasonic fields. I. Basic principles,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 39(5), 555–566 (1992). [CrossRef] [PubMed]
28.	M. Fink, “Time reversal mirrors,” J. Phys. D Appl. Phys. 26(9), 1333–1350 (1993). [CrossRef]
29.	C. Prada, L. Thomas, and M. Fink, “The iterative time reversal process: analysis of the convergence,” J. Acoust. Soc. Am. 97(1), 62–71 (1995). [CrossRef]
30.	A. J. Devaney, E. A. Marengo, and F. K. Gruber, “Time-reversal-based imaging and inverse scattering of multiply scattering point targets,” J. Acoust. Soc. Am. 118(5), 3129–3138 (2005). [CrossRef]
31.	M. Fink, D. Cassereau, A. Derode, C. Prada, P. Roux, M. Tanter, J. L. Thomas, and F. Wu, “Time-reversed acoustics,” Rep. Prog. Phys. 63(12), 1933–1995 (2000). [CrossRef]
32.	W. A. Kuperman, W. S. Hodgkiss, H. C. Song, T. Akal, C. Ferla, and D. R. Jackson, “Phase conjugation in the ocean: experimental demonstration of an acoustic time reversal mirror,” J. Acoust. Soc. Am. 103(1), 25–40 (1998). [CrossRef]
33.	G. Lerosey, J. de Rosny, A. Tourin, and M. Fink, “Focusing beyond the diffraction limit with far-field time reversal,” Science 315(5815), 1120–1122 (2007). [CrossRef] [PubMed]
34.	A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm.
35.	H. Lev-Ari and A. J. Devaney, “The time reversal techniques re-interpreted: subspace-based signal processing for multi-static target location,” in Proceedings of the 1st IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM '00), (Cambridge, MA, USA, 2000) pp. 509–513.
36.	S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am. 113(5), 2742–2753 (2003). [CrossRef] [PubMed]
37.	F. K. Gruber, E. A. Marengo, and A. J. Devaney, “Time-reversal imaging with multiple signal classification considering multiple scattering between the targets,” J. Acoust. Soc. Am. 115(6), 3042–3047 (2004). [CrossRef]
38.	C. Prada and J. L. Thomas, “Experimental subwavelength localization of scatterers by decomposition of the time reversal operator interpreted as a covariance matrix,” J. Acoust. Soc. Am. 114(1), 235–243 (2003). [CrossRef] [PubMed]
39.	P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. 34(4), 561–580 (1992). [CrossRef]
40.	E. A. Marengo, F. K. Gruber, and F. Simonetti, “Time-reversal MUSIC imaging of extended targets,” IEEE Trans. Image Process. 16(8), 1967–1984 (2007). [CrossRef] [PubMed]
41.	A. Ishimaru, “Diffusion of a pulse in densely distributed scatterers,” J. Opt. Soc. Am. 68(8), 1045–1050 (1978). [CrossRef]
42.	K. Furutsu, “Diffusion equation derived from the space-time transport equation,” J. Opt. Soc. Am. 70(4), 360–366 (1980). [CrossRef]
43.	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(12), 2331–2336 (1989). [CrossRef] [PubMed]
44.	S. Chandrasekhar, Radiative Transfer (Clarendon Press, Oxford, 1950).
45.	A. Ishimaru, Wave Propagation and Scattering in Random Media, Volume 1: Single Scattering and Transport Theory (Academic, New York, 1978).
46.	S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25(12), 123010 (2009). [CrossRef]
47.	S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef] [PubMed]
48.	R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11(10), 2727–2741 (1994). [CrossRef] [PubMed]
49.	M. Lax, V. Nayaramamurti, and R. C. Fulton, “Classical diffusion photon transport in a slab,” in Laser Optics of Condensed Matter, J. L. Birman, H. Z. Cummins, and A. A. Kaplyanskii, eds. (Plenum, New York, 1987), pp. 229–237.
50.	C. Prada, S. Manneville, D. Spoliansky, and M. Fink, “Decomposition of the time reversal operator: detection and selective focusing on two scatterers,” J. Acoust. Soc. Am. 99(4), 2067–2076 (1996). [CrossRef]
51.	A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]
52.	N. Bourbaki, Topological Vector Spaces (Springer, 1987).
53.	H. J. van Staveren, C. J. M. Moes, J. van Marie, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. 30(31), 4507–4514 (1991). [CrossRef] [PubMed]
54.	C. Bordier, C. Andraud, E. Charron, J. Lafait, M. Anastasiadou, and A. Martino, “Illustration of a bimodal system in Intralipid 20% by polarized light scattering: experiments and modelling,” Appl. Phys., A Mater. Sci. Process. 94(2), 347–355 (2009). [CrossRef]
55.	S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16(7), 5048–5060 (2008). [CrossRef] [PubMed]
56.	M. Xu, Y. Pu, and W. Wang, “Clean image synthesis and target numerical marching for optical imaging with backscattering light,” Biomed. Opt. Express 2(4), 850–857 (2011). [CrossRef] [PubMed]
57.	T. Nielsen, B. Brendel, R. Ziegler, M. van Beek, F. Uhlemann, C. Bontus, and T. Koehler, “Linear image reconstruction for a diffuse optical mammography system in a noncompressed geometry using scattering fluid,” Appl. Opt. 48(10), D1–D13 (2009). [CrossRef] [PubMed]
58.	Y. Ardeshirpour, N. Biswal, A. Aguirre, and Q. Zhu, “Artifact reduction method in ultrasound-guided diffuse optical tomography using exogenous contrast agents,” J. Biomed. Opt. 16(4), 046015 (2011). [CrossRef] [PubMed]
59.	B. W. Pogue, M. S. Patterson, H. Jiang, and K. D. Paulsen, “Initial assessment of a simple system for frequency domain diffuse optical tomography,” Phys. Med. Biol. 40(10), 1709–1729 (1995). [CrossRef] [PubMed]
60.	S. Hou, K. Solna, and H. Zhao, “Imaging of location and geometry for extended targets using the response matrix,” J. Comput. Phys. 199(1), 317–338 (2004). [CrossRef]
61.	F. K. Gruber and E. Marengo, “Reinterpretation and enhancement of signal-subspace-based imaging methods for extended scatterers,” SIAM J. Imaging Sci. 3(3), 434–461 (2010). [CrossRef]

OCIS Codes
(110.3080) Imaging systems : Infrared imaging
(170.0170) Medical optics and biotechnology : Medical optics and biotechnology
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.3830) Medical optics and biotechnology : Mammography
(170.3880) Medical optics and biotechnology : Medical and biological imaging
(170.6510) Medical optics and biotechnology : Spectroscopy, tissue diagnostics
(110.0113) Imaging systems : Imaging through turbid media
(110.6955) Imaging systems : Tomographic imaging

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: August 31, 2011
Revised Manuscript: September 29, 2011
Manuscript Accepted: October 1, 2011
Published: October 21, 2011

Virtual Issues
Vol. 6, Iss. 11 Virtual Journal for Biomedical Optics

Citation
Binlin Wu, W. Cai, M. Alrubaiee, M. Xu, and S. K. Gayen, "Time reversal optical tomography: locating targets in a highly scattering turbid medium," Opt. Express 19, 21956-21976 (2011)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-19-22-21956

Sort: Author | Year | Journal | Reset

References

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol.50(4), R1–R43 (2005). [CrossRef] [PubMed]
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(5), 426–428 (1995). [CrossRef] [PubMed]
S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol.42(5), 841–853 (1997). [CrossRef] [PubMed]
S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl.15(2), R41–R93 (1999). [CrossRef]
W. Cai, S. K. Gayen, M. Xu, M. Zevallos, M. Alrubaiee, M. Lax, and R. R. Alfano, “Optical tomographic image reconstruction from ultrafast time-sliced transmission measurements,” Appl. Opt.38(19), 4237–4246 (1999). [CrossRef] [PubMed]
B. A. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron.9(2), 199–209 (2003). [CrossRef]
Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R. L. Barbour, “Frequency-domain optical imaging of absorption and scattering distributions by a Born iterative method,” J. Opt. Soc. Am. A14(1), 325–342 (1997). [CrossRef]
S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys.20(2), 299–309 (1993). [CrossRef] [PubMed]
N. Kroman, J. Wohlfahrt, H. T. Mouridsen, and M. Melbye, “Influence of tumor location on breast cancer prognosis,” Int. J. Cancer105(4), 542–545 (2003). [CrossRef] [PubMed]
D. S. Kepshire, S. C. Davis, H. Dehghani, K. D. Paulsen, and B. W. Pogue, “Subsurface diffuse optical tomography can localize absorber and fluorescent objects but recovered image sensitivity is nonlinear with depth,” Appl. Opt.46(10), 1669–1678 (2007). [CrossRef] [PubMed]
P. Mohajerani, A. A. Eftekhar, and A. Adibi, “Object localization in the presence of a strong heterogeneous background in fluorescent tomography,” J. Opt. Soc. Am. A25(6), 1467–1479 (2008). [CrossRef] [PubMed]
A. Godavarty, A. B. Thompson, R. Roy, M. Gurfinkel, M. J. Eppstein, C. Zhang, and E. M. Sevick-Muraca, “Diagnostic imaging of breast cancer using fluorescence-enhanced optical tomography: phantom studies,” J. Biomed. Opt.9(3), 488–496 (2004). [CrossRef] [PubMed]
Q. Zhao, L. Ji, and T. Jiang, “Improving depth resolution of diffuse optical tomography with a layer-based sigmoid adjustment method,” Opt. Express15(7), 4018–4029 (2007). [CrossRef] [PubMed]
M. Alrubaiee, M. Xu, S. K. Gayen, M. Brito, and R. R. Alfano, “Three-dimensional optical tomographic imaging of scattering objects in tissue-simulating turbid medium using independent component analysis,” Appl. Phys. Lett.87(19), 191112 (2005). [CrossRef]
M. Alrubaiee, M. Xu, S. K. Gayen, and R. R. Alfano, “Localization and cross section reconstruction of fluorescent targets in ex vivo breast tissue using independent component analysis,” Appl. Phys. Lett.89(13), 133902 (2006). [CrossRef]
M. Xu, M. Alrubaiee, S. K. Gayen, and R. R. Alfano, “Three-dimensional localization and optical imaging of objects in turbid media with independent component analysis,” Appl. Opt.44(10), 1889–1897 (2005). [CrossRef] [PubMed]
M. Xu, M. Alrubaiee, S. K. Gayen, and R. R. Alfano, “Optical diffuse imaging of an ex vivo model cancerous human breast using independent component analysis,” IEEE J. Sel. Top. Quantum Electron.14(1), 43–49 (2008). [CrossRef]
A. Poellinger, J. C. Martin, S. L. Ponder, T. Freund, B. Hamm, U. Bick, and F. Diekmann, “Near-infrared laser computed tomography of the breast first clinical experience,” Acad. Radiol.15(12), 1545–1553 (2008). [CrossRef] [PubMed]
V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A.97(6), 2767–2772 (2000). [CrossRef] [PubMed]
Q. Zhu, M. Huang, N. G. Chen, K. Zarfos, B. Jagjivan, M. Kane, P. Hedge, and S. H. Kurtzman, “Ultrasound-guided optical tomographic imaging of malignant and benign breast lesions: initial clinical results of 19 cases,” Neoplasia5(5), 379–388 (2003). [PubMed]
A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt.42(25), 5181–5190 (2003). [CrossRef] [PubMed]
W. Cai, M. Alrubaiee, S. K. Gayen, M. Xu, and R. R. Alfano, “Three-dimensional optical tomography of objects in turbid media using the round-trip matrix,” Proc. SPIE5693, 4–9 (2005). [CrossRef]
B. Wu, M. Alrubaiee, W. Cai, M. Xu, and S. K. Gayen, “Optical imaging of objects in turbid media using principal component analysis and time reversal matrix methods,” in Computational Optical Sensing and Imaging, OSA Technical Digest (CD) (Optical Society of America, 2009), paper JTuC10. http://www.opticsinfobase.org/abstract.cfm?uri=COSI-2009-JTuC10
B. Wu, W. Cai, M. Alrubaiee, M. Xu, and S. K. Gayen, “Three dimensional time reversal optical tomography,” Proc. SPIE7892, 78920G (2011). [CrossRef]
C. Prada, F. Wu, and M. Fink, “The iterative time reversal mirror: a solution to self-focusing in the pulse echo mode,” J. Acoust. Soc. Am.90(2), 1119–1129 (1991). [CrossRef]
M. Fink, C. Prada, F. Wu, and D. Cassereau, “Self-focusing in inhomogeneous media with time-reversal acoustic mirrors,” in IEEE Ultrasonics Symposium Proceedings (Montreal, Que., Canada, 1989), vol. 2, pp. 681–686.
M. Fink, “Time reversal of ultrasonic fields. I. Basic principles,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control39(5), 555–566 (1992). [CrossRef] [PubMed]
M. Fink, “Time reversal mirrors,” J. Phys. D Appl. Phys.26(9), 1333–1350 (1993). [CrossRef]
C. Prada, L. Thomas, and M. Fink, “The iterative time reversal process: analysis of the convergence,” J. Acoust. Soc. Am.97(1), 62–71 (1995). [CrossRef]
A. J. Devaney, E. A. Marengo, and F. K. Gruber, “Time-reversal-based imaging and inverse scattering of multiply scattering point targets,” J. Acoust. Soc. Am.118(5), 3129–3138 (2005). [CrossRef]
M. Fink, D. Cassereau, A. Derode, C. Prada, P. Roux, M. Tanter, J. L. Thomas, and F. Wu, “Time-reversed acoustics,” Rep. Prog. Phys.63(12), 1933–1995 (2000). [CrossRef]
W. A. Kuperman, W. S. Hodgkiss, H. C. Song, T. Akal, C. Ferla, and D. R. Jackson, “Phase conjugation in the ocean: experimental demonstration of an acoustic time reversal mirror,” J. Acoust. Soc. Am.103(1), 25–40 (1998). [CrossRef]
G. Lerosey, J. de Rosny, A. Tourin, and M. Fink, “Focusing beyond the diffraction limit with far-field time reversal,” Science315(5815), 1120–1122 (2007). [CrossRef] [PubMed]
A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and MUSIC” (2000). http://www.ece.neu.edu/faculty/devaney/ajd/preprints.htm .
H. Lev-Ari and A. J. Devaney, “The time reversal techniques re-interpreted: subspace-based signal processing for multi-static target location,” in Proceedings of the 1st IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM '00), (Cambridge, MA, USA, 2000) pp. 509–513.
S. K. Lehman and A. J. Devaney, “Transmission mode time-reversal super-resolution imaging,” J. Acoust. Soc. Am.113(5), 2742–2753 (2003). [CrossRef] [PubMed]
F. K. Gruber, E. A. Marengo, and A. J. Devaney, “Time-reversal imaging with multiple signal classification considering multiple scattering between the targets,” J. Acoust. Soc. Am.115(6), 3042–3047 (2004). [CrossRef]
C. Prada and J. L. Thomas, “Experimental subwavelength localization of scatterers by decomposition of the time reversal operator interpreted as a covariance matrix,” J. Acoust. Soc. Am.114(1), 235–243 (2003). [CrossRef] [PubMed]
P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev.34(4), 561–580 (1992). [CrossRef]
E. A. Marengo, F. K. Gruber, and F. Simonetti, “Time-reversal MUSIC imaging of extended targets,” IEEE Trans. Image Process.16(8), 1967–1984 (2007). [CrossRef] [PubMed]
A. Ishimaru, “Diffusion of a pulse in densely distributed scatterers,” J. Opt. Soc. Am.68(8), 1045–1050 (1978). [CrossRef]
K. Furutsu, “Diffusion equation derived from the space-time transport equation,” J. Opt. Soc. Am.70(4), 360–366 (1980). [CrossRef]
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(12), 2331–2336 (1989). [CrossRef] [PubMed]
S. Chandrasekhar, Radiative Transfer (Clarendon Press, Oxford, 1950).
A. Ishimaru, Wave Propagation and Scattering in Random Media, Volume 1: Single Scattering and Transport Theory (Academic, New York, 1978).
S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl.25(12), 123010 (2009). [CrossRef]
S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt.34(31), 7395–7409 (1995). [CrossRef] [PubMed]
R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A11(10), 2727–2741 (1994). [CrossRef] [PubMed]
M. Lax, V. Nayaramamurti, and R. C. Fulton, “Classical diffusion photon transport in a slab,” in Laser Optics of Condensed Matter, J. L. Birman, H. Z. Cummins, and A. A. Kaplyanskii, eds. (Plenum, New York, 1987), pp. 229–237.
C. Prada, S. Manneville, D. Spoliansky, and M. Fink, “Decomposition of the time reversal operator: detection and selective focusing on two scatterers,” J. Acoust. Soc. Am.99(4), 2067–2076 (1996). [CrossRef]
A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag.53(5), 1600–1610 (2005). [CrossRef]
N. Bourbaki, Topological Vector Spaces (Springer, 1987).
H. J. van Staveren, C. J. M. Moes, J. van Marie, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt.30(31), 4507–4514 (1991). [CrossRef] [PubMed]
C. Bordier, C. Andraud, E. Charron, J. Lafait, M. Anastasiadou, and A. Martino, “Illustration of a bimodal system in Intralipid 20% by polarized light scattering: experiments and modelling,” Appl. Phys., A Mater. Sci. Process.94(2), 347–355 (2009). [CrossRef]
S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express16(7), 5048–5060 (2008). [CrossRef] [PubMed]
M. Xu, Y. Pu, and W. Wang, “Clean image synthesis and target numerical marching for optical imaging with backscattering light,” Biomed. Opt. Express2(4), 850–857 (2011). [CrossRef] [PubMed]
T. Nielsen, B. Brendel, R. Ziegler, M. van Beek, F. Uhlemann, C. Bontus, and T. Koehler, “Linear image reconstruction for a diffuse optical mammography system in a noncompressed geometry using scattering fluid,” Appl. Opt.48(10), D1–D13 (2009). [CrossRef] [PubMed]
Y. Ardeshirpour, N. Biswal, A. Aguirre, and Q. Zhu, “Artifact reduction method in ultrasound-guided diffuse optical tomography using exogenous contrast agents,” J. Biomed. Opt.16(4), 046015 (2011). [CrossRef] [PubMed]
B. W. Pogue, M. S. Patterson, H. Jiang, and K. D. Paulsen, “Initial assessment of a simple system for frequency domain diffuse optical tomography,” Phys. Med. Biol.40(10), 1709–1729 (1995). [CrossRef] [PubMed]
S. Hou, K. Solna, and H. Zhao, “Imaging of location and geometry for extended targets using the response matrix,” J. Comput. Phys.199(1), 317–338 (2004). [CrossRef]
F. K. Gruber and E. Marengo, “Reinterpretation and enhancement of signal-subspace-based imaging methods for extended scatterers,” SIAM J. Imaging Sci.3(3), 434–461 (2010). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper

Figures


Fig. 1	Fig. 2	Fig. 3


Fig. 4	Fig. 5	Fig. 6

« Previous Article | Next Article »

OSA is a member of CrossRef.

Select a Conference
Session or Paper #
Year

OpticsInfoBase

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

Time reversal optical tomography: locating targets in a highly scattering turbid medium

Article Tools

Article Navigation

Abstract

1. Introduction

2. Formalism

2.1 Diffusion approximation, perturbation method and response matrix

2.2 Point Spread Functions

2.3 Time reversal and MUSIC

2.3.1 Well-resolved targets

2.3.2 Poorly-resolved Targets and MUSIC

3. Algorithm

4. Testing TROT with simulated data

5. Testing TROT using Experimental Data

5.1 Experimental materials and methods

5.2 Data Processing and Analysis

5.3 Results

5.3.1 Single absorptive target at different depths

5.3.2 Resolving two absorptive targets

5.3.3 Single scattering target at different depths

6. Discussion

Acknowledgements

References and links

References

Acad. Radiol.

Appl. Opt.

Appl. Phys. Lett.

Appl. Phys., A Mater. Sci. Process.

Biomed. Opt. Express

IEEE J. Sel. Top. Quantum Electron.

IEEE Trans. Antenn. Propag.

IEEE Trans. Image Process.

IEEE Trans. Ultrason. Ferroelectr. Freq. Control

Int. J. Cancer

Inverse Probl.

J. Acoust. Soc. Am.

J. Biomed. Opt.

J. Comput. Phys.

J. Opt. Soc. Am.

J. Opt. Soc. Am. A

J. Phys. D Appl. Phys.

Med. Phys.

Neoplasia

Opt. Express

Opt. Lett.

Phys. Med. Biol.

Proc. Natl. Acad. Sci. U.S.A.

Proc. SPIE

Rep. Prog. Phys.

Science

SIAM J. Imaging Sci.

SIAM Rev.

Other

Cited By

Figures