In the last decade, the possibility to quantitatively reconstruct absorbing, scattering, and fluorescent inclusions within in vivo organisms attracted a great deal of interest for diagnostic purposes (e.g., tumor detection) [1], functional studies (e.g., brain oximetry) [2], and molecular imaging on small animals (e.g., pharmacological research) [3]. The general measurement scheme consists of illuminating a sample and detecting the diffused light exiting it. Then, by solving the inverse problem, based on a model of photon propagation through the biological tissue, the optical parameters in each point of the sample can be reconstructed quantitatively. These modalities are usually referred to as diffuse optical tomography (DOT) or fluorescence molecular tomography (FMT) when, respectively, the absorption/scattering or fluorescence properties are reconstructed. The performance of DOT/FMT is mainly characterized by its capability to resolve the position and shape of inhomogeneities inside the tissue, and, consequently, improve the quantification capability of their optical parameters. Previous studies have demonstrated the importance of a dense source/detector [4] and a multiple view measurement scheme [5,6] in order to increase the tomographic spatial resolution. Moreover, further data, such as spectral and temporal information, are crucial [7,8]. Temporal information provides three main advantages: (i) better disentanglement of absorption/scattering properties, (ii) temporal encoding of photon depth, and (iii) fluorescence lifetime quantification in the case of FMT. Spectral information (i.e., different excitation/detection wavelength) allows one to discriminate among tissue chromophores. Hence, DOT/FMT turns out to be a highly multidimensional problem with the drawback of a huge data set being generated. This represents a practical limitation of these techniques because of the extremely long acquisition and computational times, which are not typically compatible with clinical and preclinical needs. Hence, a reduction of the acquired data set by preserving the spatial resolution, or, more generally, the data set information content, is highly desirable.

Following this concept, different studies have recently exploited the fact that a highly scattering medium (such as biological tissue) behaves as a low-pass filter in the spatial domain. Hence, a few illumination patterns, instead of the more typical raster scanning approach, can be adopted without losing significant spatial information [8–10]. This in turn leads to a reduction of the data set dimension and, consequently, of the acquisition and computational times. Recent studies have exploited such an approach in both imaging and tomographic schemes, and detection is generally carried out by a parallel detector such as a CCD, CMOS, or gated cameras [11]. Moreover, the use of a wide-field approach (such as the case of structured illumination) allows one to illuminate the sample with high power without exceeding safety limits. This improves the signal-to-noise ratio.

Recently, a patterned detection [12], following the single-pixel camera (SPC) scheme [13], has been proposed for FMT applications, as well as for PhotoAcoustics [14]. Basically, the image of the diffused light exiting the sample is modulated spatially and subsequently focused on a single element detector. This operation is equivalent to projecting the image on an element of a base set, such as Fourier, Wavelets, or Hadamard patterns. By repeated acquisitions for different base elements, it is possible to recover the same image as would be measured in a conventional pixel basis. Due to the fact that a highly scattering medium acts as a low-pass filter in the spatial frequency domain, just a few frequencies are needed. This approach has the great advantage of exploiting the superior characteristics of a single detector (e.g., higher temporal resolution and larger spectral bandwidth) at a lower cost with respect to a parallel device. Moreover, compared with raster scanning, a further advantage is the acquisition speed given by a wide-field detection analogous to a structured illumination approach. Finally, it is worth mentioning that both structured illumination and detection open the possibility to acquire images and reconstructions with increasing spatial details by increasing the number of measurements.

Whereas the patterned detection approach has been successfully demonstrated for fluorescence optical tomography in a single view [12], a fully tomographic modality requires multiple views to reduce the ill-posedness in depth resolution, which leads to a more challenging experimental arrangement. In this work, we propose a multiple-view, time-domain compressed sensing DOT system exploiting Hadamard patterns in both the illumination and collection planes, and applicable to non-planar geometries. The system has been experimentally validated on tissue phantoms with absorbing inclusions, demonstrating both imaging and tomographic capabilities.

