1.

Introduction

Optical imaging in turbid media saw rapid development over the last few years. Applications in cancer detection,^{1, 2} the monitoring of brain activation,^{3, 4} and small animal fluorescence imaging^{5, 6, 7} were demonstrated. These developments also led to new commercial offerings. With this emergence of diffuse optical tomography (DOT), systems with an increasing number of sources and detectors were built to enhance image quality by exploiting a larger set of measures. In particular, it has been recently recognized that with dense spatial measurements, it is possible to image complex structures with high-resolution.⁸

In this work, we focus on noncontact transmission geometry with planar source and detector planes. Such geometry is commonly employed in optical mammography⁹ and small animal imaging associated with coupling liquid to mitigate high dynamical range in the detection signals.¹⁰ Previous diffuse optical tomography systems were implemented in this geometry by either positioning a series of punctual light sources and detection points at different positions on the sample surface or scanning the sources. This scan was done mechanically⁹ or optically with galvanometer mirrors⁸ providing high-density measurements. Detection in these configurations was mostly performed with either single detectors or a sensitive camera. While these acquisition strategies are experimentally straightforward to implement, they suffer from a major drawback. To acquire dense illumination and detection measurements, every light position has to be spatially scanned or time multiplexed, thereby limiting the frame rate that can be achieved. This is especially crippling in functional imaging applications where spectral information at multiple wavelengths is required.¹¹ Furthermore, the large amount of data acquired leads to massive inverse problems and necessitates extensive computation resources.

The incentive for this work is found in two simple observations: (1) even if high-density measurements were used to reconstruct complex structures in Ref. 8, only a smaller number of Fourier modes were employed for the reconstructions; (2) by using digital micromirror devices (DMDs) not only for illumination but also for detection, it is possible to generate specifically those relevant measures for reconstruction. With this approach, the number of measurements required is drastically reduced, leading to significant data compression and potential frame rate increase.

The ability to illuminate specific spatial Fourier modes was exploited in recent work under the denomination of modulated imaging.¹² By exploiting new technologies—more precisely, DMD arrays—it is possible to shine arbitrary patterns of light on tissues and use a model of light propagation to infer information about the media. For example, baseline tissue optical properties were estimated by illuminating distinct spatial Fourier modes and using diffusion theory.¹³ A recent extension included the study of time-domain spatial mode propagation¹⁴ and two-layer optical parameter estimation.¹² All these initial works were not tomographic but “sectioning” methods, as the tomographic information was encoded in the frequency and not retrieved through an inverse problem. Only recently has tomography been explored.¹⁵ Furthermore, “pure tones” or spatial frequencies were used for illumination, diminishing significantly the acquisition rate. The technical reason for the reduced acquisition rates is that mirrors in DMDs can only be switched “on” or “off” and gray levels are produced by applying a pulse-width modulation (PWM) of the mirrors to create the intensity effect. While it is estimated that with current micromirror devices, the whole array can be controlled at rates exceeding $10\mathrm{kHz}$ , the frame rate in PWM mode is decreased to around $60\mathrm{Hz}$ .

In these previous applications, DMD technology was used only for illumination in combination with a camera for detection. The wide-field acquisition capabilities of the camera allowed the computation of all detection Fourier modes with a single frame. However, acquisition frame rates are limited by the integration time of the camera—typically, $100\mathrm{ms}$ —due to the low-light nature of diffuse optical applications. Conversely, it is possible to acquire sequentially all the detection Fourier modes directly with a detection channel comprised of a DMD and a single detector such as a photomultiplier tube (PMT) or an avalanche photodiode (APD). Thanks to the increased level of light, measures can be acquired at increased frame rates to obtain a complete set of measures for high-resolution DOT at unprecedented speeds.

In this work, a dual DMD acquisition technique to perform diffuse optical imaging is experimentally demonstrated for the first time in phantoms. Previous work on diffusion theory is used to guide an approximation of the Fourier measures by simple “on-off” mirror patterns, thereby increasing the DMD switching speed (from $60\mathrm{Hz}\text{to}10\mathrm{kHz}$ ). We show that by selecting a subset of measurements, it is possible to get sufficient data for high-fidelity three-dimensional (3-D) image reconstruction at real-time frame rates.

2.

Model

In the diffusion approximation to the radiative transport equation, the propagation of light in tissues (valid if ${\mu }_a⪡{\mu }_s^{\prime }$ ) is modeled by a diffusion equation for the fluence $\varphi$ :^{16, 17, 18}

By assuming translational symmetry in the $x\text{-}y$ directions, define the Fourier transform in the spatial coordinates that are “in-plane,” i.e., ${\rho }_{s,d.}$ , as

