1.

Introduction

Optical imaging based on diffusive near-IR (NIR) light has shown promise in diagnostic and functional studies of tissue.¹ For example, it appears to be possible to optically differentiate breast cancer through determination of tissue parameters such as the oxygenation and the concentrations of hemoglobin, lipids, and water in breast tissue.² Another area of growing importance is optical imaging of the brain, where functional mapping has been able to identify regions of brain activation.³ However, NIR diffusive light imaging suffers from relatively low spatial resolution and poor depth discrimination. Previously proposed methods to enhance spatial resolution include time gating⁴ in the time domain, using high frequencies of intensity modulation⁵ or two-element phased-arrays⁶ in the frequency domain, or identifying optimal wavelengths in continuous-wave approaches.⁷ The issue of depth discrimination has been addressed at different levels by introducing off-axis detection,⁸ by applying two-layer or multilayer models,^{9, 10} or by full-fledged solutions of the inverse imaging problem.^{11, 12} Here, in an effort to develop a robust approach to diffuse optical imaging, we propose a multielement phased-array method that does not rely on any assumptions concerning boundary conditions, spatial features of tissue inhomogeneities, or any other kind of a priori information. We previously reported on the enhancement of spatial resolution afforded by this phased-array method,¹³ and here we present a depth discrimination study to demonstrate the potential of this method in tomographic sectioning of turbid media.

2.

Methods

2.1.

Phased-Array Intensity Associated with 2-D Arrays of Source/Detector Elements

The proposed phased-array approach consists of introducing arbitrary amplitude and phase factors by postprocessing the continuous-wave intensities associated with all individual source-detector pairs. Because of the asymmetrical source-detector configuration associated with a source array and a single detector, which is ultimately the reason for the depth discrimination capabilities, the idea is to use both a source array and a detector array. If one considers a single source (detector) element, the intensities associated with this specific source (detector) and with the elements of the array of detectors (sources) can be grouped according to the directions along which subsets of the detectors (sources) are aligned. If we label these directions as $d_j$ , the phased-array intensity associated with this specific source (detector) element is defined as

Eq. 1

$$I_{\mathrm{PA}}=\underset{d_j}{\max }\left\{\sum _i{A}_i^{\left(d_j\right)}\frac{I_i^{\left(d_j\right)}}{I_{0i}^{\left(d_j\right)}}\cos \right[{\alpha }_i^{\left(d_j\right)}\left]\right\},$$

where the intensity $I_i^{\left(d_j\right)}$ is normalized by the background intensity $I_{0i}^{\left(d_j\right)}$ , and $A_i^{\left(d_j\right)}$ and ${\alpha }_i^{\left(d_j\right)}$ are the amplitude and phase factors, respectively.¹³ The maximum over the various directions is taken to enhance the sensitivity to directional structures, which are most strongly detected by source/detector arrays that are arranged along a direction that is perpendicular to them.

For a proof of principle demonstration, and to better illustrate the features of the proposed phased-array approach, we present here a specific case of a 2-D array of seven continuous-wave light sources and a single optical detector. The seven light sources are electronically multiplexed to time share the optical detector. The bottom view of the planar array of illumination points is shown in Fig. 1 . This arrangement consists of the combination of three three-element source arrays along three different directions $(d_1,d_2,d_3)$ , where the central source element is in common for the three directions. For each direction $d_j$ , the phased-array intensity $\left(I_{\mathrm{PA}}^{\left(d_j\right)}\right)$ is defined by the expression¹³

Eq. 2

$$I_{\mathrm{PA}}^{\left(d_j\right)}=\frac{I_1^{\left(d_j\right)}}{I_{01}^{\left(d_j\right)}}+\frac{I_3^{\left(d_j\right)}}{I_{03}^{\left(d_j\right)}}-2\frac{I_2^{\left(d_j\right)}}{I_{02}^{\left(d_j\right)}}(j=1,2,3),$$

