1. INTRODUCTION

Optoacoustic tomography in living small animals was first demonstrated more than a decade ago [1,2], while new acquisition geometries are being constantly proposed to optimize imaging performance in targeted applications [3–7]. In general, to efficiently record signals generated by absorbers with multiple sizes and orientations and reproduce accurate and quantitative optoacoustic images, the tomographic data should be ideally captured with the broadest tomographic coverage and highest number of detectors surrounding the imaged object [8]. In reality, physical considerations related to signal-to-noise ratio (SNR) as well as geometrical and digital throughput limitations often lead to sub-optimal solutions based on dimensionality reduction. For instance, signals measured by a spherically focused transducer within its depth of focus can be assumed to approximately represent axial (depth) profiles of the optical absorption. Two- (2D) or three-dimensional (3D) imaging is then enabled by raster scanning the transducer and stacking the acquired signals [9–12]. Reduction of imaging dimensionality is also the basis of the widespread cross-sectional optoacoustic imaging approaches. In this case, optoacoustic signals are commonly collected with cylindrically focused transducers [13–15], which can be further engineered to form arrays with sufficient angular coverage for accurate reconstruction of 2D slices for each excitation laser pulse [16,17]. In this way, real-time optoacoustic imaging has been shown to enable the study of a number of dynamic biological processes in single cross sections of small animals [18–20].

Acquisition of whole-body images can be further facilitated by translating 2D or linear imaging arrays along their elevational direction [21,22], in which case the spatial resolution has been shown to be severely compromised by the focusing characteristics of the detection array as well as out-of-plane signal artifacts [23]. Most importantly, the existing systems for 3D whole-body small animal imaging do not support real-time volumetric imaging capabilities [3,4,24,25], making them inadequate for concurrent visualization of dynamic processes.

In this work, we describe a novel small animal tomographic imaging approach based on rotational scanning of a spherical matrix detection array. While the array is quickly translated and rotated around the imaged animal to attain sufficient tomographic coverage for reconstructing high-quality images at the whole body level, volumetric images are simultaneously rendered from selected regions (e.g., whole organs) in real time and are automatically co-registered with the whole-body data.

2. MATERIALS AND METHODS

A. Experimental Setup

Layout of the experimental system is depicted in Fig. 1. For optimal acquisition of tomographic data from multiple angles, a spherical matrix transducer array (Imasonic SaS, Voray, France) was rotated around the imaged sample and perpendicularly translated along its elevational direction, essentially performing a spiral-like motion. The azimuthal and elevational positions of the array can be independently controlled by means of two motorized stages (IAI Inc., Japan) within a range of up to 15 cm and 360°, respectively. The custom-made holder attached to the rotation stage further enables manually setting the radial location of the array with respect to the rotation axis. The imaging samples were immersed in a water tank to ensure optimum ultrasound transmission. The water temperature was kept constant during each experiment with a computer-controlled heating stick.

 

Fig. 1. Layout of the experimental system. The mouse is mounted vertically in a heated water bath and subsequently scanned by translating and rotating a spherical transducer array along the $z$ and $\phi$ axes, as marked by the dashed arrows. Optoacoustic signals are excited using nanosecond-laser pulses. The time-resolved optoacoustic signals generated by each laser pulse are recorded at 256 individual detector locations of the spherical array, time-sampled by the data acquisition system (DAQ), and saved to PC for further processing and image reconstruction.

Download Full Size | PDF

Optoacoustic signals were generated with a short-pulsed ($\sim 10\mathrm{ns}$) optical parametric oscillator (OPO)-based laser source (SpitLight, Innolas Laser GmbH, Krailling, Germany) tunable in the 690–950 nm range, even though the wavelength tuning capability was not exploited in the current study, which was solely performed at 800 nm. The pulse repetition frequency (PRF) was set to 100 Hz. The excitation light was guided to the imaged object through an opening in the center of the matrix array using a custom-made fiber bundle (CeramOptec GmbH, Bonn, Germany). In this way, a Gaussian illumination profile was created at the tissue surface with a full width at half-maximum (FWHM) of approximately 10 mm. The per-pulse energy at the fiber output was kept below 15 mJ. With an approximate diameter of the illuminated area of 1 cm, the fluence at the tissue surface corresponds to approximately $19\mathrm{mJ}/{\mathrm{cm}}^2$.