As mentioned in Sec. 1, a few issues arise with this strategy when focusing on this equation. First, the pure Fourier measures are slower to generate with DMDs due to the necessity of generating gray levels. Until now in modulated imaging, sinusoidal measures were used, but these carry a temporal cost due to the PWM of some mirrors and integration time of the camera.¹⁵ Second, the optimal spatial patterns described by Eq. 7 rely on translational invariance, which in some cases will not be valid. Last, in situations where the diffusion theory does not hold, this approach is not adequate since the forward model used is inaccurate.

However, these issues can be mitigated by the following observations: While the diffusion theory and spatial invariance are not always true, it is possible to build a forward model that takes into account those two issues. For example, a Monte Carlo simulation of the fluence induced by a pattern is simple to compute by convolution, taking into account both the issue of spatial invariance (by modeling an arbitrary medium) and the diffusion theory approximation. While the analytical inversion formula requires pure Fourier modes that are longer to measure, the DMDs can be used to get measures that are approximates of the Fourier modes. Formula 7 is no longer available in that case, but it provides guidance as to which measures are important in the inversion. The decrease of the number of measures (by using approximate Fourier modes) can then be exploited by using a pseudo-inverse (that can be precomputed) to provide a real-time inversion.

These issues are investigated in the following, both in simulations and with experimental work on phantoms. We first use Monte Carlo simulations to generate measures by convolution in a simple homogeneous volume and use perturbation theory to compare different illumination/detection schemes for tomography. Once identified, we study the effect of noise on the reconstructions in a given scheme. Last, we provide experimental evidence of the feasibility of the approach proposed in this paper by using liquid phantoms.

3.

Simulations

The availability of a large number of illumination and detection patterns ( $2\times {10}^{12}$ combinations given that each array has ${10}^6$ mirrors) makes it impossible to explore the entire measurement space. Since our goal was to provide a proof of principle for this dual DMD technique, we considered and compared two types of patterns. Many others can be explored, but doing so requires an extensive study. To model light propagation, Monte Carlo simulations were performed for a single source in a homogeneous medium.²⁰ Assuming perpendicular translation invariance, the resulting fluence was then convolved with each illumination pattern to create the fluence distribution for each case. A sensitivity matrix was computed by doing the same for the detector configuration.

Using this forward model, a numerical phantom was created and simulated measures generated. These data were used for the ensuing inversion, which was done with a regularized pseudo-inverse. Forward and inverse problems were performed on distinct grids to minimize the “inverse crime.” Figure 1 shows the two types of patterns used in simulations. In each case, the same pattern sets were used for both illumination and detection ( $N=36\times 36$ patterns). The total number was chosen to achieve a frame-rate of $2\mathrm{Hz}$ , which was the initial target but can be chosen arbitrarily. In simulations, a numerical phantom of rectangular geometry of $9\times 7\times 2\mathrm{cm}$ was used. Optical parameters were set to ${\mu }_a=0.01∕\mathrm{mm}$ and ${\mu }_s^{\prime }=1∕\mathrm{mm}$ to mimic the experimental configuration described in the following. Two numerical inclusions were added to the phantom consisting of two bars, each of $5\mathrm{mm}$ width and situated at $z=8\mathrm{mm}$ . The contrast of the inclusions was set to 100%, i.e., ${\mu }_a=0.02∕\mathrm{mm}$ .

Fig. 1

Top left: Phantom used for the simulations. Top right: Two types of patterns for illumination and detection. Bottom: Reconstructions for each pattern type. Top middle: Cross sections and quantitative values of the differential absorption across a $z$ -slice. The volume for the reconstruction is divided in 16 $z$ -slices, each having $1.25\mathrm{mm}$ thickness; even $z$ -slices are shown.

Figure 1 compares reconstructions between the two illumination/detection patterns: checkerboard and moving bars. Taking a line profile across the $x$ axis leads to a FWHM of $9\mathrm{mm}$ for both bars, but the contrast for pattern 1 was 50% lower than that of pattern 2 (moving bars). For inclusions having 100% contrasts, the reconstructed value was 40%, coherent with the observed spatial extent of the solution. Expectations were that the checkerboard pattern would perform better based on the analytic block diagonalization of Eq. 7 and a singular value analysis (Fig. 2 ). However, it was the moving bars that enabled the best reconstructions for the chosen phantom. The numerical noise introduced by the discretization and the associated truncation performed during the computation of the pseudo-inverse may explain the slightly better performance by the moving bars pattern.

Fig. 2

Singular values for the checkerboard (circles) and moving bars patterns (crosses).

To evaluate sensitivity to noise, simulations were performed with 0, 5, and 10% noise added to the simulated measures with respect to the mean. In each case, the regularization parameter was optimized based on the mean-squared distance between the reconstruction and the gold-standard. Results are presented in Fig. 3 for the moving bars (pattern 2).

