The physical properties of transducers used in optoacoustic tomography significantly influence image quality [1]. Modeling and characterization of the total impulse response (TIR) of a transducer element provide a convenient way to capture the combined effects of the spatial impulse response (SIR) due to the aperture size and the electrical impulse response (EIR) arising from the detection bandwidth of the transducer elements. Nevertheless, TIR characterization in a high-resolution field of view (FOV) remains an arduous task. In the case of a handheld optoacoustic probe with limited view, the complexity of this challenge is further compounded by the acoustic speed mismatch between the coupling medium and the sample. Modeling and characterization of a single transducer element can be sufficient in the case of tomographic imaging where a single element is rotated around or translated along the sample to collect signals. However, if transducer arrays are used for simultaneous data recording from multiple channels, as is often done in pre-clinical and clinical optoacoustic imaging systems, the assumption of identical impulse responses is generally not valid. In practice, the EIRs of elements vary along an array, and the importance of characterization of such variability for quality control [2] and identification of defective or weak elements has previously been highlighted in [3]. Failing to account for such individual EIRs from each transducer element in reconstruction algorithms may degrade image quality and resolution.

Previous TIR correction attempts for full view small animal optoacoustic imaging systems [1] and limited view clinical handheld systems [4] used one single EIR for all transducer elements in the array. The EIR derived by a simple deconvolution procedure was not robust to measurement noise. More importantly, the EIRs obtained for multiple array elements were averaged to obtain this single smooth EIR [4]. It is, however, important to develop a robust EIR derivation method for TIR characterization of each transducer element in an array-based handheld optoacoustic tomographic system. Towards offering better image quality and isotropic resolution, we set out to develop a method for individual EIR derivation. We present a robust individual TIR characterization method for transducer arrays using only sparse measurements throughout the FOV. We obtain individual EIRs by formulating the EIR derivation problem as a system of linear equations that can be solved using a least square solver. The EIRs of all the transducer elements are then combined with the simulated SIR to generate the individual synthetic TIR (isTIR) model.

In this Letter, we first illustrate the robust method to characterize the individual transducer response. We then use the isTIR in the model-based reconstruction framework to demonstrate (a) improvements in image quality and isotropic resolution using a physical microsphere phantom, and (b) improvements in resolution of reconstructed images of clinical optoacoustic scans.

TIR model in handheld optoacoustic tomography. For currently emerging limited view handheld transducers in clinical applications, the forward TIR model [4] incorporating transducer properties and refraction of acoustic waves at the interface between sample and coupling medium can be written as

Here, $h_{{r_e}}^{{\rm{aEIR}}}$ is the experimentally derived approximate EIR (aEIR) of the transducer element located at ${r_e}$ and $m_{{r_e},r^\prime}^R$ consisting of the numerically modeled terms $h_{{c_c},{c_t},{r_e},r^\prime}^{{\rm{SIR}}}$, the SIR of the transducer element, and $p_{{c_c}}^N$, the optoacoustic response of a radial absorber. The intensity of the reconstructed initial pressure at location $r^\prime $ is denoted by ${f_{{r^\prime}}}$. The time of flight from the source to the detector is denoted by ${n_R}$. It was assumed that the optoacoustic response remains intact after refraction. Velocity dispersion analysis is provided in Fig. S1 of Supplement 1.

isTIR: measurements on a sparse grid in the FOV. To obtain the left-hand side in (1) for an overdetermined system of linear equations, we performed multiple measurements of a microsphere, placing it at $P$ different locations on a sparse grid in the FOV, as shown in Figs. 1(a) and 1(b).

Modeling SIR and pixel response. We numerically modeled the term $m_{{r_e},r^\prime}^R$ in (2) by discretizing the active surface of each transducer element into ${{50}}\;{{\unicode{x00B5}{\rm m}}} \times {{50}}\;{{\unicode{x00B5}{\rm m}}}$ sub-apertures and calculating the SIR using Field-II [5]. To account for refraction, we assigned the source location with the coordinates of the virtual point [4] and set the acoustic speed to 1397 m/s, corresponding to the speed of sound in a coupling medium (heavy water in our case). The sampling rate in Field-II was set to 40 MHz, the same as the sampling rate of our data acquisition. Then the geometry of the array (60 mm radius covering an angle of 145°) was modeled in MATLAB, and the coordinates of each transducer element were passed to the Field-II program to obtain the corresponding SIR.

Derivation of individual aEIR. We used the same target microsphere and could hence drop the term ${f_{{r^\prime}}}$ in (1) and construct a system of equations for $q$th transducer element using the $P$ measurements [as shown in Figs. 1(a) and 1(b)] at grid locations $p = 1,2,\ldots P$ as

