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–5] 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] and direct linear inversion of 3-D matrices [7], are time consuming. The iterative methods [7,8] 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], 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–8]. 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]. 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]. 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]. 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]. 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–17]. 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–21].

In this article we report on the development of a time reversal optical tomography (TROT) [22–24] 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–28]. The theoretical and numerical techniques involved in time reversal have been used for applications involving both acoustic waves and electromagnetic waves (radar) [28–33].

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,35]. 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]. They also generalized the theory which was based on distorted wave Born approximation (DWBA) to account for multiple scattering between the targets [37]. 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,35]. TR-MUSIC provided higher spatial resolution than the conventional TR imaging, especially in the case where targets were not well resolved [34,35,38].

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,40]. 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–43] of the radiative transfer equation (RTE) [44,45]. 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,47]

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, $\Delta \varphi =\varphi -{\varphi }_0$, is found, where $\varphi$ is the light intensity distribution measured on the sample boundary with targets embedded in the scattering medium and ${\varphi }_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 $-\Delta \varphi$, to describe the transport of light from different sources through the embedded objects to the array of detectors [22,36].

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

For a homogeneous background medium, the rank R of matrix K, is equal to the dimension of the source array vector space ${\mathcal{G}}_s$ spanned by $g_s(r_m)$, and also equal to the dimension of the detector array vector space ${\mathcal{G}}_d$ spanned by $g_d(r_m)$, where ${\mathcal{G}}_s\subseteq C^{N_s}$ and ${\mathcal{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

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,50]. So the signal evaluated at r is [34]

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,34,38,50,51]. For frequency-domain data, $T_{SDDS}=K^{†}K$, and $T_{DSSD}={(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_{SDDS}=K^{T}K$, and $T_{DSSD}=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,\cdots ,N_s$) and u_i ($i=1,\cdots ,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,\cdots ,M$, and ${\lambda }_j\approx 0$, when $j=M+1,\cdots ,N_s$ for T_SDDS and $j=M+1,\cdots ,N_d$ for T_DSSD. The eigen system $\{v_j,u_j,{\lambda }_j>0\}$, $j=1,\cdots ,M$, is related to the targets. The TR matrix T_SDDS can be written as [22,34]

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^{th}$and $m{\text{'}}^{th}$ targets ($m\ne m\text{'}$) are well resolved, defined by the conditions: $H_s(X_m,X_{m\text{'}})\approx 0$ and $H_d(X_m,X_{m\text{'}})\approx 0$, i.e. the GFVs at $X_m$ and $X_{m\text{'}}$ are orthogonal, $<g_d(X_m),g_d(X_{m\text{'}})>=H_d(X_{m\text{'}},X_m)={\Vert g_d\left(X_m\right)\Vert }^2{\delta }_{mm\text{'}}$. So we have

The eigenvectors are

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

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,\cdots ,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,36,51]. 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 $\left\{v_j{}^{*}\right\}$ ($j=1,\cdots ,N_s$) and $\left\{u_j{}^{*}\right\}$ ($j=1,\cdots ,N_d$) span the spaces $C^{N_s}$ and $C^{N_d}$ associated with the source and detector arrays, respectively. While $\left\{v_j{}^{*}\right\}$ and $\left\{u_j{}^{*}\right\}$, with ${\lambda }_j>0$, form the signal subspaces on the source and detector arrays, ${\mathcal{S}}^s=\left\{v_j{}^{*}\right\}$ and ${\mathcal{S}}^d=\left\{u_j{}^{*}\right\}$ ($j=1,\cdots ,M$), respectively; $\left\{v_j{}^{*}\right\}$ and $\left\{u_j{}^{*}\right\}$, with ${\lambda }_j\approx 0$, form the noise subspaces, ${\mathcal{N}}^s=\left\{v_j{}^{*}\right\}$ ($j=M+1,\cdots ,N_s$) and ${\mathcal{N}}^d=\left\{u_j{}^{*}\right\}$ ($j=M+1,\cdots ,N_d$), respectively. Thus $C^{N_s}={\mathcal{S}}^s\oplus {\mathcal{N}}^s$ and $C^{N_d}={\mathcal{S}}^d\oplus {\mathcal{N}}^d$ [36,51]. Since the dimensions of the signal subspaces ${\mathcal{S}}^s$ and ${\mathcal{S}}^d$ and of the GFV spaces ${\mathcal{G}}_s$ and ${\mathcal{G}}_d$ are all equal to M, ${\mathcal{G}}_s\equiv {\mathcal{S}}^s$ and ${\mathcal{G}}_d\equiv {\mathcal{S}}^d$ [51]. The GFVs, $g_s(X_m)$ and $g_d(X_m)$, $m=1,\cdots ,M$, are linear combinations of v_j* and u_j*, $j=1,\cdots ,M$, respectively. Therefore, $g_s(X_m)\in {\mathcal{S}}^s$ and $g_d(X_m)\in {\mathcal{S}}^d$, $m=1,\cdots ,M$, associated with m^th target are orthogonal to $v_j{}^{*}\in {\mathcal{N}}^s$ ($j=M+1,\cdots ,N_s$) and $u_j{}^{*}\in {\mathcal{N}}^d$ ($j=M+1,\cdots ,N_d$), respectively:

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

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,36,51], $Q_s(X_p)$and $Q_d(X_p)$ can be calculated as

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,51] associated with the source array,

When scattering targets are concerned, the GFVs ${\partial }_{\alpha }g$ ($\alpha =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 }_{\alpha }g$, ($\alpha =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$.

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

View Table | View all tables in this article

 

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

Download Full Size | PDF

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,54] 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.

Download Full Size | PDF

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 $\Delta \varphi$ 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, $-\Delta \varphi$. 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]. 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]. 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]. 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]. 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]. 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–30,34,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.

Download Full Size | PDF

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

View Table | View all tables in this article

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

Download Full Size | PDF

Table 3. Positions of one target located at different depths

View Table | View all tables in this article

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

Download Full Size | PDF

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

View Table | View all tables in this article

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

Download Full Size | PDF

Table 5. Positions of one scattering target located at different depths

View Table | View all tables in this article

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,60,61]. These interesting problems for further study are currently being pursued.