where the amplitude and phase factors of Eq. 1 along each direction $d_j$ are set to be $A_1=A_3=1$ , $A_2=2$ , ${\alpha }_1={\alpha }_2=0$ , and ${\alpha }_2=\pi$ . As shown by Eq. 1, the phased-array intensity associated with this seven-element source array and one detector is obtained by taking the maximum of the three-element phased-array intensities along the three directions:

Eq. 3

$$I_{\mathrm{PA}}=\max [I_{\mathrm{PA}}^{\left(d_1\right)},I_{\mathrm{PA}}^{\left(d_2\right)},I_{\mathrm{PA}}^{\left(d_3\right)}).$$

In this source-array single-detector configuration, the detector is placed across the center source at a distance of $6\mathrm{cm}$ , which is representative of the source-detector separations used in typical noninvasive, diffuse optical studies tissues.

Fig. 1

Bottom view of the planar array of illumination fibers, where $d_1$ , $d_2$ , and $d_3$ indicate three directions, and the numbers indicate different fibers in each direction. The distance between two adjacent fibers along any of the three directions is $1\mathrm{cm}$ .

The spatial distribution of the sensitivity function can be considered in the cases of a single source-detector pair and a three-element, linear phased-array, $S$ and $S_{\mathrm{PA}}$ , respectively. In a highly scattering medium the sensitivity functions $S$ and $S_{\mathrm{PA}}$ can be defined as

Eq. 4

$$S=\mid \mathrm{\Delta }I(y,z)∕\mathrm{\Delta }I(0,0)\mid ,$$

Eq. 5

$$S_{\mathrm{PA}}=\mid \mathrm{\Delta }{I}_{\mathrm{PA}}(y,z)∕\mathrm{\Delta }{I}_{\mathrm{PA}}(0,0)\mid ,$$

where $\mathrm{\Delta }I(y,z)$ and $\mathrm{\Delta }{I}_{\mathrm{PA}}(y,z)$ are the changes in single source-detector intensity and phased-array intensity, respectively, caused by a point-like absorbing object at position $(y,z)$ . The $y$ axis is oriented along the $d_1$ direction, as defined in Fig. 1. The spatial distribution of the sensitivity function in the plane defined by the detector and by the linear source-array (plane $x=0$ ) was obtained by using first-order perturbation theory of the diffusion equation in the infinite-medium geometry,¹³ and it is shown in Fig. 2 . The region of sensitivity of the single source-detector intensity is symmetrical with respect to the midplane $z=0$ (which is perpendicular to the plane of Fig. 2). By contrast, in the phased-array case, the sensitivity close to the source array is much higher than that in the region close to the detector. The reason why the contrast is higher close to the source array is that the intensities detected from the individual sources become increasingly similar for objects closer to the single detector, and they tend to cancel each other out in the phased-array intensity of Eq. 2. By using this asymmetric character of the phased-array method properly, it is possible to discriminate objects at different depths in the medium.

Fig. 2

Spatial distribution of the sensitivity function on the $y\text{-}z$ plane for (a) one source and one detector located at $(y_s,z_s)=(0,-3)$ and $(y_d,z_d)=(0,3)$ , respectively, and (b) a linear, three-source array consisting of sources located at $(y_1,z_1)=(-1,-3)$ , $(y_2,z_2)=(0,-3)$ , and $(y_3,z_3)=(1,-3)$ . In both panels, the sensitivity function is defined relative to the value at $(y,z)=(0,0)$ (i.e., the center point of the plane). The optical coefficients used in Eq. 3 to calculate these sensitivity functions using diffusion theory are ${\mu }_{a0}=0.03{\mathrm{cm}}^{-1}$ and ${\mu }_{s0}^{\prime }=14{\mathrm{cm}}^{-1}$ .

2.2.

Experimental Setup