As convolution with the function $m_{q,p}^R$ is a linear operation, we can write the system of Eq. (3) as a linear system using the corresponding Toeplitz matrices that we denote by $Tm_{q,p}^R$:

 

Fig. 1. Method of derivation of the individual transducer response. (a) Schematic of the transducer array with grid locations where signals from a microsphere were measured. (b) Photograph of the setup with microsphere placed in the FOV of the optoacoustic handheld scanner. (c) Individual EIRs derived for different elements of the transducer array. (d) Sensitivity profile across different transducer elements.

Download Full Size | PDF

Individual synthetic TIR. Equipped with the individual aEIR for all $Q$ elements, we combined it with the modeled SIR to generate the full isTIR forward model as

If the location of the pixels at $r^\prime $ in the FOV is indexed by $(i,j)$, where $i = 1,2,\ldots ,M;j = 1,2,\ldots,N,$ for ${\rm{M}}\times {\rm{N}}$ pixels in the FOV, then (7) can be denoted as an isTIR forward model matrix $M = [{m_{q,(i,j)}}]$. Ultimately, arranging the pixel responses (6) can be expressed as a linear system $s = Mf$, which can be solved to reconstruct the image $f$.

Image reconstruction using isTIR: characterization and experimental setup. The generated optoacoustic waves are detected using a cylindrically focused concave array (Imasonic, France) of 256 piezoelectric transducer elements with an angular coverage of $145^\circ$. The elevation radius of the curvature of each element is 65 mm, while the array radius is 60 mm. The center frequency of the transducer is approximately 4 MHz with a ${-}{{6}}\;{\rm{dB}}$ bandwidth of 50% as characterized by the manufacturer in transmit/receive mode. Details of the imaging system can be found in [4]. For characterization of the handheld scanner, a two-axis translation stage was built using motorized linear translation stages from Thorlabs, USA. The stages were programmed to perform a raster scan along a predefined grid in the FOV of the scanner. We performed measurements using a single microsphere phantom immersed in water as shown in Fig. 1(b). These measurements were used to demonstrate the isotropic resolution improvement. We also recorded clinical optoacoustic data by scanning arm cross-sections of healthy volunteers to visualize vasculature under the skin. Ultrasound gel was used as a coupling medium. We received informed consent prior to performing the scans.

Framework for image reconstruction. We reconstructed the initial pressure ${f_{{\rm{sol}}}}$ of all scans using a model-based reconstruction algorithm. We added Tikhonov regularization to counteract the two main causes for the ill posedness of the problem: simple Tikhonov regularization to reduce limited view noise and a Laplacian-based regularization term to suppress sub-resolution noise:

The regularization parameter for Tikhonov regularization, ${\lambda _I}$, was chosen using the L-curve. After fixing ${\lambda _I}$, an appropriate parameter for Laplacian regularization, ${\lambda _L}$, was chosen by visual inspection of the images. The optimization in (8) converges after approximately 50 iterations (see Fig. S3 in Supplement 1). Images were reconstructed in shift-invariant function spaces [6,7], which are defined by a Gaussian-shaped pixel model and a discretization grid with either 25 µm or 50 µm resolution, which allowed us to visualize reconstructed images in arbitrary resolution without pixilation artifacts. The zoomed images were generated using evaluation of the reconstructed images at a discretization of 12.5 µm. Because we aimed to correct for individual transducer responses to enhance the isotropic resolution, it was therefore necessary to address the effects of limited view, e.g., negative values. We therefore used non-negative constrained inversion using the projected conjugate gradient method [8] to reconstruct images to obtain physically meaningful optoacoustic contrast.

Characterization of isotropy and resolution. Figure 2(a) shows the reconstructed image of the microsphere using the sTIR model with an average transducer response, while Fig. 2(b) shows the same microsphere reconstructed using the isTIR. We observed that the inclusion of individual transducer element responses eliminates the artifacts around the microsphere and tends to improve the isotropic shape of the sphere. The radius of the red dashed circle in the frequency domain plots, shown in Figs. 2(c) and 2(d), corresponds to the resolution, which is approximately 113 µm in the lateral direction. The lateral resolutions measured using the full width half maximum (FWHM) of the reconstructed microsphere [shown in Figs. 2(a) and 2(b)] using sTIR and isTIR are 125 µm and 113 µm, respectively. Comparing Figs. 2(c) and 2(d), we can observe that the uniformity of the disc outlined by a red dashed line is enhanced by isTIR correction, demonstrating an improvement in the isotropy and concentration in reconstruction of a microsphere. The black dashed lines shown in Figs. 2(c) and 2(d) mark the limited view sectors that arise from limited angular coverage of the handheld scanner. We observed that the isTIR correction attempts to smoothen the effect of limited view by filling in the blind sectors such that the transition is less abrupt, which implies reduced streak artifacts [9]. Image improvement using isTIR was analyzed by placing the microsphere at ${{6}} \times {{7}}$ grid locations of full FOV (see Fig. S4 in Supplement 1). Figures 2(e) and 2(f) display the improved images of the microsphere using isTIR compared to sTIR, at depths of 5 mm, 11 mm, and 17 mm from the surface of the handheld scanner. The resolution improvement, defined as ${{\rm{FWHM}}_{{\rm{sTIR}}}} - {{\rm{FWHM}}_{{\rm{isTIR}}}}$, at various depths is displayed in Fig. 2(g). The mean resolution improvement was 21 µm (marked using a dashed magenta line).

 

