1.1.

Background

Photoacoustic Imaging (PAI) is a promising hybrid modality capable of providing image contrast by optical absorption, but with tissue penetration and spatial resolution at depths that are superior to pure optical imaging methods.¹ PAI has been employed in vivo to visualize blood-vessel structures,^{2, 3, 4} detect breast tumors,^{5, 6, 7} estimate oxygenation levels,^{8, 9, 10} perform functional imaging,^{3, 11} and track molecular probes.¹²

The technique is based on volume irradiation with a short laser pulse followed by preferential energy deposition at the optically absorbing structures within the volume. These structures heat up slightly, but so rapidly that the condition of thermal confinement is met and thermoelastic expansion leads to a transient bipolar pressure wave that propagates outward toward the surface.¹³ The frequency composition of the pressure transient (i.e., the “photoacoustic wave”) is typically in the ultrasonic range,¹⁴ and is detectable by wide-band ultrasound detector(s) placed in acoustic contact with the surface of the volume. Based on the time-domain surface measurements, the distribution of the photoacoustic (PA) sources (structures) inside the volume can be reconstructed using a back-projection algorithm, several of which have been suggested and adapted from medical imaging.^{15, 16, 17, 18, 19}

The photoacoustic process is intrinsically three-dimensional. This is a product of the optical radiation scattered in tissue, which diffusely illuminates a volume, rather than at a single point, line, or plane. Furthermore, the PA waves are of a spherical nature and propagate in all directions from their point of generation. Several approaches have been suggested for 3-D PAI and are differentiated by detection scheme. These can be divided into scanning methods, where a single detector is scanned along the detection surface in two dimensions,^{2, 3, 4, 20, 21} and staring methods, where a 2-D array of detectors is used and no scanning is necessary.^{6, 22, 23} A combined scanning–staring approach has been suggested as well, where a linear array or combination of single detectors is mechanically scanned.^{7, 24, 25, 26, 27}

1.2.

Motivation and Approach

We were motivated to develop a method for 3-D PAI that can be used for small animal research with eventual translation to clinical applications. We based the development of the imaging system on the following design goals: (1) real-time image acquisition capability, (2) easy placement and accessibility to the scanned object, (3) simple design (i.e., low detector count), and (4) potential for transportability to permit integration with other imaging modalities for multimodality hybrid imaging and validation [e.g., magnetic resonance imaging (MRI) and x-ray computed tomography (CT)]. Review of the literature suggested that the methods developed to date lack in one or more of these design features. Mechanically scanned approaches require long times for data acquisition, and therefore cannot be used in real time. Depending on the scanning geometry, some also lack accessibility to the subject under investigation and are too cumbersome to be considered candidates for transportability. They are, however, usually inexpensive and simple to set up, since only a single transducer and detection channel is used. Staring methods do not require mechanical scanning and hence, in principle, can acquire all the required data during a single laser pulse, which has the potential to create 3-D photoacoustic images in real time. Staring methods also have open access to the imaging region, and some have been built as transportable systems. ^{6, 7} However, all published implementations of the staring method have incorporated a dense array of detector elements, which is difficult to manufacture. Furthermore, real-time 3-D imaging with dense detector arrays is prohibitively expensive due to the high channel count required for simultaneous data acquisition from all detector elements. A combined scanning–staring approach can reduce system cost compared to a pure staring approach; however, imaging speed is still limited by the need to scan the detector array mechanically.

Based on the design goals and the implementations described in the literature, we chose to build a system based on a staring 2-D detector array due to its inherent speed and sample accessibility. In order to reduce channel count, a sparse detector array geometry was investigated, where instead of packing detection elements as close together as possible, the elements were spaced to cover a large area but with a low packing ratio. The spread of the detection elements over a large area was necessary to provide a wide range of viewing angles of the volume to be imaged, which is known to improve the boundary reconstruction of PA objects in the field of view.²⁸ The geometry of the sparse array was annular to accommodate backward-mode illumination of the object through the center of the annulus, which was amenable to easy object placement and accessibility.

Due to the small number of detectors, and the limited viewing angles they cover, a simple radial back-projection of the PA signals into volume space would have resulted in streaking image artifacts nearby the reconstructed location of the PA sources. For this reason, we chose to adapt an iterative image reconstruction algorithm, originally developed by Paltauf, ¹⁷ to handle data collected with the sparse detector array. In addition, we developed a characterization method²⁹ to measure the 3-D sensitivity distribution of the annular detector array, which was used to constrain the iterative image reconstruction algorithm.

2.1.

The 3-D Photoacoustic Imaging System

