Single-pixel camera photoacoustic tomography.

Since it was first demonstrated more than a decade ago, the single-pixel camera concept has been used in numerous applications in which it is necessary or advantageous to reduce the channel count, cost, or data volume. Here, three-dimensional (3-D), compressed-sensing photoacoustic tomography (PAT) is demonstrated experimentally using a single-pixel camera. A large area collimated laser beam is reflected from a planar Fabry–Pérot ultrasound sensor onto a digital micromirror device, which patterns the light using a scrambled Hadamard basis before it is collected into a single photodetector. In this way, inner products of the Hadamard patterns and the distribution of thickness changes of the FP sensor—induced by the photoacoustic waves—are recorded. The initial distribution of acoustic pressure giving rise to those photoacoustic waves is recovered directly from the measured signals using an accelerated proximal gradient-type algorithm to solve a model-based minimization with total variation regularization. Using this approach, it is shown that 3-D PAT of imaging phantoms can be obtained with compression rates as low as 10%. Compressed sensing approaches to photoacoustic imaging, such as this, have the potential to reduce the data acquisition time as well as the volume of data it is necessary to acquire, both of which are becoming increasingly important in the drive for faster imaging systems giving higher resolution images with larger fields of view.

Introduction

Photoacoustic tomography (PAT) is a hybrid imaging technique based on the use of laser-generated ultrasound within soft tissue that has been demonstrated in a wide variety of applications in preclinical research and clinical medicine., When a short pulse of near infrared (NIR) light is absorbed by chromophores within soft tissue, it gives rise to a pressure increase that propagates through the tissue as an ultrasound pulse and can be detected at the surface. Following the measurement of these photoacoustic signals on the tissue surface, an image of the initial pressure distribution can be reconstructed.

Photoacoustic signals are broadband, containing frequencies up to and above 10 MHz, at which the wavelength is . The classical approach to sampling—spatial and temporal—follows the Shannon–Nyquist theorem, which states that a band-limited signal can be recovered exactly if the sampling rate is at least twice the maximum frequency present in the signal. This suggests a large number of detectors are required for PAT, because the measurement plane should subtend a large solid angle at the imaging target in order to avoid limited-view artifacts, e.g., a aperture sampled at spacing results in more than 70,000 detection points. However, a dense array of many tens of thousands of small elements can be expensive and difficult to fabricate. One alternative is to use a smaller number of detectors and scan them across the surface, but this has the drawback of reducing the imaging frame rate. Another alternative is compressed sensing (CS), also called compressive sampling. The idea behind CS is that, under certain conditions, it is possible to reconstruct a target image accurately from fewer samples than determined by the Nyquist rate. The first requirement is that the target is known to be of low spatial (or more generally spatiotemporal) complexity, e.g., that it is sparse in a given basis. The second requirement is that the set of measurements contains nonredundant information about the target at all relevant scales. With these two requirements satisfied, it is often possible to reconstruct an accurate image of the target with only a fraction of the data that would be acquired with Nyquist sampling. The original and classic demonstration of this CS paradigm is the single-pixel camera.,

In this paper, a single-pixel camera is used to measure time-varying photoacoustic signals reaching a planar Fabry–Pérot (FP) ultrasound sensor, thereby facilitating compressed-sensing PAT. The FP sensor comprises a polymer film spacer sandwiched between a pair of dichroic mirrors. Any changes in the optical thickness of the spacer will modulate the optical reflectivity. When the wavelength of the interrogating laser light is tuned correctly, the reflected optical power will be proportional to the acoustic wave modulating the spacer thickness. A focussed laser beam is usually used to read out these acoustic pressure waves point-by-point. In this way, it is possible to synthesize arrays of many tens of thousands of detection points with small element sizes. This approach has been shown to give high-resolution images with high contrast. However, as mentioned above, the need to scan results in slow data acquisition. To acquire the data more quickly, one possibility is to use an interrogation system that can read out multiple points simultaneously, but it comes with high equipment cost and is technically challenging to implement. Another choice would be to use a camera, e.g., a CCD or CMOS camera, to record the signal at many pixels simultaneously,, but this requires a separate measurement for every time point. In this paper, the FP sensor was interrogated using wide-beam illumination and the reflected light was patterned on reflection from a digital micromirror device (DMD) before being collected into a single photodetector. With this system, it is possible to measure the spatial integral of the product of the pattern and the acoustic field at the sensor, i.e., projections of the field onto arbitrary spatial patterns. When choosing the set of patterns it is necessary to ensure that they will capture the data in such a way as to facilitate maximizing the image quality for a given size of dataset during the image reconstruction. Here scrambled Hadamard patterns are used. This approach has previously been demonstrated for real-time ultrasound detection; here, that work is extended to PAT.