For the optical measurements, we used seven multiplexed (i.e., turned on and off in sequence) laser diodes emitting at $690\mathrm{nm}$ from a NIR tissue spectrometer (OxiplexTS, ISS, Inc., Champaign, Illinois). These laser diodes were coupled to multimode $400-\mathrm{\mu }\mathrm{m}-\mathrm{diam}$ optical fibers, whose emitting ends were arranged in a 2-D source array according to the layout of Fig. 1. The collection optical fiber, a fiber bundle with an internal diameter of $3\mathrm{mm}$ , was facing the central source fiber at a distance of $6\mathrm{cm}$ . The illumination and collection optical fibers were immersed in a highly scattering liquid sample (2% milk diluted 50% by volume in water, with added black India ink) having an absorption coefficient of $0.030\pm 0.002{\mathrm{cm}}^{-1}$ and a reduced scattering coefficient of $14\pm 1{\mathrm{cm}}^{-1}$ . Three $0.63-\mathrm{cm}-\mathrm{diam}$ black cylinders were placed between the source array and detector with their axes parallel to the plane of the source array, and at different depths or distances from the array. The experimental setup and the cylinders are shown in Fig. 3 . Specifically, one cylinder (number 2 in Fig. 4 ) was oriented along the $y$ axis and placed at the midplane, i.e., equidistant $\left(3.0\mathrm{cm}\right)$ from the source array and the collection fiber. A second cylinder (number 1 in Fig. 4) was oriented along the $x$ axis and placed $1.0\mathrm{cm}$ above the midplane. The third cylinder (number 3 in Fig. 4) was oriented at $40\mathrm{deg}$ with respect to the $x$ axis and placed $1.0\mathrm{cm}$ below the midplane. We conceptually divide the medium between the source array and the detector into three layers and we refer to these layers (and the cylinders in them) as top, middle, and bottom, respectively. We performed a 2-D tandem scan of the source fiber array and the collection fiber in the $x\text{-}y$ plane, where the sources and the detector were scanned together as a rigid body. The scanning speed along the $x$ axis was $0.26\mathrm{cm}∕\mathrm{s}$ , while the optical sampling rate was $14.3\mathrm{Hz}$ , so that we acquired one optical data point (for each of the seven source elements) every $182\mathrm{\mu }\mathrm{m}$ along the $x$ axis. In more detail, the acquisition of intensity was performed by time multiplexing the seven sources, where each source was on for about $10\mathrm{ms}$ , which corresponded to $26\mathrm{\mu }\mathrm{m}$ along the scanning coordinate. By assigning the detector position to the $x$ -scanning coordinate, the offset introduced by the continuous scanning in the intensity measurements from different source-detector pairs was at most $182\mathrm{\mu }\mathrm{m}$ , which did not affect the content of our results. The spatial sampling density along the $y$ axis was $1{\mathrm{mm}}^{-1}$ , which means that we performed linear scans along the $x$ axis every $1\mathrm{mm}$ along the $y$ axis.

Fig. 3

Side view of the sample and experimental setup showing the seven-element, 2-D source array, the single detector, and the three cylinders embedded at different depths within the highly scattering medium (2% milk, water, and black India ink): PMT, the photomultiplier tube. The numbers 1, 2, and 3 label three cylinders at different depths along the $z$ axis. If we imagine dividing the medium between source and detector in three layers, cylinder 1 is in the top layer, cylinder 2 is in the middle layer, and cylinder 3 is in the bottom layer.

Fig. 4

Top view of the scanning area. The numbers 1, 2, and 3 label three cylinders at different depths along the $z$ axis. Cylinder 1 is in the top layer; cylinder 2 is in the middle layer, and cylinder 3 is in the bottom layer (Fig. 3). The dashed square represents the real scanned area. The scan began from the left lower corner.

3.

Results