The spherical matrix array, whose design is described in detail elsewhere [26], consists of 256 individual piezocomposite detection elements. Each detection element has an approximate area of $9{\mathrm{mm}}^2$ and 100% bandwidth around a central frequency of 4 MHz. The recorded time-resolved optoacoustic signals were digitized at 40 mega samples per second with a custom-made high-speed data acquisition system (Falkenstein Mikrosysteme GmbH, Taufkirchen, Germany) that was triggered with the $Q$-switch output of the laser. For each scanning position, the data recorded by all the 256 channels was saved in the computer for further processing. Data acquisition and motor positioning was computer-controlled using MATLAB (R2013a).

B. System Characterization

To accurately calibrate the relative position and orientation of the array with respect to the rotation axis, a gauge phantom consisting of a single 100 μm polyethylene microsphere (Cospheric, Santa Barbara, California) embedded in agar was scanned before each experiment. Specifically, the radial position, lateral shift, and axial rotation angle of the spherical array were determined by considering 35 angular locations (projections) of the array for each elevational ($z$ axis) position (see Fig. 1). Following the initial calibration, two phantoms were imaged to precisely characterize the spatial resolution of the system across the entire field of view. The first phantom consisted of 100 μm polyethylene microspheres (Cospheric) randomly distributed in a 25 mm agar cylinder. The phantom was scanned along 50 projections covering 360° and five translational positions along the $z$ axis with 2 mm steps. To improve accuracy, the acquired signals were averaged 100 times for each scanning position. To evaluate dependence of the spatial resolution on the radial position of the array, two different scans were performed with the focus of the spherical array positioned at a distance of 5.5 or 8.7 mm from the rotation axis. The second phantom consisted of two black 125 μm diameter surgical sutures (Ethicon, USA) arranged to form a cross along the translational scanning axis. In this case, the array’s focus was placed at a distance of 8.7 mm from the rotation axis and the phantom was scanned over 360° with 50 equidistant projections and translation range of 2 cm (2 mm step size). The acquired signals were again averaged 100 times for each position.

Both phantoms were prepared by solving agar powder (1.3 g/100 ml) in deionized water. Intralipid 20% by volume (1 ml/100 ml) was added to the agar solution to achieve more uniform excitation light distribution within the volume containing the absorbing spheres.

C. Signal Processing and Image Reconstruction

The image reconstruction process was conducted in two steps. First, the acquired time-resolved signals were deconvolved with the impulse response of the transducer elements and bandpass filtered between 300 kHz and 6 MHz using a second-order Butterworth filter. The processed data from the 256 detection elements of the array were subsequently used to reconstruct 3D images individually for each scanning position. The images were reconstructed on an anisotropic Cartesian grid of points centered in the focal point of the transducer array using a GPU implementation of a three-dimensional backprojection algorithm [27]. All the individual 3D images were combined by added them up to the final image volume in a second reconstruction step. This was done by using the 3D position data of the spherical array with respect to the rotational axis. For all phantom experiments, the reconstruction grid consisted of 33 μm sized voxels, while a coarser 100 μm grid was used for reconstructing in vivo data from mice. The two-step approach ensures that the real-time data from a smaller field-of-view can be seamlessly projected onto the whole-body anatomical image by using the same reconstruction framework. The reconstructed data was displayed using maximum intensity projections (MIPs) while thresholding was applied to minimize background noise in the final image. All data management and analysis was carried out on an Intel Core i7 machine (8-core, 2.8 GHz) with 64 GB memory using MATLAB (R2013a). The proposed frame sorting approach via k-means was implemented using the available MATLAB (R2013a) routine, whereas the backprojection code was implemented in OpenCL on an AMD Raedon HD 7900 graphics card. A typical processing and reconstruction time for a whole-body volumetric mouse scan was in the 10 min range for the 100 μm voxel size.