The experimental setup is sketched schematically in Fig. 1. The sample is illuminated by a pulsed structured light while detection is carried out by either a time-resolved (TR) SPC or a continuous-wave (CW) parallel detector. The sample is placed on a rotational stage to allow different view acquisitions. By means of an acousto-optic tunable filter, light pulses at 650 nm are selected spectrally from a picosecond pulsed supercontinuum (repetition rate of 80 MHz) laser source (SuperK Extreme, NKT). Structured illumination is carried out by a digital micromirror device (DMD Discovery kit 1100, Vialux), which spatially modulates the light, and an objective lens ($f=50\mathrm{cm}$) to create an image over an area of $3\mathrm{cm}\times 3\mathrm{cm}$ of the sample. The diffused light, exiting the sample over an area of about $2\mathrm{cm}\times 2\mathrm{cm}$, is imaged by a lens ($f=60\mathrm{cm}$) on a second digital micromirror device (DMD) (DMD Discovery 4100, Vialux). A flip mirror allows us to image the DMD plane on either a low-noise 16-bit cooled CCD camera (Versarray 512, Princeton Instruments) or a single element detector. The latter consists of a long-working-distance objective ($10\times /0.25$), which focuses the light reflected by the second DMD onto an optical fiber of 1 mm diameter. The light exiting the fiber is finally detected by a photomultiplier (PMT) (HPM-100-50, Becker and Hickl) connected to a time-correlated single-photon counting (TCSPC) board, which samples the temporal profile of the diffuse light. The system is fully computer-controlled by a homemade LabView software, enabling automated acquisition of the whole data set (illumination/detection patterns, sample rotation, and acquisition). The sample is a homogeneous cylindrical tissue mimicking a phantom ($Ø=20\mathrm{mm}$, height 50 mm) made of epoxy resin, ${\mathrm{TiO}}_2$ (as a scatterer), and a toner (as an absorber). By means of a time-resolved spectroscopy system [15], the optical parameters were measured: absorption coefficient (${\mu }_a$) of about $0.01{\mathrm{mm}}^{-1}$ and reduced scattering coefficient (${\mu }_s^{\prime }$) of about $1{\mathrm{mm}}^{-1}$. Two holes, drilled into the sample ($Ø=1.6\mathrm{mm}$), allowed us to insert either solid or liquid absorbing inclusions. In particular, to better simulate a realistic perturbation, three solutions of calibrated ink and Intralipid have been prepared [16], giving ${\mu }_s^{\prime }\sim 1{\mathrm{mm}}^{-1}$ (the same as that of the background) and ${\mu }_a$ of $0.05/0.05{\mathrm{mm}}^{-1}$ (Exp 1), $0.05/0.1{\mathrm{mm}}^{-1}$ (Exp 2), and $0.05/0.03{\mathrm{mm}}^{-1}$ (Exp 3) for inclusions A and B.

 

Fig. 1. Experimental setup.

Download Full Size | PDF

Initially, images have been acquired by means of the CCD camera on the detection side and a low-cost camera on the illumination side, to register the illumination/detection area over the sample. Then 360 shadows of the object (every 1°) have been acquired to create the mesh [6]. It is worth emphasizing that precise calibration is critical to achieving an accurate simulation of the forward problem, which in turn is a prerequisite to obtain a high-quality tomographic reconstruction.

Measurements have been performed on the phantom with and without the absorbing inclusions. The acquisition procedure is carried out by a complete 360° rotation of the sample with steps of 45° (eight views). For each view, $8\times 8$ ordered Walsh–Hadamard (WH) patterns, covering an area of $1.3\mathrm{cm}\times 1.3\mathrm{cm}$ on the sample, have been used for both illumination and detection. Each WH pattern consists of two states ($-1$ to $+1$). Hence, two positive patterns (ranging from 0 to $+1$), complementary to one another, have been acquired and properly subtracted to obtain the desired WH pattern. The acquisition time for each pattern is 1 s with 800 kHz as maximum count rate to fulfill the single-photon statistics. This last parameter is, indeed, the limiting factor on the overall acquisition time, which is about 25 min. A full-pixel image can be recovered by applying a fast WH inverse transform to the detected data [17].

For the reconstruction of the absorption map in the volume, the following objective function has been constructed:

First, measurements have been carried out using black solid rods as inclusions to demonstrate the imaging capability of the system and to estimate the number of patterns to be used in the tomographic reconstruction. In particular, time-resolved data acquired by the SPC have been integrated over time to obtain CW data and compared with the CCD images. An example of the images acquired by the CCD and the ones based on the SPC (by spatially modulating the detection) is shown in Fig. 2. We observe a good agreement between the two images, which is improved by increasing the number of adopted patterns as reported in Fig. 2, where the root-mean-square error (RMSE) is reported as a function of the WH pattern order. In particular, we do not observe a significant improvement for a WH pattern order higher than 8. Moreover, it is possible to observe that the RMSE plot for the inhomogeneous phantom presents a higher variability among the different views with respect to the homogeneous case. It is worth stressing that the number of required patterns strongly depends on the optical parameters/shape of the sample and the position/dimension of the inclusions.

 

Fig. 2. Comparison between CCD and SPC images outcoming from the homogeneous (first row) and heterogeneous (second row) phantoms. First column: image measured with the CCD. Second column: image measured by the SPC using $16\times 16$ WH patterns for a single view. Third column: RMSE as a function of WH order for eight views.

Download Full Size | PDF