Figure 5 shows optical images based on the normalized intensity [Fig. 5(a)] and its spatial second derivative [Fig. 5(b)] collected from the center source element. The image of the spatial second derivative is obtained by taking the maximum among the second derivatives in four directions ( $x$ , $y$ , and two other mutually perpendicular directions oriented at $45\mathrm{deg}$ with respect to $x$ and $y$ ). The step of derivation was chosen as $4\mathrm{mm}$ in all the directions. Figures 5(c) and 5(d) show the phased-array images with the phased array on the top [Fig. 5(c)] and on the bottom [Fig. 5(d)]. The spatial second derivative improves the spatial resolution in diffuse optical imaging,¹⁴ and the phased-array approach proposed here takes advantage of this feature because of the formal similarity between the phased array intensity of Eq. 2 and the discrete spatial second derivative.¹³ The second derivative enhances the spatial features in the image and exploits some key features of photon migration not only for edge detection but also for the discrimination of single versus multiple defects. However, without off-axis measurements, the images of Fig. 5 are 2-D projections and do not provide any depth information. Figures 5(c) and 5(d) demonstrate the different sensitivity of the phased-array intensity to structures at different depths within the medium. Such selective sensitivity is exploited by taking the difference of the phased-array intensities measured with the phased array on the top and on the bottom $(I_{\mathrm{PA}}^{\mathrm{top}}-I_{\mathrm{PA}}^{\mathrm{bottom}})$ , as shown in Fig. 6(a) . Positive values (white) in Fig. 6(a) correspond to top structures, whereas negative values (black) correspond to bottom structures. The sum of the phased-array intensities $(I_{\mathrm{PA}}^{\mathrm{top}}+I_{\mathrm{PA}}^{\mathrm{bottom}})$ shown in Fig. 6(b) identifies all detected structures, regardless of their depth.

Fig. 6

(a) Difference $(I_{\mathrm{PA}}^{\mathrm{top}}-I_{\mathrm{PA}}^{\mathrm{bottom}})$ and (b) sum $(I_{\mathrm{PA}}^{\mathrm{top}}+I_{\mathrm{PA}}^{\mathrm{bottom}})$ of the phased-array intensities measured when the phased array was on the top and on the bottom, respectively.

Fig. 5

Imaging based on (a) the intensity from the center-source (normalized to the peak value), (b) its spatial second derivative, (c) the phased-array intensity with the phased-array on top, and (d) the phased-array intensity with the phased array on the bottom. On images (c) and (d) we superimposed the segments of the scan defining the $c$ -axis of Fig. 7. The segment of the scan parallel to the $y$ axis was taken along the line $x=10\mathrm{mm}$ and $50<y<90\mathrm{mm}$ , while the segment of the scan parallel to the $x$ axis was taken along the line $y=50\mathrm{mm}$ and $20<x<100\mathrm{mm}$ .

Fig. 7

Principle diagram to separate the intensity of three layers: (a) phased-array intensity when the phased array is above the top layer; (b) phased-array intensity when the phased array is below the bottom layer; (c) difference of (a) and (b); (d) dot curve represents the sum of (a) and (b), and the solid curve is the middle layer [calculated using Eq. 8] with some side contributions for top and bottom cylinders.