Several experimental demonstrations of CS in PAT have been reported since Provost and Lesage first proposed it. The systems used have included ring arrays or circular scanning systems, systems employing integrating line detectors, linear arrays or line scans, and two-dimensional (2-D) arrays. In all of these studies, the sensors are restricted to subsampling the acoustic field at a set of points (or lines in the case of the integrating line detectors). In contrast, this paper is concerned with measurements made with a 2-D planar sensor interrogated with patterns. As well as a variety of detection schemes, a number of approaches to reconstruction from sparse data have been proposed for PAT: principal component analysis, sparsifying transforms,, deep learning,,,,, and variational approaches that minimize a functional,,,,, including joint motion estimation. Here a variational minimization approach will be taken. (For clarity, the term “compressed sensing” has been used in the photoacoustic imaging literature to refer to 2-D photoacoustic imaging using patterned excitation light. This is difficult to extend to 3-D imaging and is quite a different idea from the patterned acoustic sensing described here.)

Compressed Sensing Photoacoustic Tomography with Patterned Detection

In a PAT experiment, a laser pulse is used to illuminate the target volume, and where the light is absorbed it generates an initial acoustic pressure distribution . The aim is to image this initial acoustic pressure distribution. Because tissue is elastic, the initial pressure distribution excites acoustic waves, which propagate through the tissue to the sensor on the tissue surface. The acoustic field at the FP sensor at time can be denoted by , where is the acoustic propagation operator. To describe how the FP sensor facilitates CS, it is useful to start by considering a point-by-point interrogation scheme, in which there are detection points on the sensor. A complete dataset, , which can be described as a collection of measurements at different times , where represents the measurements at the detection points at a single time , and is the scalar acoustic pressure amplitude at a single point with coordinates at time . If, for a single time , the pressure field can be represented sparsely in a basis , then we can write , or where the columns of are basis functions and are the corresponding coefficients. If measurements are recorded at all points, as they are with a complete point-by-point scan, it would be possible to calculate the sparse coefficients using an inner product . This is analogous to the case of image compression. In the CS approach, the coefficients are obtained directly from measurements. Each measurement in the pattern-interrogation scheme is an integral of the field weighted by a pattern. In other words, the measurements are the set of amplitudes given by where each is a measurement pattern. The idea behind CS is to assume that is sparsely represented in a basis and use measurement patterns incoherent to it. The incoherence of the basis and the measurement patterns is crucial. If one knew, ahead of time, which sparse coefficients represent the solution, it would be possible to coherently measure the relevant coefficients directly using the respective basis functions as the measurement patterns, . In practice, however, the basis in which the data will be sparse is rarely known in advance, so CS proceeds using patterns that, through their incoherence, “equally” sense all the basis vectors and afterward use sparse recovery to extract the sparse coefficients from these measurements.