Note that while the particular method employed for rendering whole-body images consists of adding up individual small volumes in the image domain, the entire mouse volume could be alternatively reconstructed by using all the projection data concurrently. However, given the particular illumination and scanning geometry, both reconstruction approaches will result in practically identical images. Indeed, only a small volume is illuminated for each projection, and the single projection data contains no information from unilluminated areas.

D. In Vivo Validation

The in vivo imaging capacity was tested with two 6- and 7-month-old female athymic nude-Foxn1nu mice (Harlan Laboratories LTD, Switzerland). The mice were imaged in full compliance with the institutional guidelines of the Helmholtz Center Munich and with approval from the Government District of Upper Bavaria. Their eyes were covered with vet ointment (Bepanthen, Bayer AG, Germany) to ensure protection from the laser light as well as to prevent eye dehydration during the scanning. The focal point of the spherical array geometry was set at a distance of 8.4 mm from the center of rotation, corresponding to $\sim 5\mathrm{mm}$ average depth from the surface of the mouse.

Mouse imaging was performed entirely non-invasively with a self-developed stainless steel mouse holder consisting of two circular plates separated by a thin rod (Fig. 1). Mice were anaesthetized with isoflurane (2%–3% by volume with 0.9 l/min gas flow) through a custom-made breathing mask, while the head was positioned above the water surface at all times. The water temperature in the tank was kept at 34°C to stabilize body temperature during data acquisition. Depending on the number of signal averages at each scanning position, acquisition of a full-view (360°) whole-body mouse dataset over a length of 10 cm takes from 5 to 20 min.

3. RESULTS

A. Spatial Resolution Characterization

Figure 2(a) shows a representative MIP along the $z$ axis of the reconstructed phantom volume. One may recognize that the microspheres are accurately reconstructed across the entire volume. Yet, lower intensity values are attained in the central region of the phantom, which is attributed to heterogeneous light distribution through the turbid phantom. Multiple spheres at different radial distances from the rotational axis were then selected and the spatial resolution performance was determined in cylindrical coordinates as a function of the radial $e_r$, azimuthal $e_{\phi }$, and elevational $e_z$ positions. For this, a FWHM was measured for every reconstructed sphere individually [Figs. 2(b)–2(d)]. The values for $e_r$, $e_{\phi }$, and $e_z$ were then determined as the root square difference of the measured FWHM and the diameter of the spheres [5]. The measured resolution as a function of the radial position of the spheres is displayed in Figs. 2(e) and 2(f). It was found that the radial resolution remains constant throughout the imaged volume in the range between 150 and 250 μm. Similar behavior is observed for the azimuthal resolution, although stronger dependence exists on the position of the center of the spherical array geometry relative to the rotation axis. Elevational resolution exhibits the most significant variations across the imaged volume while also strongly depending on the position of the array with respect to the axis of rotation. Thus, the latter must be carefully selected depending on the total size of the imaged object. Anisotropic spatial resolution across such a large field of view is generally expected due to the relatively high directivity of the individual elements of the ultrasound array along the scanning trajectory.

 

Fig. 2. Spatial resolution within the imaged field of view was determined using an agar phantom containing randomly distributed 100 μm diameter absorbing spheres. (a) MIP of the reconstructed phantom volume along the $z$ axis. The focus of the spherical array geometry was set at a distance of 8.7 mm from the rotational axis. An exemplary region of interest is defined by the orange box. (b)–(d) A zoomed-in image of a reconstructed sphere, as indicated by an orange box in panel (a); three MIPs are shown from orthogonal directions. (e) and (f) show dependence of the radial, azimuthal, and elevational resolution on the radial position from the rotational axis. The focus of the spherical array was set at 5.5 mm and 8.7 mm from the rotational axis in (e) and (f), respectively.

Download Full Size | PDF

Results of the second phantom experiment are shown in Fig. 3. The reconstructed volume was cropped to fit the dimensions of the cross. The sutures appear to have a clear round shape, while no imaging artifacts associated to volume stitching or incomplete tomographic coverage are noticed. Three cross sections through the reconstructed volume were subsequently analyzed [Figs. 3(b)–3(d)]. For the cross sections in 3(b) and 3(d), the FWHM suture diameter was in the range of 266 to 333 μm, in general agreement with the range of resolution values measured in the first phantom experiment. The effective resolution could be alternatively estimated by finding the closest slice to the intersection in which the sutures could still be separated. With this approach, the distance between the two intensity maxima corresponding to the two sutures is associated to what is commonly known as the Sparrow resolution criterion, and it was determined to be 230 μm. A rotational 3D visualization of the reconstructed phantom volume can be further found in Visualization 1.

 