The two-step depth discrimination imaging approach described by Eqs. 6 to 8 is illustrated in Figs. 7 and 8 . Figure 7 shows a line trace of the phased-array intensities with the phased array on top [Fig. 7(a)] and on the bottom [Fig. 7(b)], their difference [Fig. 7(c)], and their sum [Fig. 7(d)]. The line trace was obtained by joining two segments of the scans, one parallel to the $y$ axis and one parallel to the $x$ axis for the purpose of clarifying the principle of the algorithm. More precisely, the segment of the scan parallel to the $y$ axis was taken along the line $x=10\mathrm{mm}$ and $50<y<90\mathrm{mm}$ , while the segment of the scan parallel to the $x$ axis was taken along the line $y=50\mathrm{mm}$ and $20<x<100\mathrm{mm}$ [these segments are shown in Figs. 5(c) and 5(d)]. Therefore the $c$ axis in these figures corresponds to the coordinate along this line obtained by joining two perpendicular segments of the scan that crosses the three cylinders (top, middle, and bottom). The threshold $\sigma$ is chosen just above the noise level. In Fig. 7(c), any peaks above the positive threshold are assigned to the top layer, whereas peaks below the negative threshold are assigned to the bottom layer. Any peaks in Fig. 7(d) at locations that do not correspond to top or bottom layers are assigned to the middle layer. The images that result from this procedure are shown in Fig. 8 for the top [Fig. 8(a)], bottom [Fig. 8(b)], and middle [Fig. 8(c)] layers. The top and bottom cylinders are isolated quite clearly, while the middle cylinder is partially obscured at the sections where it overlaps with the other two cylinders. We also note that the area where the top and bottom cylinders cross is assigned to the middle layer because, by the way the algorithm is currently devised, anything that appears symmetrically on the top and bottom layers will cancel out in the difference $I_{\mathrm{PA}}^{\mathrm{top}}-I_{\mathrm{PA}}^{\mathrm{bottom}}$ and would be erroneously assigned to the middle layer.

Fig. 8

Depth-resolved images (a) of the top layer, (b) of the bottom layer, and (c) of the middle layer.

4.

Discussion

The top and bottom cylinders as imaged in Figs. 8(a) and 8(b) lack a section that corresponds to the area where the two cylinder projections overlap on the $x\text{-}y$ plane. This overlapping area is assigned to the middle layer image [Fig. 8(c)]. This is a limitation of the proposed algorithm: whenever we have symmetrical structures with the same optical properties, belonging to the top and bottom layer, they will be absent in the top and bottom layer and will be erroneously assigned to the middle layer. In this case (which might be rather unlikely in real tissues), one way to have an idea that we are in the presence of an artifact is to look at the raw phased-array images (that is, those before subtraction or addition of the intensities), where the overlapping region appears with high contrast in both the phased-array-on-top image and phased-array-on-bottom image. Even in this case, however, we cannot exclude the presence of a single defect placed in the midplane instead of two defects placed symmetrically in the top and bottom layers. However, if this would be the case, the contrast of the single defect should be much higher than the defects placed in the top and bottom planes. We are currently working on an improvement of the algorithm to achieve a better discrimination of overlapping structures [that is, having the same $(x,y)$ coordinate and located at different depths]. We notice that the current method should be able to coarsely discriminate those structures whenever they are not symmetrical with respect to the middle plane or have different contrast with respect to the background medium. One way to improve its performance and to better discriminate directional structure might be to use the information of the single phased-array intensities (i.e., along different directions) and compare this information with the maximum and minimum intensities. In this perspective, it might be useful to also choose a direction-dependent value of $\sigma$ . The other obvious way to improve the performance of the current algorithm is to take another projection perpendicular to the previous one (the $y\text{-}z$ plane of Fig. 4). Finally, we notice that the value of $\sigma$ does not have a one-to-one correspondence with the location of structures at different depths because the difference of phased-array intensities when the array is on top and on bottom depends both on the optical contrast between defects and background and on the absolute values of the background optical properties. Nevertheless, by varying the value of $\sigma$ we can inverse or decrease the thickness of the middle layer, thus providing a variable geometrical sectioning capability.

5.

Conclusion

We presented a multielement phased-array approach to diffuse optical imaging that is based on postprocessing cw data and that has the potential to discriminate the depth of multiple absorbing defects in highly scattering media. This multielement phased-array approach can enhance spatial resolution similarly to the way that the spatial second-derivative does. Furthermore, it can obtain depth discrimination because of a variable sensitivity as a function of depth. A simultaneous implementation of 2-D arrays of light sources and optical detectors would enable the most effective implementation of the proposed approach to depth discrimination. We also discussed possible ways to improve the current performance of depth discrimination for the case of overlapping structures.