Fig. 3

Simulations with 0, 5, and 10% noise around the mean measurements in the frame. Line profile of the reconstructions across the bars at $z=8\mathrm{mm}$ .

While spatial aliasing was present in all reconstructions, degradation of the reconstruction appeared only at high noise levels (10%). These simulations were used to find the optimal regularization parameter for noise levels similar to that measured experimentally (1 to 2%; see the following). These latter values were used for the following experimental reconstructions.

In all of the preceding simulations, the reconstructions were performed both on an eight-core computer ( $32\text{-}\mathrm{Gb}$ RAM, $2.83\mathrm{GHz}$ ) and a graphical processor unit (GPU; GForce 280GTX $1\mathrm{Gb}$ ) system to evaluate the feasibility of real-time 3-D imaging. On both computational platforms, reconstructions based on precomputed pseudo-inverse were executed in a few milliseconds. For the inverse problem investigated herein, 1296 pattern pairs and 131,072 voxels, the multicore and the GPU platform led to $540\text{-}\mathrm{ms}$ and $220\text{-}\mathrm{ms}$ inversion time, respectively.

4.

Experiments

The experimental setup is depicted in Fig. 4 . The illumination consisted of a single LED ( $635\mathrm{nm}$ , $400\mathrm{mW}$ ) collimated with an $\mathrm{f}=75\text{-}\mathrm{mm}$ biconvex lens. The collimated beam was aligned with the first micromirror matrix (DMD, TI Discovery 3000 board), and the pattern generated by the DMD was then imaged on the surface of the object through a $75\text{-}\mathrm{mm}$ focal length magnifying lens. The average surface illuminated was $100{\mathrm{cm}}^2$ , leading to an estimated $4\mathrm{mW}∕{\mathrm{cm}}^2$ fluence, well below ANSI standards. The detection side was symmetrical: the light coming out of the turbid medium was imaged on a second DMD with another $f=75\text{-}\mathrm{mm}$ biconvex lens. The detection patterns generated by the second mirrors then redirected the light in the direction of a collimating lens adjusted to image the width of the DMD onto the surface of a photomultiplier tube. Data from the PMT, as well as the trigger signals controlling the DMD boards, went through a homemade adaptive electronic circuit to a DAQ card (National Instruments PCI6009; transistor-transistor logic (TTL) levels differ between the DAQ and the DMD boards).

Fig. 4

Experimental setup.

For the experiments, moving bars were chosen based on the simulations results. In the following pilot study, a complete acquisition frame consisted of 1296 measurements (36 emission patterns by 36 detection patterns), as in the preceding simulations. This number is not limited, however; it is possible to include larger numbers of source detector pairs by reprogramming the DMD boards. DMDs took approximately $125\mu \mathrm{s}$ to generate a pattern. Once the pattern was generated, 10 measurements from the voltage outputs of the PMT were taken by the acquisition board, sampling at $200\mathrm{k}\mathrm{S}∕\mathrm{s}$ , adding a time lapse of $50\mu \mathrm{s}$ . Last, software averages and communication between computer and DMD boards took $210\mu \mathrm{s}$ , for a total of approximately $385\mu \mathrm{s}$ by pattern. Thus, a full frame of 1296 measures was acquired in $0.3\mathrm{s}$ leading to a frame rate of $2\mathrm{Hz}$ including overhead. Since all patterns were known ahead of time, a precomputed pseudo-inverse with Tikhonov regularization was used to reconstruct the images in real time. The data was sent by TCP/IP from the acquisition computer to a second reconstruction computer with little transmit time given the intrinsic data compression afforded by this illumination-detection approach. Data was reconstructed in $220\mathrm{ms}$ and displayed in real-time on this second computer.

The experimental acquisitions were done on a controlled phantom with optical properties set to ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and a ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ (dilution of Intralipid 10% with India ink). An imaging chamber with Plexiglas windows $(150\mathrm{mm}\times 100\mathrm{mm}\times 20\mathrm{mm})$ was used to hold the liquid. Since the reconstruction algorithm required a baseline measure, a set of frames was taken on the homogeneous liquid, phantom. Following this, the inclusion was inserted in the liquid, and further frames were taken. In all cases, data were measured at the $2\text{-}\mathrm{Hz}$ frame rate and stability of the system was first assessed. The mean noise level of the measurements was observed to be between 1 to 2% for a $2\text{-}\mathrm{Hz}$ frame rate. Then, the system was run continuously following an initial frame, and no drift was observed while maintaining the $2\text{-}\mathrm{Hz}$ rate over periods exceeding $30\min$ . Figure 5 shows the baseline and phantom frames with their noise.

Fig. 5