The 3-D photoacoustic imaging system [Fig. 1a ] used backward-mode illumination, where the laser beam illuminated the imaging volume from the same side as the ultrasound detectors. The ultrasound detectors were sparsely arranged around the laser entrance window and pointed toward the volume of interest. The electronic signal generated by the detectors was captured in parallel by a dedicated data acquisition (DAQ) system and was sent to the computer, where it was processed by an iterative image reconstruction algorithm based on back-projection. The imaging system was built from the hardware components described here.

Fig. 1

Overview of the sparse photoacoustic imaging system. (a) Imaging mode, where the object to be imaged was placed on top of and in acoustic contact with the detector array and the laser illumination was delivered through the window in the center. PA signals were acquired once by the DAQ system for all detectors and were sent to the computer for image reconstruction. (b) Setup for characterization, where a point photoacoustic source was scanned in a tank of water across the volume of interest and sensitivity and coverage maps were generated. PA signals were acquired in parallel by the DAQ system for all detectors at each scan point and were sent to the computer for analysis and display. (c) Top view (left) and cross-section side view (right) of the sparse annular detector array. The array consisted of $14\mathrm{PV}\mathrm{dF}$ film detectors, $3\mathrm{mm}$ in diameter, arranged on a $30\text{-}\mathrm{mm}$ -diam ring, an optical window in the center of the ring (for laser transmission), and an annular delay line with a wedge (shaded) directly on top of each detector. The delay line refracted the acoustic waves arriving at the detectors so that an effective focal zone was created above the center of the array, at a height of $30\mathrm{mm}$ (marked by the crosshair on the side view). The characterization scan volume was marked by the dashed box on the side view. The Cartesian coordinate system shown here was used consistently throughout the text. The array plane was located at $z=40\mathrm{mm}$ .

2.2.

Iterative Image Reconstruction Algorithm

We used an iterative image reconstruction algorithm first proposed by Paltauf ¹⁷ The main assumption of the model was that the PA pressure profile measured at each detector could be described as a linear superposition of individual PA pressure profiles generated by spherical sources located on a 3-D grid of voxels. In the forward model, each spherical source generates a pressure profile given by:

The velocity potential $\left(\mathit{VP}\right)$ is obtained by integrating the pressure profile with respect to $\tau$ and has the form:

For a given detector located at position ${\mathit{r}}_i$ , the $\mathit{VP}$ calculated for a PA source located at ${\mathit{r}}_j$ can be written as:

The back-projection model is based on the assumption that each $\mathit{VP}$ time point ${\mathit{VP}}_i\left(t_k\right)$ is a result of a linear superposition of sources that lie on a shell around detector $i$ whose distance from the detector is in the range $c\cdot t_k-R_s<r_{ij}<c\cdot t_k+R_s$ .

The amount of signal back-projected into each voxel $j$ is given by