Once the measurements at each time step have been recorded, the challenge is to reconstruct the initial acoustic pressure distribution . If a full set of amplitudes have been recorded, one can obtain the original acoustic field at the FP sensor as and use standard PAT image reconstruction techniques for this scanning geometry, e.g., time-reversal. However, when only a subset of the data has been recorded, , there are two options: two-step schemes first reconstruct from compressed measurements, solving a problem akin to basis pursuit in CS (see e.g., Ref. 27, which assumes sparsity of in a Curvelet basis), and then use a standard photoacoustic inversion technique, whereas one-step schemes reconstruct the image directly from the compressed data. Reference 29 contains a detailed description of the one-step approach used here, namely an accelerated proximal gradient-type algorithm to solve the following minimization problem: where denotes the total-variation regularization. -Wave was used to compute both the acoustic operator and its adjoint. The regularization parameter was chosen via manual inspection. In the experiments reported here, the measurement matrix was chosen as a scrambled Hadamard matrix. This choice is both theoretically and practically appealing as explained below. The Hadamard transform is a matrix that can be recursively defined as and (note that this means that a Hadamard matrix can be written as a matrix with entries that are 1 or , multiplied by a normalization factor). The Hadamard transform is orthogonal and self-inverse, i.e., . A scrambled Hadamard matrix is formed by permuting the columns and rows of the matrix , to give , where and are the column and row permutation matrices, respectively. Notice that the scrambled Hadamard transform amounts to applying the column permutation to the vector , , applying the Hadamard transform to the permuted vector and subsequently permuting the rows. Thus, the cost of the scrambled Hadamard transform is of the same order as that of the Hadamard transform, which in turn can be realized with a fast algorithm, . Similarly, the inverse scrambled Hadamard transform is equivalent to applying inverse row permutation to , , applying the Hadamard transform to the permuted vector and subsequently applying the column permutation to the vector. Each row of the matrix represents a measurement pattern, so for CS, where , the first rows out of are selected. The selection of rows yields an underdetermined matrix with desirable properties for CS, similar to those of a random Gaussian/Bernoulli matrices.

Experimental Measurements

The experimental setup is shown in Fig. 1(a). The FP sensor (aluminum coatings, PPXC spacer, thickness of 40 μm) was illuminated by a 20-mm diameter expanded interrogation beam (Santec TSL-510 tunable laser source connected to a IPG Laser GmbH Erbium Fiber Amplifier EAD-4-L), Fig. 1(b). The reflected beam from the sensor was then redirected to the DMD (ViaLUX V-7000 DLP 0.7″ XGA array, pitch size of 13.68 μm) by the polarized beam splitter (Thorlabs CM1-PBS254). The optical reflected beam from the DMD was collected by lens () and focused into a photodetector (InGaAs Hamamatsu G8376-03 and a customized transimpedance amplifier configuration with DC- and AC-coupled outputs). A digitiser (NI PCI-5114) was used to acquire time series signals from the photodetector’s output. The system’s sampling rate and bandwidth were set at 50 and 20 MHz. The excitation laser was a -switched fiber-coupled Nd:YAG laser (Continuum Minilite-20, repetition rate of 20 Hz). The excitation beam diameter was and the total pulse energy was set to . Since the FP sensor used in this study was not transparent, the system was operating in forward mode.

Fig. 1

(a) Experimental setup: L1, lens 1; L2, lens 2; LP, linear polarizer; PBS, polarizing beam splitter; , quarter waveplate; PD, photodetector; DMD, digital micromirror device; (b) Fabry–Perot sensor with a wide interrogation beam; and (c) illustrations of scrambled Hadamard patterns.

Due to the periodic distribution of the micromirrors, a DMD acts as a 2-D diffraction grating when laser light is reflected from it. The reflected light from the FP sensor is, therefore, diffracted at the DMD into many diffraction orders that, within the NIR range, are well separated. The detector is positioned so that it collects light from only the strongest order. The details of how to optimize the DMD arrangement have been discussed elsewhere. In this setup, the strongest order at 1580 nm is at about 50 deg with the incident angle of about 26 deg. With this arrangement, the DMD was used to pattern (spatially sample) the reflected beam from the sensor, before it was passed through an integrating lens and onto the photodetector. The scrambled Hadamard patterns, Fig. 1(c), were sequentially displayed on the DMD, and a time series signal was recorded at the photodetector for patterns. Because light intensities cannot be negative, it was not possible to implement Hadamard patterns experimentally so (0, 1) patterns were used and the mean value of all the time series was subtracted from the set. Also, to avoid saturating the photodiode, the data corresponding to the all-1 pattern were constructed from two half-1 half-0 patterns. The active area on the DMD was chosen to be (). It was positioned to align with the most uniform region on the FP sensor so only one wavelength was required to interrogate the sensor. A scrambled Hadamard operator for was used for these experiments. Each image pixel corresponded to a group of micromirrors on the DMD, so the effective sensing element size was . The DMD and the digitizer were synchronized and triggered by the -switched laser.