Fig. 3. Volumetric reconstruction of the phantom containing crossed sutures. (a) MIP of the reconstructed volume along the $y$ axis. The focal point of the matrix array was set to 8.7 mm from the rotational axis. (b)–(d) Single cross sections through the reconstructed volume at $z$ axis positions indicated by the pink lines in (a). The pink line in the inset in (c) indicates distance between the two maxima used to determine the resolution according to the Sparrow criterion.

Download Full Size | PDF

B. In Vivo Imaging

During collection of whole-body in vivo data, heart beat and respiration may cause significant motion, which in turn may lead to image blurring and loss of contrast [28,29]. This issue is commonly addressed by gated data acquisition [30]. In the current study, no gating was employed. Instead, 50 frames were collected at each scanning position and a clustering approach was applied. Considering a normal respiration frequency of 1–2 Hz in mice [28], the 50 frames typically correspond to less than two full breathing events at the 100 Hz laser pulse repetition rate, thus about 10–13 frames are effectively affected per each respiratory cycle. It was found that these frames exhibit low frame-to-frame correlation [Fig. 4(a)] as compared to non-respiratory frames. Consequently, they were removed using a k-means sorting algorithm applied to the thresholded correlation matrix of the frame series. An exemplary case of a frame series acquired from the mouse heart and subsequently sorted into two bins can be found in Visualization 2. Figure 4(b) shows a MIP image of the reconstructed whole-body mouse volume—note that, in this case, the 50 frames acquired for each scanning position were averaged before reconstruction. Figure 4(c) shows instead the equivalent image after applying the frame rejection criterion, clearly demonstrating a significant contrast enhancement and better visibility of smaller blood vessels after rejecting the respiratory frames.

 

Fig. 4. (a) Exemplary plot of the correlation matrix of a series of 50 volumetric optoacoustic image frames acquired at a PRF of 100 Hz. Frames with respiratory motion can be identified by low values of the correlation coefficient, whereas frames with no respiratory motion show significantly higher correlation. (b) The 3D reconstruction acquired from a living mouse (MIP along the $y$ axis) where all 50 frames at each single position were averaged before the stitching. (c) The same dataset reconstructed by using only the non-respiratory frames. Scale bar corresponds to 1 cm.

Download Full Size | PDF

Figures 5(a)–5(c) display the same reconstructed volume from different view angles (left backside, backside, and right backside, respectively) with major anatomical structures labeled. Main organs like the kidney or the spleen as well as the spine and the surrounding vascular system can be clearly identified. A different view from the front and left sides is shown in Figs. 5(d) and 5(e), respectively. Here the heart, liver, and thoracic vessels are mainly visible. Fine vascular structures around the organs can also be visualized. A video showing rotational 3D views of the reconstructed data is further available in Visualization 3. Note that the MIP views are naturally emphasizing the superficial signals due to their higher intensity, thus concealing signals from deeper areas. A fly-through movie of the cross-sectional slices is further available (Visualization 4), where deeper structures are visible also in areas close to the center of the mouse.

 

Fig. 5. Whole-body 3D image acquired from a living intact mouse. MIP views are shown from the (a) left-back, (b) back, (c) right-back, (d) front, and (e) right-front sides. (Legend): 1, left atrium; 2, cardiac ventricles; 3, liver; 4, spleen; 5, thoracic vessels; 6, kidney; 7, spine; 8, brown adipose tissue. A 3D rotational video as well as fly-through video of the cross-sectional images are further available in Visualization 3 and Visualization 4. Scale bar corresponds to 1 cm.

Download Full Size | PDF