Fig. 2. Improvement in isotropy and resolution in phantom images with inclusion of the individual transducer response. (a), (b) Images of a microsphere located approximately at the center of the FOV of the handheld scanner shown in Fig. 1 using sTIR and isTIR, respectively; scale bar, 1 mm. (c), (d) Log of magnitude 2D Fourier transform of the reconstructed images in (a), (b) using sTIR and isTIR, respectively. (e), (f) Images of microspheres located at three depths reconstructed using sTIR and isTIR, respectively. (g) Boxplot showing the resolution improvement using isTIR compared to sTIR at various depths.

Download Full Size | PDF

Improvement in clinical image reconstruction. Figures 3(a) and 3(b) show reconstructed images of the first clinical scan using sTIR and isTIR models, respectively. Similarly, Figs. 3(c) and 3(d) show the reconstructed images of the second scan using sTIR and isTIR models, respectively. Comparing the first and second columns of Fig. 3, one can clearly see the contrast improvement in most of the vascular structures. To highlight the improvement in resolution, we illustrate the high-resolution zoomed images of two small vessels marked with a dashed line and solid line boxes for sTIR and isTIR reconstructions. Comparing the zoomed-in images of the longitudinal and cross-sectional blood vessels, we again observed that isTIR improves the contrast and sharpness of the resolved vascular structures.

 

Fig. 3. Improvement in clinical images with inclusion of individual transducer response. (a), (b) Images of the arm of volunteer 1 using sTIR and isTIR, respectively. (c), (d) Images of the arm of volunteer 2 using sTIR and isTIR, respectively. Excitation wavelength, 800 nm; scale bar, 5 mm.

Download Full Size | PDF

We herein proposed a robust method to determine the individual EIR of transducer elements in an array and combined it with a refraction-based SIR model to obtain the isTIR. We hypothesized that inclusion of individual EIRs of transducer elements could improve the image resolution isotopically and formulated a system of linear equations using multiple measurements of TIRs at different locations of the FOV to derive the individual EIR for each transducer element. This was feasible as the EIR is independent of the location (distance and direction) of the source relative to the transducer. The proposed method not only overcomes the tedious characterization process of the entire handheld scanner, but also offers a robust characterization of individual transducer responses. The adverse effects of considering an average transducer response [1,4] over the whole array has been highlighted in this work during the reconstruction of a single microsphere using the previously reported sTIR model [4]. The robust characterization of the individual aEIR revealed an inhomogeneity in responses among the elements of the transducer array. The isTIR-model-based reconstruction, which incorporates these individual responses, resulted in a higher level of isotropy in the achieved resolution as depicted from the 2D Fourier transforms of the reconstructed images of a single microsphere. In comparison to the previously reported sTIR model (based on the average transducer response), the isTIR model (based on the individual transducer response) demonstrates improvements in image quality across different locations of the FOV using physical microsphere phantoms embedded in agar. Additionally, we also demonstrated improvements in image resolution of small vessels using reconstruction of clinical scans recorded from healthy volunteers.

The primary advantage of individual sTIR correction over the existing average sTIR correction was observed to be isotropic resolution improvements, which increase the overall image quality. However, the limited view artifacts were still present as observed from the 2D Fourier transform of the microsphere image. Further sTIR improvements on spectral unmixing can be explored in the future. Additionally, the effects of light fluence and ultrasound attenuation can also be integrated into the forward model as a future scope.

In summary, the proposed method to characterize the individual transducer response in the context of a handheld optoacoustic system produces images with superior quality and resolution, is robust in nature, and can be easily extended to any handheld tomography system with any coupling medium with known optoacoustic properties. This can ultimately facilitate high-quality image reconstruction and increase the clinical diagnostic value of handheld transducer arrays.