Top: Sequence of patterns used for both illumination and detection. Middle: Measured data $(N=36\times 36)$ for baseline $\left(I_0\right)$ and phantom $\left(I\right)$ with noise in the same graph $(<2\mathrm{\%})$ in both cases. Bottom: Schematic of the rod inclusion in the liquid phantom for the reconstructions.

Using a second DMD for detection seems to provide faster acquisitions than previously obtained in the literature. In previous work using a CCD,²⁰ $100\mathrm{ms}$ was necessary to integrate light information in similar conditions as the ones presented here. Given that only 36 illumination patterns would be necessary in this case, since detection patterns would be computed in software, the estimated acquisition time for an equivalent data set is $3.6\mathrm{s}$ . This is approximately 7 times slower than the method presented in this study but is also dependent on the camera and optical configuration.

Figure 6 shows the reconstructions for distinct inclusion positions and a comparison of these reconstructions with data simulated for the same phantom positions. In this case, the phantom consisted of a $3\text{-}\mathrm{mm}$ -diam graphite rod (completely absorbing). For the three positions, the FWHM of the reconstructed inclusion was estimated to be $6.25\mathrm{mm}$ for the ones positioned near the boundaries and $7.5\mathrm{mm}$ for the inclusion located at the center of the liquid. These results are consistent with a simulation done numerically for the same inclusion size and the associated reconstruction (similar noise levels were used in the simulation). Contrast is observed to be similar in both simulations and experimental data, but it is difficult to infer any meanings to the values since the inclusion was completely absorbing and a Born approximation was used for reconstruction. A mild aliasing phenomenon is observed and is due to the regularity of the chosen measurement patterns for emission and detection. This aliasing is expected to disappear if a set of incoherent measures are used akin to work performed in compressed sensing in recent years.²¹ Resolution in both simulations and experiments are similar, as are the simulated and experimental measurement frames (data not shown).

Fig. 6

Three different phantoms used and the associated reconstructions with a regularized pseudo-inverse for the experimental data and a simulation modeling the phantom for comparison. The choice of barpatterns introduces a small aliasing in the images related to the positions of the measurements.

5.

Conclusion

In this work, simulations and experimental results were presented to validate the feasibility of a fast tomographic imager using structured illumination and detection. Simulations demonstrated the possibility to perform reconstructions with high resolution and accurate quantification without using a large number of measures. The reconstructions’ sensitivity to noise was also evaluated, and the acquisition scheme was shown to be stable in the presence of noise equivalent to that found in experiments. The use of a small number of measures enabled us to perform a full pseudo-inverse on a desktop computer. With precomputation of the pseudo-inverse, a frame rate of $2\mathrm{Hz}$ , including 3-D image reconstruction, was achieved on a GPU. These results were further exploited to provide ultra-fast reconstruction display to assess our experiments in real time. While not as accurate as the reconstructions found in, Ref. 8 the reconstructions were comparable to the literature in resolution. It is expected that by increasing the number of patterns used, the resolution of the reconstructions will increase and may be able to achieve that found in Ref. 8. However, the frame rate will diminish, and further work is required to assess the best patterns to optimize the resolution to frame rate trade-off. We believe that the ability of performing real-time 3-D images is the main benefit of this system, and applications to dynamical imaging of optically thick tissue such as pharmacokinetics will be the subject of further studies.²²

A more extensive characterization of the reconstructions could have been performed, but since only simple illumination-detection patterns were explored, it is not clear that it provides a quantitative evaluation of the system. The goal of this paper was to provide a proof of principle for the dual DMD approach and estimate achievable frame rates and noise levels. A separate study to identify the optimal illumination-detection strategy is needed, but it is beyond the scope of this work. With these optimal measures, the trade-off between frame rate and image quality can be evaluated more accurately. Another drawback of this work is the use of precomputed pseudo-inverses to perform image reconstructions. It is not clear that this is the optimal reconstruction method, but it does provide a means to visualize data rapidly.

Besides illuminating only the relevant modes, there are significant advantages to this new approach when compared to other techniques. First, by using wide-field illumination, it is possible to increase the light intensity while preserving the ANSI limit for tissue exposition since the limit is per area. Here data were acquired with a simple LED, but using a more powerful source, one should be able to perform fluorescence tomography. Second, although higher total light intensity is used, it is spatially distributed, and at any point in the volume imaged, the light intensity is lower than in local illumination schemes. This has the advantage of decreasing bleaching effects for potential fluorescence applications while still increasing signal-to-noise ratio (SNR). Last, compared to previous reports, when using a DMD for both illumination and detection, it is possible to limit the measurements to only those necessary and recollimate the light detected to spread it in distinct spectral components, providing spectral measurements without time penalty. We hope to revisit these issues in future work.