Finally, the real-time volumetric imaging capacity of the imaging system was investigated by positioning the effective field of view of the spherical matrix array around the heart region. Two-hundred consecutive volumetric frames of the beating mouse heart were acquired at this position with a pulse repetition frequency of 100 Hz without applying signal averaging. Figures 6(a)–6(d) show MIPs of the reconstructed heart images for consecutive time points superimposed onto the previously acquired whole-body anatomical image. The high volumetric imaging speed of the system enables motion-artifact-free differentiation of the different states of the heart cycle on a beat-by-beat basis. The actual heart motion can be best perceived in Visualization 5, which shows the entire series of 200 frames slowed down to 25 frames per second.

 

Fig. 6. Real-time overlay of the beating mouse heart onto a whole body anatomical image of the same mouse. Different phases of the cardiac cycle are shown: $t=0\mathrm{ms}$, ventricular systole; $t=30\mathrm{ms}$, ventricular filling; $t=90\mathrm{ms}$, ventricular diastole; $t=110\mathrm{ms}$, ventricular ejection. The full image series is available in Visualization 5. Scale bar corresponds to 1 cm.

Download Full Size | PDF

4. DISCUSSION AND CONCLUSION

A new approach to volumetric optoacoustic small animal imaging at multiple spatial and temporal scales was presented in this work. The suggested imaging system is based on the combination of rotation and translation of a spherical matrix array transducer array around the imaged object, whereas, for each position of the array, a 3D image covering an approximate volume of $10\mathrm{mm}\times 10\mathrm{mm}\times 10\mathrm{mm}$ can be rendered in real time. All the individual volumetric reconstructions are subsequently combined in an extra reconstruction step to yield a high-quality whole-body image, which in turn can serve as an excellent anatomical reference readily co-registered with the real-time sequence of images acquired from a limited region (e.g., an organ).

The actual imaging performance of the suggested scanning approach for whole-body imaging is determined by the effective angular coverage provided by combining all the measurement locations. In this sense, spherical acquisition geometry represents arguably the most convenient approach for collecting tomographic data with finite-size detection elements having strongly directional sensitivity [3,25]. Our scanning approach is less sensitive to directionality of the individual detection elements, which is evinced by the accurately reconstructed shapes of the blood vessels and other anatomical structures with different sizes and orientations. The results further suggest that the radial position of the detector with respect to the axis of rotation can be adapted to optimize the resolution and imaging performance in a given region of interest.

Simultaneous acquisition of signals with a spherical array of transducers has been previously shown to offer unique capabilities for real-time optoacoustic imaging of dynamic biological events in limited volumes [29,31,32]. Here it was further demonstrated that continuous 3D imaging of the mouse heart is possible with frame rates of 100 Hz. We also exploited the fast volumetric image acquisition rate to correct for motion-associated artifacts in in vivo experiments, thus leading to a consistent enhancement in the spatial resolution performance and contrast of the rendered images. The proposed implementation allows for sacrificing the high temporal resolution performance to gain larger fields of view, thus facilitating better interpretation of the functional data captured in real time.

To this end, a simple reconstruction approach based on stitching of the individual volumes rendered with a backprojection formula has resulted in whole-body mouse images of unprecedented quality. Clearly, more advanced reconstruction schemes, such as model-based inversion algorithms [33], can be attempted instead for further image quality enhancement, e.g., by efficiently accounting for speed of sound variations in the object [34], scattering or dispersion in acoustically mismatch regions, or the effects of limited bandwidth and finite size of the ultrasound detection array [35–37]. Additional hardware developments may lead to better performance, e.g., via continuous (unaveraged) signal collection to reduce the total acquisition time required for the whole-body scans [5]. Alternative scanning protocols could be used to optimize the trade-off between effective field of view and scanning time [4], where the volumes captured in real time can be further optimized via different spherical array designs. Finally, the spatial resolution in optoacoustic tomography scales with the frequency bandwidth of the ultrasound detection elements, so can also be adjusted to optimize the trade-offs among acquisition speed, spatial resolution, and field of view.

In summary, we report on a novel concept for small-animal optoacoustic tomographic imaging that can deliver whole-body images of unparalleled quality while retaining real-time volumetric imaging capability within selected regions at the whole organ scale. The newly introduced system is ideally suited for concurrent data acquisition across multiple spatial and temporal scales, thus offering powerful capacities for simultaneous anatomical, functional, and molecular imaging with optoacoustics.