Experimentally, the measured PA signals were rectified and filtered with a time-domain moving average to compute ${\mathit{VP}}_i^{\exp }\left(t\right)$ . A first estimate of the master image was formed by setting all voxel intensities to zero. The master image was forward-projected to obtain ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ . The difference ${\mathit{VP}}_i^{\exp }\left(t\right)-{\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ was calculated and back-projected into volume space. The resultant difference image was added to the master image, and the outcome was forward-projected to obtain the new ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ . The iterative process was repeated until either of two stop criteria was met: (1) a set number of iterations was reached, or (2) a relative change in ${\mathit{VP}}_i^{\exp }\left(t\right)-{\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ between successive iterations dropped below a given threshold.

The sensitivity map of the detector array (see the following) was interpolated and used as the weighting function $A_{ij}$ for the forward- and back-projections. This map attenuated $\mathit{VP}$ estimates from voxels where detector sensitivity was low and enhanced $\mathit{VP}$ estimates from voxels where detector sensitivity was high.

All images were reconstructed with the following algorithmic parameters: image volume of $20\mathrm{mm}\times 20\mathrm{mm}\times 20\mathrm{mm}$ , voxel size of $1\mathrm{mm}\times 1\mathrm{mm}\times 1\mathrm{mm}$ , Cartesian coordinate system as defined in Fig. 1c, and 500 iterations. The relative voxel intensities were displayed on a linear grey scale.

2.3.

Sensitivity Mapping of the Detector Array

Kruger have reported the use of a scanned PA point source for calibrating the spatial response of a linear detector array along two orthogonal axes.³⁰ However, to account for the response profile of each detector element and the overall sensitivity profile of the array in the entire imaging volume, we raster-scanned a PA point source in three dimensions [Fig. 1b] to generate volume maps of detector sensitivity and coverage.²⁹ The point source was constructed from a $400\text{-}\mu \mathrm{m}$ -core optical fiber whose tip was polished and coated with black paint. At each point of the scan, the source was motionless and the laser was fired through the fiber, which resulted in a pressure transient generated at the coated tip due to the photoacoustic effect. Each element of the detector array received the outward-propagating pressure wave with a delay time and viewing angle that depended on the source–element separation geometry. Each detector element produced an analog voltage proportional to the pressure profile, which was amplified, digitized, and stored. For each laser pulse, part of the beam was directed toward a photodiode and recorded synchronously with the PA signals to monitor fluctuations in laser power.

The scan was performed in water for acoustic coupling between the source and the detector array. The volume under investigation was a $25\mathrm{mm}\times 25\mathrm{mm}\times 25\mathrm{mm}$ cube, extending from $(x,y,z)=(0,0,0)$ to $(x,y,z)=(25,25,25)$ . The volume was chosen to cover the area of the optical window in the $xy$ plane and to cover the focal zone of the annular array in the $z$ direction [see Fig. 1c]. The acquisition points were spaced $5\mathrm{mm}$ apart in the $x$ , $y$ , and $z$ directions for a total of 216 points (i.e., $6^3$ ) that spanned the volume. The data sets acquired from the array elements at each scan point were fed into the analysis software, where the location and strength of each PA peak was detected, and the effects of laser power variation and acoustic spherical divergence were taken into account. In a separate experiment, the angular emission profile of the point source was measured by rotating the source at a fixed point above one array element and recording the relative change in PA signal. The resulting curve was interpolated and incorporated into the calculations of the sensitivity maps. The software generated a sensitivity map for the detector array that depicted the relative strength of the PA signal detected by each element at each scan point. The sensitivity map was used as a weighting function by the image reconstruction algorithm.

2.4.

Element Coverage Mapping of the Detector Array

The analysis software also generated a coverage map, which depicted the number of elements in the array that received signals above background noise for each voxel. Due to the limited angular response of the detectors and the sparse geometry of the array, it was expected that there would be voxels where not all detectors had enough sensitivity to detect the point source. Visualization of the regions of maximum coverage gave an indication of where the best imaging quality was to be expected, since for voxels where coverage was good, the reconstruction algorithm would be able to incorporate data from a multitude of angles.

2.5.

Phantom Development

To validate the accuracy of the 3-D image reconstruction, we developed a 3-D phantom that had a distribution of PA sources with known spatial locations. Instead of building a physical phantom of small absorbing objects, we used a composite of the PA signals collected during the sensitivity scan to build up a synthetic phantom of point sources. This approach had two advantages: (1) the position of each source was accurately known and referenced to the same coordinate system as the imaging volume, and (2) the PA signal intensity was constant regardless of position, since the source was illuminated through the optical fiber and the signals were corrected for variations in laser energy. This approach provided a means to generate a variety of synthetic phantoms with user-defined distributions of point sources. Although the synthetic sources were accurately localized and uniform in intensity and therefore well suited as test objects for the 3-D PAI system, they did not represent a real imaging condition with external illumination and simultaneous detection of multiple point absorbers. To test 3-D PAI on real objects with external illumination, we first used the coated fiber tip as a point absorber and illuminated it in backward mode through the laser access window instead of through the fiber directly [Fig. 1a]. We then imaged an extended object in the form of a $0.9\text{-}\mathrm{mm}$ -diameter graphite rod, which was suspended above the center of the detector ring at a height of $30\mathrm{mm}$ . All phantom experiments were carried out with distilled water as the acoustic coupling medium.

3.1.

Characterization of the Sparse Array

Figure 2a and 2b show typical PA signals acquired at one of the characterization scan points. Figure 2c represents the calculated coverage map. The coverage peaked at the center of each $xy$ plane, as expected from symmetry, but there was also a focusing effect along the $z$ axis. The largest lateral extent was observed in the plane $z=10\mathrm{mm}$ , and the coverage diminished toward both higher and lower $z$ values. The coverage distribution in 3-D, hence, took a form similar to an ellipsoid. The center of the ellipsoid was located directly above the array center and at a distance of $30\mathrm{mm}$ from the plane of detection. Its short axis, oriented laterally, and long axis, oriented along the $z$ axis, were 15 and $30\mathrm{mm}$ long, respectively. Figure 2d shows the calculated sensitivity map, which describes the relative sensitivity of each array element to each voxel in the volume under investigation. The map displays the $xy$ distribution for each $z$ plane (map rows) and for each detector element (map columns). In the $z$ direction, the combined effect of the delay line and the detector response profile can be appreciated: starting at the bottom row (i.e., the highest $z$ value and the $xy$ plane closest to the detectors), the most compact distribution is observed. Moving away from the detectors, the distribution becomes extended in the $xy$ plane. The location of the center of the distribution, however, also shifts away from the detector location for increasing distances from the detector plane. For example, at $z=25\mathrm{mm}$ , the distribution for detector 13 peaked at the bottom-left corner, which was closest to the detector location; at $z=0\mathrm{mm}$ , it peaked on the opposite side, at the top-right corner. These 3-D sensitivity distributions can be explained by the angular response profile of each detector/wedge combination. We can imagine each detector’s sensitivity field as a beam of light (e.g., from a spotlight) and correlate the beam divergence with the angular profile’s FWHM and the pointing angle with the refracted angle of the delay line. It is then easy to visualize that when the beam enters the imaged volume, it is compact, and when it leaves at the farthest $z$ plane, it is broad. Now, thinking of a multitude of light beams, it is evident that the region where they meet (the region of maximum coverage) will have the highest intensity and an ellipsoidal shape.

Fig. 2

Characterization of the detector array. (a) Typical PA signals obtained from the 14 array detectors and photodiode when the PA source was placed at $(x,y,z)=(10,15,10)$ . The signal from the photodiode appears as a sharp negative peak at time zero, and the 14 PA signals, arriving at the detectors at a delay time of $20\text{to}25\mu \mathrm{s}$ , are outlined by the dashed box. (b) A magnified view of the region of interest outlined in panel (a). Three detector channels are highlighted for better visualization (channel 1—thick solid line; channel 3—solid triangles; and channel 13—open circles). Signals from the same three detector channels are presented in Figs. 3, 4, 5. (c) Volume coverage map where each pixel represents a point from the 3-D characterization scan. The 3-D volume is presented as separate 2-D slices along the $z$ axis. The gray scale represents the number of detectors for which there was sufficient detection sensitivity, on the range of 3 to 14 (all detectors). (d) Sensitivity map where each image contains the sensitivity distribution in the $xy$ plane, for the indicated $z$ plane and detector number. The gray scale represents relative units, with white being most sensitive. All dimensions reported in mm.

Fig. 3

Reconstruction of a point source located at $(x,y,z)=(10,15,10)$ . In panels (a) to (c), readings from a point in the calibration scan were used to synthesize the source. In panels (d) to (f), the coated fiber tip was illuminated from below in backward mode. (a) Measured velocity potentials (VPs). Only three detector channels are presented for clarity. See Fig. 2b caption for details. (b) Estimated VPs from the reconstruction algorithm. The three curves presented correspond to the same channels as in panel (a). (c) Three orthogonal slices of the reconstructed volume, showing the location of the reconstructed point source. Note that the source was brightest in the $x=9\mathrm{mm}$ plane. (d) Measured VPs for the source illuminated in backward mode. Note the “missing” signal on channel 1 indicated by the solid curve. (e) Estimated VPs. Note that the signal “reappeared” on channel 1. (f) Corresponding volume slices of the reconstructed image. The gray scale indicates relative voxel intensity on the range of [0,1]. All dimensions reported in mm. Axes tick marks are $5\mathrm{mm}$ apart. Arrows indicate axis origin and direction.

Fig. 4

(a) to (c) Reconstruction of three point sources arranged vertically on the $z$ axis: $(x,y,z)=(10,15,5)$ , (10,15,10), and (10,15,15). Note that the sources were brightest in the $x=9\mathrm{mm}$ plane. (d) to (f) Reconstruction of three point sources arranged horizontally on the $y$ axis: $(x,y,z)=(10,5,10)$ , (10,10,10), and (10,15,10). The gray scale indicates relative voxel intensity on the range of [0,1]. All dimensions reported in mm. Axes tick marks are $5\mathrm{mm}$ apart. Arrows indicate axis origin and direction.

Fig. 5

Reconstruction of a backward-mode illuminated graphite rod held horizontally above the center of the detector array and at a height of $30\mathrm{mm}$ . (a) Measured velocity potentials (VPs). (b) Estimated VPs output from the reconstruction algorithm. (c) Three orthogonal slices of the reconstructed volume, showing the location of the reconstructed object. The gray scale indicates relative voxel intensity on the range of [0,1]. (d) Photograph of the object taken from the perspective of the laser (from below). The photograph was registered to the same coordinate system. All dimensions reported in mm. Axes tick marks are $5\mathrm{mm}$ apart. Arrows indicate axis origin and direction.

3.2.

3-D Images of Phantoms

The reconstruction of a single synthetic point source is shown in Figs. 3a, 3b, 3c . Fig. 3a displays three of the VPs calculated from the measured PA signals $\left[{\mathit{VP}}_i^{\exp }\left(t\right)\right]$ . In Fig. 3b, three of the ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ , obtained after 500 iterations ( $4\min$ run time on the PC) are shown. The magnitude, shape, and arrival time of the peaks in ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ are well reproduced. Figure 3c presents three orthogonal sections of the reconstructed volume, showing that the point source was reconstructed at the expected location, but with broadened dimensions. It is interesting to note that the broadening was stronger in the $x$ and $y$ directions, and the edges were better defined in the $z$ direction.

Similarly, reconstruction of a single, externally illuminated point source is presented in Figs. 3d, 3e, 3f. In this experiment, only 10 of the 14 detector elements recorded signals above the noise level. Although the exact cause of the signal dropout could not be determined, it may have been related to asymmetry in the PA wave generated by the source. Nevertheless, the algorithm was able to reconstruct the point source at the proper position even with signals from an incomplete set of detectors. Comparison of the ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ in Fig. 3e with ${\mathit{VP}}_i^{\exp }\left(t\right)$ in Fig. 3d demonstrates that some detector signals that were not present in the ${\mathit{VP}}_i^{\exp }\left(t\right)$ (e.g., solid curve) were reconstructed in ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ .

The same VP signal that “reappeared” in Fig. 3e is also present in Fig. 3b, and at the same arrival time, indicating that the reconstruction could accurately compensate for the missing data. Image artifacts, in the form of scattered dim voxels, appear in most of the images. However, these artifacts seem to be more pronounced in images where ${\mathit{VP}}_i^{\mathrm{est}}\left(t\right)$ did not reproduce the shape and magnitude of ${\mathit{VP}}_i^{\exp }\left(t\right)$ well.

Next, we reconstructed an array of three synthetic sources arranged along the $z$ axis [Figs. 4a, 4b, 4c ] and three synthetic sources along the $y$ axis [Figs. 4d, 4e, 4f]. As can be appreciated from the orthogonal slices, the sources along the $z$ axis were better defined and better separated than the sources along the $y$ axis. This was to be expected from the detector arrangement relative to the imaged volume: a separation of $5\mathrm{mm}$ in the $z$ axis created a time-of-flight difference approximately two-times longer than the same separation along the $y$ axis. Therefore, for the sources distributed along the $z$ axis, the reconstruction algorithm was presented with well-separated peaks, which were better reconstructed and resulted in sharper edges.

In both Fig. 4c and Fig. 4f, the intensity of the reconstructed objects varied with location in the image even though the original PA source used to synthesize them was constant in intensity. The objects seemed to reconstruct brighter when closer to the plane of the detectors [Fig. 4c] and when closer to the axis of the detector ring [Fig. 4f]. A potential explanation for the inhomogeneity in voxel intensity could be related to the detection geometry. For example, objects closer to the detectors are better separated in time, and hence better reconstructed. On the other hand, on-axis objects will be covered by a greater number of detectors, which will also result in better reconstruction.

Last, the reconstruction of the graphite rod is presented in Fig. 5 . Comparing Fig. 5a with Fig. 3d, we find that the PA signals recorded on some of the detectors were substantially stronger with the rod than with the point source, while on other detectors, they were weaker. This was likely due to the different propagation of the PA waves from a point source (isotropic) versus a line source (perpendicular to the line). Detectors on the ring that view the rod perpendicularly to its length will hence detect stronger PA signal than detectors that view the rod from a direction parallel to its length. Interestingly, we observe again that even when signals were absent at some detectors, they were generated by the reconstruction algorithm at their expected location. This is another indication of the robustness of the algorithm. The images in Fig. 5c reproduce the rod shape quite reliably. The cross section along the $z$ direction, as can be seen from the $xz$ and $yz$ slices, is only one voxel wide, which agrees well with the rod diameter of $0.9\mathrm{mm}$ . The in-plane orientation of the rod, as indicated by the tilt in the $xy$ slice, is also in very good agreement with the photograph in Fig. 5d. It is worthwhile to note that the PA images of the rod were obtained with a single laser shot and no averaging.

Although the image reconstruction algorithm had robust qualities, some image artifacts were found, appearing as bright voxels near the boundary of the imaging volume [lower-left corner in the $xy$ slice of Fig. 5c] similar to those observed for the synthetic sources.