Research on photoacoustic tomography (PAT) got prosperous development for its being promisingly characterized with noninvasive and nonionizing diagnose of breast cancer, arthritis, and relevant disease. PAT combines the metrics of both ultrasound imaging and pure optical imaging technique, providing high ultrasonic resolution and high optical contrast images^[1–5]. Due to the peculiarity that the optical absorption characteristic of blood has a strong relationship with the hemoglobin content, functional imaging as well as structural imaging can also be realized by PAT, making this imaging modality extremely potential in clinical application^[6].

The basic principle of PAT is that a tissue is irradiated with short nanosecond laser pulses, and then the absorbed energy may result in a thermo-elastic expansion and subsequent contraction of irradiated volume that generates time-trace photoacoustic (PA) waves, which can be acquired by scanning small-aperture ultrasound detectors over a surface that encloses the source under study. The recorded PA wave can then be reconstructed to spatially resolve the initial absorber distribution and concentration via PA reconstruction algorithms^[7–9]. However, biomedical tissue is a highly scattering medium for electromagnetic waves in the optical spectral range, and the propagation ultrasound waves are extremely attenuated before received by ultrasound sensors. Furthermore, the dose of laser beam exposed on the biomedical tissue has to be limited under 20mJ/cm2 for safety operation. Thus, in clinical trials, ultrasound transducer can only receive weak PA signals with low signal-to-noise ratio (SNR) which degraded the final reconstructed PA images^[10–14]. This Letter proposes an adaptive multi-sample-based approach to enhance the SNR of PA signals and in addition, detailed information of PAT that used to be buried in the noise and artifacts can be distinguished.