In order to explore the imaging capability of the proposed method, the relative contrast, calculated as the difference between heterogeneous and homogeneous images divided by the homogeneous one, provided by one solid inclusion, in both the CCD and SPC images, is reported, for eight different views, in Fig. 3. In particular, three cases are reported: (i) the sample is illuminated with $16\times 16$ ordered WH patterns while the detection side has a uniform square pattern (Fig. 3, first row); (ii) the sample is illuminated with a uniform square pattern while the detection side is spatially modulated by $16\times 16$ ordered WH patterns (Fig. 3, second row); (iii) CCD images by using a uniform square illumination pattern are also reported (Fig. 3, third row). Cases (ii) and (iii) show good agreement; in particular, the presence of the absorbing inclusion can be clearly observed when it is located, during sample rotation, closer to the detector (see the vertical blue bar in the image at 0° or 45° for OUT). On the contrary, for the other views, the inclusion cannot be observed because of the scattering. In case (i), there is no correspondence between the images acquired by the CCD and SPC. In particular, we observe that, by modulating the illumination, we can better observe the inclusion for views where it is closer to the illumination source (see the vertical blue bar in the image at 180° or 225° for IN). In fact, the SPC approach measures the integral of the signal; then the imaging capability is not influenced by the scattering events followed by photons after impinging on the inclusion as occurs for the CCD [20]. These examples demonstrate the imaging capability of the proposed method and, in particular, the importance of the choice of illumination/detection patterns according to the view, the sample (shape and optical parameters) and inclusions, for both imaging and reconstruction.

 

Fig. 3. Relative contrast for different rotation views (columns) in the case of modulated illumination (first row), modulated detection (second row), and CCD images by using uniform square illumination (third row). Images in rows 1 and 2 were obtained by a WH inverse transform.

Download Full Size | PDF

As a first demonstration of the tomographic capability of the proposed system, three tomographic reconstructions by using an early gate of the time-resolved profile have been carried out. For each view, a single constant illumination pattern was used and detection was performed with $8\times 8$ WH patterns while the temporal gate has been chosen corresponding to the rising edge of the TR signal, here resulting in a time window of 500 ps length. The homogeneous measurements have been used to scale the inhomogeneous phantom data to match the magnitude of the forward temporal point-spread functions (TPSFs). The mesh used for the forward problem has 120,000 elements and 1016 TPSFs have been generated (127 WH patterns for eight views) and sampled in 156 temporal steps of 8 ps length. The computational time for the forward problem is about 10 s on a machine with 10 Dual Intel Xeon processors of 2.3 GHz each, an Nvidia Tesla K-40 GPU, and 64 GB of RAM memory. First, the eight TR measurements of the homogeneous phantom with planar illumination have been used to retrieve the background optical properties. For this purpose, a Levenberg–Marquardt fitting procedure with TOAST as a forward solver has been used, obtaining ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.88{\mathrm{mm}}^{-1}$ showing a good agreement with the properties measured using the spectroscopy system mentioned above. The reconstruction has been carried out on a regular grid of 85,731 points containing the whole cylindrical mesh with a voxel size of $0.5\mathrm{mm}\times 0.5\mathrm{mm}\times 1\mathrm{mm}$. The Jacobian has been calculated using the adjoint method [21] deploying the fast Fourier transform for fast implementation of the temporal convolution. The Gauss–Newton algorithm has been terminated after three iterations, after which the reconstruction ceased to improve. The overall reconstruction time was about 3 h.

Figures 4(a), 4(c), and 4(e) show the tomographic reconstructions of ${\mu }_a$ at different vertical slices. Due to the limited field of view of both illumination and detection, only a part of the cylinder can be reconstructed (about 16 mm centered at about 14 mm from the top). We observe good reconstruction quality pertaining to both the localization and relative contrast of the two inclusions. By fitting the reconstructed inclusions with a three-dimensional Gaussian function, we obtained a total contrast of about 0.3 times the truth for all cases. Moreover, in Figs. 4(b), 4(d), and 4(f), normalized line profiles across the inclusions at the $z=13\mathrm{mm}$ plane are shown for all the experiments. We observe a worsening of the localization for inclusion B of Exp 3 probably due to the reduced contrast produced by the low-absorbing solution poured in it. Finally, in order to quantify the localization capability of the reconstruction, the center of mass (COM) for each inclusion has been computed on a region twice larger than the inclusions (see Table 1).

Table 1. COM Coordinates in Millimeters

View Table

 

Fig. 4. Reconstruction of ${\mu }_a$ using an early gate. Slices are shown from $z=19\mathrm{mm}$ (top of the cylinder) to $z=4\mathrm{mm}$ with a step of 1 mm. The $z=0\mathrm{mm}$ plane is at the middle of the cylinder. Panels (a), (c), and (e) show reconstructions (values are in ${\mathrm{mm}}^{-1}$) for Exp 1, Exp 2, Exp 3, respectively. Panels (b), (d), and (f) show the normalized profiles of the absorption perturbations along line $\mathrm{\Gamma }$ connecting the centers of the two inclusions on the plane at $z=13\mathrm{mm}$.

Download Full Size | PDF

In conclusion, a fully tomographic time-resolved DOT system based on sampling in the spatial frequency domain (both illumination and detection spaces) and multiple view acquisition has been proposed and validated on a tissue-mimicking phantom, demonstrating state-of-the-art reconstruction quality. Moreover, the imaging capability of the system has been validated in CW by comparing SPC with standard CCD acquisition, showing the importance of the choice of illumination/detection patterns for imaging purposes. Future work will be devoted to the optimization of the data set (choice of illumination/detection patterns, number of views, and temporal gates) and system improvements (detection efficiency, calibration procedure) in order to strongly reduce the acquisition time, while preserving or even increasing the information content. In particular, adaptive basis scan approaches will be investigated [22].