The two phantoms shown in Fig. 2, a knotted artificial hair and a twisted polymer ribbon, were used in this study. The phantoms were immersed in 1% Intralipid with a reduced optical scattering coefficient of and were positioned above the sensor and 4 mm below the Intralipid surface. The diameter of the hair was about and the ribbon width was .

Fig. 2

Phantoms for photoacoustic imaging experiments (the black bar in both images measures 1 mm): (a) artificial hair ( in diameter) and (b) twisted black polymer ribbon ().

Results and Discussion

Photoacoustic images of the phantoms were recovered using Eq. (3) for different degrees of compression, and the results are presented in Figs. 3 and 4. Figure 3 shows the 3-D images of the knotted hair. The images were obtained using 100%, 50%, 20%, and 10% of the scrambled Hadamard patterns. It is instructive that with the omission of 50% of the data it is still possible to recover an image of similar quality to that recovered using the full data, and the main features of the targets are recovered successfully even when only 10% of the data are used. Figure 4 shows the slice images of the polymer ribbon.

Fig. 3

CS-PAT reconstructions () of the knotted hair with different levels of compression, visualized as maximum intensity projections from the top, i.e., (top row) through the sensor plane and (bottom row) from the side. The sensor plane is located on the top of the side view; the first six depth slices below it are set to zero to prevent sensor noise from dominating the maximum intensity projection from the top.

Fig. 4

CS-PAT reconstructions of the twisted polymer ribbon with different levels of compression visualized by a single slice through the ribbon, orthogonal to the sensor plane (located on top of the slice view).

In the accelerated proximal gradient-type algorithm used here for the reconstructions, the loss of data was compensated for through the use of the total variation regularization. However, there are many other minimization algorithms and regularization strategies described in the literature that could be employed to tackle this CS inverse problem. In particular, because the scrambled Hadamard basis is close to ideal for CS, the data obtained with this system allows an investigation into which basis is the best for reconstructions with partial data. It remains a focus of future work to determine which approach is optimal for photoacoustic imaging systems such as this.

Despite the advantages a CS approach has over full-data scanning systems, it still requires multiple sequential measurements. The data acquisition speed is, therefore, currently limited by the pulse repetition rate of the photoacoustic excitation laser (although the pulse repetition rates of photoacoustic lasers are gradually increasing). Devices using multichannel systems with arrays of detectors that can detect simultaneously do not suffer from this limitation. However, they are usually limited in their bandwidth and the cost and complexity typically becomes prohibitive beyond a thousand or so channels.

The images obtained here are not as high quality as published images obtained using a FP sensor with point-wise interrogation. This is due to a lower signal-to-noise ratio. There are principally two factors that affect the image quality for both patterned and point-wise detection: (1) the signal-to-noise ratio, and (2) the effectiveness of the CS, i.e., the degree to which the undersampling can be ameliorated in the reconstruction. It is the second point that is fundamental here. The signal-to-noise ratio, while important practically, is not a fundamental constraint, and it could in future be enhanced, for example, with improved sensor fabrication, averaging, or using higher power illumination and a photodiode with a greater dynamic range.

The sensitivity of FP sensors that are interrogated point-by-point can be maximized using high-finesse cavities and by tuning the wavelength to the optimal bias point on the interference fringe at each point. With the wide-beam, single-wavelength, illumination that this system requires, it is not possible to tune the interrogation wavelength to the optimal wavelength for each point, and so for this proof-of-principle experiment, a low-finesse cavity was used to ensure that every point in the field of view could be interrogated with some, if not optimal, sensitivity at the one wavelength. FP sensors with more uniform thickness over cm-sized areas are being developed, which will allow higher finesse cavities to be used in future, and therefore, larger signal-to-noise ratios to be achieved with this system.

Conclusions

3-D CS PAT using a single-pixel camera has been demonstrated experimentally. Scrambled Hadamard patterns from a DMD were used to sample the photoacoustic field as detected by a planar FP ultrasound sensor. Photoacoustic images of knot and ribbon phantoms were obtained with compression rates as low as 10%. CS will become important in PAT as the demand grows for fast, high-resolution systems with large fields of view.