The PA reconstruction is an inverse problem of the source pressure. We assume that a tissue with inhomogeneous microwave absorption but a relatively homogeneous acoustic property and the heat diffusion’s effect on the thermoacoustic wave can be ignored. For cases where the scanning radius in a circular scan configuration is much greater than the PA wavelengths, the optical absorption p0(r) within the sample at a given position r is given as^[1]p0(r)=−r02Cp2πc4β∫θ0dθ01t∂p(r0,t)∂t|t=|r−r0|/c,(1)where c is the speed of sound; β is the isobaric volume expansion coefficient; Cp is the heat capacity; r0 is the detector position with respect to the imaging cente; p(r0,t) is the PA signals detected at each scanning angle θ0. ∂p(r0,t)/∂t can be calculated through inverse Fourier transform^[15]∂p(r0,t)∂t=FT−1{−iωp(r0,ω)WΩ(ω)},(2)where FT−1 donates the inverse Fourier transform; ω is angular frequency; WΩ(ω) is a window function to depress the high-frequency noise, which is characterized by a cutoff angular frequency Ω. In our work, we chose Hanning window, which has smooth edge shape and less calculation complexity WΩ(ω)={0.5+0.5cos(πωΩ),if|ω|<Ω0,otherwise.(3)

We can discern from Eq. (1) that the SNR of the received PA signals is of great importance to the image quality of PAT. Normally, using single-sampled PA signal can only build a low-SNR PA image. Since the PA signal and the random noise are uncorrelated in each measurement procedure, we can apply averaging frames using multi-sample method to enhance the SNR of reconstructed PA images. However, after applying multi-sample method, with the improvement of SNR, in real biomedical tissue, the spatial resolution and detailed information of reconstructed PAT image is hardly increased (as shown in a subsequent paragraph). This abnormal phenomenon can be explained by the schematic in Fig. 1.

Due to the uncertainty accumulated by the whole experiment system, there always exists tiny time shift between different PA frames. While the interval time between two sequenced pulsed lasers is precisely fixed (max{|T1−T2|}<1ns), the ultrasound signals received are not perfectly aligned (max{|T2−T4|}>0.1μs).

By calculating the cross correlation between the two PA frames p1(t) and p2(t), the time shift ts can be determined. For continuous signals, the cross correlation is defined as p1(t)*p2(t)=∫−∞∞p1*(τ)p2(t+τ)dτ,(4)where p1*(t) is the complex conjugation of p1(t). With the cross correlation^[16], ts can be determined according to the time offset between zero and the highest peak. For discrete signals p1(n) and p2(n), their cross correlation is R^p1p2(m)={1L−|m|∑n=0L−m−1p1(n+m)p2(n)m≥0R^p2p1(−m)m<0,(5)where L is the maximum length of signal p1(n) and p2(n). However, using Eq. (5) can only get a maximal precision of shift time ns between p1(n) and p2(n) at 1/fsample, where fsample is the sample frequency of the receiving system. With N times interpolation on p1(n) and p2(n), the precision of shift time ns will increase N times. In our work, shift time ns is calculated by applying cubic spline interpolation^[17] on receiving signals before using Eq. (5). With correct ns, received PA signals are aligned and averaged as p^(r,Lint)=1Ns∑i=1Nspi_shift(r,Lint),(6)where Lint is the interpolated length of PA signal; Ns is the number of sampled frames; pi_shift(r,Lint) are the shifted PA frames. Finally, N times sample is applied on p^(r,Lint) to generate PA signal p^(r,n) for imaging p^(r,n)≡p^(r,Lint)|Lint=nN.(7)

Our experiment system integrates a tunable optical parametric oscillator laser (Vibrant B, Opotek Inc., Carlsbad, CA, USA) pumped by the second harmonic output of a Nd:YAG pulsed laser (Brilliant B, Quantel, Bozeman, MT, USA) as the illumination source. The laser system is turned to 720 nm which gives the maximum output energy of 60 mJ/pulse. The pulsed laser beam with a repetition rate of 10 Hz and pulse width of 5.5 ns illuminates the imaging domain. A Verasonics system (Model V3, Redmond, WA, USA) was programmed to work with a computer. To build real-time PA images, we applied our hardware-based optimization method in experiments^[18]. In addition, we applied our linear delay compensation optimization method (similar to the concept of beamforming) to build full-range focused PA images^[19]. The schematic diagram of the imaging system is shown in Fig. 2.

The first experiment is conducted with a phantom containing four human hairs. These hairs are fixed on a U-shaped plastic holder. A metal shank is also attached on it so as to make it submerge in water. We use the Philips CL15-7 transducer probe (liner array, 128 elements, center frequency at 9 MHz with 0.178 mm element pitch) to detect the PA signals. The probe is vertical to the hairs in order to get the cross section image of the four hairs. Figure 3(a) shows the model of our phantom. Figure 3(b) is the reconstructed PA image using one sample frame which contains strong background noise and artifacts. Figures 3(c) and 3(d) are the reconstructed PA images with simple multi-frame average method and our proposed multi-frame method, respectively. Both results come from the same 10 sampled frames. In this work, we applied 3× interpolation to rectify the random shift between sample frames. It is apparent that the noise and artifacts in Fig. 3(d) is effectively depressed and the four human hair sections can be easily discerned from the background.

To quantify the effect of our method, we compare the peak SNR (PSNR; as shown in Table 1) of the reconstructed PA images. For all four human hairs, when computing their PNSRs, the background is set to a 50pixels×50pixels region far away from the region of interest. In Table 1, PNSR_S is the PSNR of the PA image reconstructed using one sample frame, PSNR_N is the PSNR of the PA image reconstructed using simple multi-sample averaging method, and PSNR_M is the PSNR of the PA image reconstructed using our multi-sample-based approach. From Table 1 we can see that 6 dB quality improvement is achieved after applying our proposed method.

Another experiment is conducted with a volunteer’s index finger, which is fixed to avoid moving artifacts in the process of PAT. The laser beam is parallel to the finger while vertical to the CL15-7 probe. Figure 4(a) is the reconstructed PA image using one sample frame. Figures 4(b) and 4(c) are the reconstructed PA images of the same target using simple multi-frame average method and our proposed method, respectively. We can see that, as a whole, Fig. 4(b) presents better image quality than Fig. 4(a) especially in the regions where we mark in the colored line (Regions A–C). Figure 4(c) is even better than Fig. 4(b) in presenting clear infrastructures. To quantify the quality of reconstructed PA images, PSNR comparison among different methods is listed in Table 2 which indicates that our proposed method has overwhelming advantages over other ones. Furthermore, we additionally include a B-mode ultrasound image [shown as Fig. 4(d)] of the same target for comparison. This in vivo experiment proves that reconstructed PA image with our proposed method strictly match the biomedical infrastructure (tendon, bone, joint, periosteum, and inner structure of tendon) to those in ultrasound image. To enhance the noise and artifacts, these images are displayed using nonlinear gray level.

In conclusion, we adaptively utilize repeatedly sampled PA signals to yield optimized biomedical PA images which not only suppress the noise and artifacts, but also enhance the detailed information in reconstructed images. We hope our improved multi-sample-based approach will benefit future clinical trials and other in vivo research.