Adaptive dual-speed ultrasound and photoacoustic computed tomography.

Full-ring dual-modal ultrasound and photoacoustic computed tomography has unique advantages of nearly isotropic spatial resolution, complementary contrast, deep penetration, and full-view detection. However, the imaging quality may be deteriorated by the inaccurate sound speed estimation. Automatic determining and compensation for sound speed has been a long-standing problem in image reconstruction. Here, we present new adaptive dual-speed ultrasound and photoacoustic computed tomography (ADS-USPACT) to address this challenge. The system features full-view coverage (360°), high-speed dual-modal imaging (10-Hz), automated dual sound speed correction, and synergistic high imaging quality. To correct the sound speed, we develop a two-compartment method that can automatically segment the sample boundary and search for the optimal sound speed based on the rich ultrasonic pulse-echo signals. The method does not require the operator's intervention. We validate this technique in numerical simulation, phantom study, and in vivo experiments. The ADS-USPACT represents significant progress in dual-modal imaging.

Introduction

Photoacoustic computed tomography (PACT) has generated increasing interest in medical imaging because of its advantages, such as molecular and functional contrasts and high resolution [1], [2], [3], [4], [5], [6]. Ultrasound (US) imaging provides rich anatomical information and has the merits of low cost, great portability, high biosafety, and deep penetration [7], [8]. Using one set of the ultrasonic transducer and data acquisition system, PACT and ultrasonography can be seamlessly combined for co-registered images with complementary contrasts [9], [10]. Dual-modal US/PA imaging has demonstrated attractive prospects in early cancer diagnosis [11], [12], [13], [14], imaging-guided biopsy, and surgical navigation [15], [16], [17]. Especially, full-ring dual-modal US/PA imaging offers full view angle and fewer streak artifacts and has been widely used in various preclinical and clinical imaging [18], [19].

However, due to inaccurate speed of sound (SoS) selection in reconstruction, full-ring US/PA tomography may suffer from feature splitting or fringing artifacts, severely lowering the imaging quality [3]. The SoS may vary a lot among different tissues. Various approaches that estimate and compensate for either the SoS map or dual SoSs have been developed to improve the imaging quality. An accurate SoS map can be iteratively calculated from PA features [20], but this method might be influenced by non-uniform illumination, heterogeneous light attenuation, and an insufficient signal-to-noise ratio [21], [22]. Transmission-mode US computed tomography (TUCT) can also estimate the SoS map [19], [23], but the iterative computation is time-consuming and may suffer from instability. A simple but effective alternative is a dual-speed method, where the field of ultrasonic propagation is divided into two compartments. Each compartment is assigned a specific SoS value. This method can significantly reduce the computational burden and achieve better imaging quality than using one SoS in the entire imaging field. However, the dual SoSs method has two technical challenges. One is how to efficiently segment the two compartments, i.e., the sample and the surrounding water. The other one is how to determine the optimal SoS values [2]. Many reported dual SoSs methods need manually segment the two compartments or visual tuning the SoS values and, thus, are time-consuming and vulnerable to subjective interference.

Here, we present adaptive dual-speed US and PA computed tomography (ADS-USPACT) that can automatically segment the two compartments and determine the optimized SoS values. The system uses a 256-element ring transducer and provides full-view US and PA imaging. To correct the SoS, we implement an adaptive dual-speed image reconstruction algorithm. The sample is automatically segmented from the surrounding water, the sound speed in the sample is calculated from the US imaging, and the sound speed in the water is determined from a temperature sensor. This method is fully automated and does not need any operator intervention. Furthermore, because the US image has rich features and is not affected by non-uniform light illumination, the US-assisted segmentation and SoS optimization are more robust than the PA feature-driven methods. The ADS-USPACT can acquire co-registered US and PA images at a 10-Hz frame rate, which mitigates motion artifacts. We conduct numerical simulation, phantom experiments, and small animal imaging to demonstrate the imaging abilities.

Materials and methods

ADS-USPACT imaging system

As illustrated in Fig. 1a, the ADS-USPACT system accommodates a customized 256-element ring array ultrasonic transducer (Doppler Electronic Technologies Corp, China). The central frequency is 6.25 MHz, the one-way (for the PA, 6 dB) bandwidth is 76.8%, and the two-way (for the US, 6 dB) bandwidth is 58.4%. Each element is cylindrically focused in the elevational direction (radius: 40 mm; numerical aperture: 0.2; pitch: 0.98 mm; elevation size: 14 mm; or see Supplementary Fig. S1). A 256-channel data acquisition (DAQ) system (Vantage-256, Verasonics Inc.) is one-to-one mapped with the transducer. Each channel can be controlled flexibly in transmission or receiving mode and has an independently programmable amplifier up to 54 dB. Both the US and PA signals are sampled at 25 MHz with a 14-bits resolution. For PA imaging, a Q-switched Nd: YAG laser (Spectra-Physics, Santa Clara, CA, USA) excite the sample with 5–8-ns pulses at a maximum frequency of 20 Hz. The laser beam is coupled into a 14 fiber bundle, and the distal end of each branch is four 32 × 1-mm2 rectangles. The distal ends are fixed around the US transducer array and are oriented at 60° relative to the imaging plane (Fig. 1a). The arrangement achieves an approximately uniform illumination for the sample. We use a 1064-nm wavelength because of high pulse energy, and low scattering [24].

Fig. 1

Adaptive dual-speed ultrasound and photoacoustic computed tomography (ADS-USPACT) platform. (a) Schematic of the imaging system and time sequence for signal acquisition. (b) Close up of animal holder in (a). The holder is used for small animal trunk imaging. (c) Diagram of the US data transmission and acquisition for US imaging. The red dot represents the sequentially activated transmission element, and the green dots on the same side as the transmitting element represent 128 elements that receive the pulse-echo signals. (d) Diagram of PA data acquisition. All 256-element receive PA signals after one laser pulse excitation. The fan region in (c) and rectangle region in (d) enclosed with white solid lines define the reconstruction regions. DAQ, data acquisition system; Rcv, receive; Tx, transmit; TM, thermocouple.

We use a synthetic transmit aperture (STA) approach to collect US signals. All elements are sequentially excited with a 1-cycle sinusoidal pulse and a 250-μs interval. Following each transmission event, only 128 elements on the same side as the emitting element receive the back-scattered US signals (Fig. 1c). The rest elements are electronically turned off. For the PA imaging, all 256 elements receive signals simultaneously within 50 μs after each laser pulse (Fig. 1d). The time interval between the end of US imaging and the next PA imaging is determined by the pulse repetition rate of the laser. The DAQ will trigger the laser Q-switch for the starting of PA imaging at 10-Hz (Supplementary Fig. S2). The time sequence in Fig. 1a shows the interleaved US and PA signals acquisition. The ADS-USPACT system achieves a 10-Hz US/PA frame rate.

Adaptive dual-speed image reconstruction

Ultrasonic-assisted automatic sample segmentation

We develop an ultrasonic-assisted sample segmentation method that can automatically separate the sample region from the coupling medium. The segmentation algorithm involves a 4-step procedure as shown in Fig. 2a. In step 1, at each US transmitting/receiving event, we extract an A-line signal received by the transmitting element. By the end of one US image acquisition, we group the 256 Aline signals into a transmission (Tx) map. In step 2, we detect the first US arrival positions using an improved Akaike information criterion (AIC) picker algorithm [25]. The first arrival positions represent the boundary between the sample and the coupling medium. In the algorithm, a temporal window with a length N = () is determined on each A-line signal. Here is the ring array radius, is the sound speed of the coupling medium (water), and is the sampling frequency. The AIC values are calculated aswhere () divides the temporal window into and . is the variance function, is the mean value in the range of to . Then we calculate an Akaike weight at the position aswhere is .

Fig. 2

Schematic diagram of adaptive dual-speed ultrasound and photoacoustic image reconstruction. (a) Flow chart of dual-speed image reconstruction. (b) Close-up image after the AIC picker, showing the first arrival positions. (c) The sound speed in water changes with temperature. (d) Golden section search (GSS) method for fast SoS searching. The blue solid line represents uniform-step searching results, while the sparse red dots show the GSS results. AIC, Akaike information criterion; amp., amplitude; Exp., experiment; MSFM, multi-stencil fast marching method; recon., reconstruction; SoS, speed of sound; ToF, time of flight.

In step 3, we determine the first arrival positions on the transmission channel as

The distance () from the transmission element to the nearest point on the sample surface is . Therefore, we can collect 256 pairs of coordinates on the tissue boundary by coordinate transformation aswhere is the azimuth angle of the transmission channel . In step 4, the 256 pairs of coordinates form a closed loop to segment two compartments.

Adaptive dual-speed assignment

The speed of sound (SoS) value in the coupling medium (Water) is decided by the temperature [26]. ADS-USPACT uses a thermocouple to measure the water temperature. Fig. 2c shows the relationship between the SoS value and the temperature of the water. When we use the correct SoS value in image reconstruction, the US signals can be coherently summed up. To determine the correct sound speed in the sample, we search the maximum coherence factor (CF) of the US signals at different sample SoS values. Firstly, we use a multi-stencil fast marching (MSFM) method to calculate a series of time of flight (ToF) maps under different sample SoS values (Step 5) [27]. The MSFM method considers the acoustic refraction in the propagation path and calculates the shortest ToF between the transducer element and the reconstructed pixel by solving the eikonal equation [28], [29],where is the ToF at the reconstructed pixel (). Secondly, we reconstruct US images with different ToF maps and compute the total CF values of the reconstructed images (Step 6). The CF value is determined fromwhere is an assumed sample SoS. and are the receiving and transmission element indexes. is the number of transducer elements. is the delayed channel data calculated from the ToF map. Finally, the optimal SoS is determined fromwhere is the total CF summation (CFS) of all pixels, is an interpolation function.

We develop a golden section search (GSS) approach to accelerate Steps 5 and 6. Compared with the uniform-step searching method, GSS can significantly reduce the number of calculations of ToF maps and CF. We show that when searching SoS from 1470 m/s to 1620 m/s, the GSS method only needs to calculate 16 SoS values to achieve 1-m/s accuracy (Fig. 2d), which is 4.75 times more efficient than the uniform-step method (Supplementary Movie 1).

Movie S1

Video showing operating mechanism of the Rehab-TENG. During flexion motion of right arm, TENG-1 is in operation, while during extension motion of right arm, TENG-2 is in operation.

US and PA image reconstruction

After segmenting the sample and determining the correct SoSs, we retrieve the optimized ToF and dual-speed map (Step 7). The US and PA images are reconstructed using a weight-based delay-and-sum beamformer (Step 8) [7]. The US reconstruction at each transmission is confined to a fan-shaped region (Fig. 1c). The maximal reconstructed region for PA is 30 × 30 mm2 (Fig. 1d), which is determined by the measured spatial resolution of the system (Supplementary Fig. S3). Both the axial and tangential resolutions of US and PA are smaller than 380 µm. The spatial resolution was measured by imaging a 20-μm-diameter tungsten wire at different positions. Before the beamforming, the US and PA data are pre-processed using a 3rd-order Butterworth bandpass filter (0.05–8.5-MHz). After reconstruction, we process the US image using envelope detection and logarithmic compression to enlarge the display range. For a fair comparison, we only process the bipolar PA images without vascularization filtering and contrast-stretching [3], [19]. All image reconstruction steps are implemented in MATLAB (2019b, MathWorks, USA) on a computer (Inter Core i7 @ 2.60 GHz, 16 GB of RAM, NVIDIA GeForce RTX 2060).

Dual-modal ultrasound and photoacoustic simulation

We simulated adaptive dual-speed US/PA imaging using the k-Wave toolbox [30]. The simulated element number and parameters of the transducer are the same as the experimental one. In the simulation, we neglected the aperture effect to ensure that the blurry errors are only caused by inaccurate SoS. The coupling medium has a 1500-m/s SoS and 1000-kg/m3 density. As shown in Fig. 3a, a numerical phantom is in a circular region with a 30-mm diameter and an SoS of 1538 ± 39.4 m/s (mean ± standard deviation). The maximal acoustic impedance is 1.7 MRayl. For US imaging, we sequentially transmitted 1-cycle sinusoidal pulses from the 256 elements. The receiving sequence is the same as the implementation of the hardware (see Section 2.1). In the PA imaging simulation, a group of absorbers is uniformly distributed by setting the equal spacing along the x-axis and y-axis in the circular region. Each absorber has a diameter of 0.1 mm and a unit initial pressure (100 Pascal). The sampling frequency is set at 53.2566 MHz to guarantee numerical stability in the ultrasound simulation. The computational grid consists of 1200 × 1200 pixels and the pixel size is 0.1 mm. We add the acoustic attenuation () in the simulation. All the computation is implemented on the graphics processing unit.

Fig. 3

Comparison between single- and dual-SoS reconstructions for simulation and phantom experiments. (a) Numerical phantom. Comparison of (b) reconstructed US image and (c) reconstructed PA image for numerical phantom. The single SoS is set at 1517 m/s, and the automatically determined dual speeds are 1500 m/s for water and 1556 m/s for the sample. (d) Photograph of a chicken breast phantom. (e) Comparison of reconstructed PA image for chicken breast phantom using a manually selected single SoS (1508 m/s) and adaptive dual SoSs (1488.2 m/s in water and 1546 m/s in the biological tissue). (f) Zoom-in images in (e). (g) Comparison of the cross-sectional line profile of the absorber near the boundary in (c) using the single SoS and dual SoSs reconstruction. (h) Quantitative comparison of the diameter of pencil lead in the center (P1) and near the boundary (P2) using the single SoS and dual SoSs reconstruction. amp., amplitude; GT, ground truth; FEP, fluorinated ethylene propylene.

Phantom imaging

We conducted phantom experiments to validate and characterize the dual-speed US/PA imaging. We placed chicken breast tissue into an irregular container. The container was made of fluorinated ethylene propylene (FEP) with a wall thickness of 0.2 mm. One pencil lead (0.7 mm) was inserted into the chicken breast, while another (0.3 mm) was placed near the edge (Fig. 3d). We set the water temperature to 21.9°. The water, FEP, and chicken breast tissue induced acoustic impedance mismatch. The cavity in the chicken breast further increases acoustic heterogeneity. The laser wavelength of 1064 nm was used for PA imaging. The laser fluence for phantom imaging was approximately 9.6 mJ/cm2.

Small animal imaging

Eight-week-old nude mice (BALB/c mouse, ~26 g) were used for in vivo imaging. Throughout the experiment, the mouse was positioned in a customized animal chamber (Fig. 1b) in an upright position and was maintained under ~2% vaporized isoflurane at 0.8 L/min. The chamber was immersed in deionized water and its position can be motorized along the elevational direction by a stepper motor. The water temperature was maintained at ~30 °C and monitored with a thermocouple. The wavelength for mice imaging was 1064 nm and the laser fluence on the skin was approximately 11 mJ/cm2, well below the ANSI limit of 100 mJ/cm2. All animal procedures have been approved by the animal ethical committee of the City University of Hong Kong.

Results

Simulation results

We compared the images reconstructed with one SoS and dual SoSs for the numerical phantom (Fig. 3a). Here, the single SoS was manually determined via visual assessment of the PA image sharpness. This method is commonly employed in PA image reconstruction [22]. We selected the SoS values of 1517 m/s to keep most absorbers in focus (top row in Fig. 3c). However, some absorbers close to the boundary are still split (Fig. 3g), which is unavoidable when maintaining the center area in focus. The central absorbers are out of focus if trying to maintain absorbers close to the boundary in focus (Supplementary Fig. S4). The fringing effect is mainly induced by the large acoustic impedance mismatch between the biological tissue and the coupling medium. Using the dual-speed reconstruction, we can improve the image quality significantly (bottom row in Fig. 3c). Almost all absorbers are well reconstructed and in focus. The single SoS result underestimates the sample size (Fig. 3b), while the dual-speed approach can greatly reduce this error.

Phantom imaging results

We determined the SoS (1488.2 m/s) in coupling water from the measured temperature. The optimal SoS (1546 m/s) in the chicken breast tissue is automatically determined from the US coherent factor. The imaging results are consistent with the simulation. The absorbers near the center are correctly resolved under both single and dual SoSs reconstructions (Fig. 3e, f, and h). However, using one SoS, the pencil lead near the boundary is split and deviation from true diameter seriously (0.399 mm versus 0.3 mm), and we can observe fringing artifacts. In comparison, the dual-speed method corrects these artifacts (Fig. 3h).

In vivo animal imaging results

We demonstrated ADS-USPACT of mice at the cross-section of the liver in the abdominal cavity. To maximize PA penetration, we used a 1064-nm wavelength where optical attenuation in tissue is relatively low. In the single SoS reconstruction, we manually select two different sound speeds by visual assessment of the sharpness of prominent features. One sound speed (1508.9 m/s) was selected to keep the main boundary vessels in focus, and another sound speed (1518 m/s) can keep the central vascular features in focus. The single SoS reconstruction results for US and PA imaging are shown in Fig. 4a and b. To suppress the background noises in the US and PA images, we segmented the biological tissue boundary using a snake-based active contour algorithm [31]. The method demands user input and an initial guess. Therefore, it is time-consuming and is susceptible to subjective variation. In contrast, the dual-speed algorithm not only provides automatic contour segmentation but also enhances the dual-modal imaging quality. We compared the two different single SoS results with the dual-speed results. The single SoS results cannot achieve global improvement in imaging quality, i.e., either the central vessels are blurred or the boundary has fringing effects. The dual-speed method can keep both the boundary (R3 in Fig. 4c) and the central vessels (R3′ in Fig. 4d) in focus. Internal organs, such as the liver, can be identified. Using the automatic segmentation mask (Fig. 4e), we can suppress the background noises for the US image (Fig. 4f) and PA image (Fig. 4g) without the operator’s intervention. We overlayed the US and PA images in Fig. 4h, which shows co-registered high-quality dual-contrast images of the internal organs. We continuously acquire 12-second dual-modal images of the mice and automatically reconstruct the images (Supplementary Movie 2). This shows that automated dual-speed reconstruction can process large datasets with high efficiency. The overlayed results show the system provides accurate registration independent of breathing and heart beating (Supplementary Movie 3). Both results demonstrate the stable performance of ADS-USPACT.

Fig. 4

Comparison between single SoS- and dual-SoS reconstruction for in vivo experiments. (a) Dual-modal US/PA image reconstruction using one SoS (1508.9 m/s). The SoS is manually selected to make the main boundary vascular features in focus. (b) Dual-modal US/PA image reconstruction using another SoS (1518 m/s). The SoS is selected to make the internal vascular features in focus. (c) and (d) zoomed-in views of the regions highlighted in solid boxes in (a), (b), and (g). (e) Automatic two-compartment segmentation and dual-speed results. (f) US image and (g) PA image reconstructed with two optimized SoSs. (h) Overlayed US and PA image. The red arrows indicate the defocused US or PA features. AA, abdominal aorta; IVC, inferior vena cava; PV, portal vein; SC, spinal cord.

Movie S2

Video showing operating mechanism of the Rehab-TENG. During flexion motion of right arm, TENG-1 is in operation, while during extension motion of right arm, TENG-2 is in operation.

Movie S3

Video showing operating mechanism of the Rehab-TENG. During flexion motion of right arm, TENG-1 is in operation, while during extension motion of right arm, TENG-2 is in operation.

Discussion

ADS-UPACT provides high-quality US/PA imaging with co-registered dual-contrasts in deep penetration. The system features adaptive image processing, including fully automatic sample segmentation and dual SoSs image reconstruction. The US-assisted segmentation method is independent of mice size and physiological status, for example, respiration. After completing the tissue surface segmentation, it is convenient to model acoustic propagation for image reconstruction. The adaptive dual-speed optimization avoids subjective uncertainty and is more efficient. Moreover, the US-based method is not interfered with by light distribution and thus is more robust than PA-based methods [22]. At present, it takes about 14 min to localize the optimized SoS values using the golden section search method (10 times MSFM calculation) and reconstruct the US/PA images when the sound speed search interval is set from 1500 m/s to 1600 m/s (Supplementary Table 1). If the mice are fixed in the elevational position, the water temperature reaches equilibrium, and the physiological status is stable, then the ToF map can be calculated only once to achieve real-time display. Via shortening the search interval, using a fast golden section search, and implementing the algorithm on faster computers, ADS-UPACT can fully achieve real-time reconstruction and display in the future.

The ADS-UPACT system shows excellent imaging quality mainly from three aspects. First, the ring transducer offers a full view angle. Second, the system has a video-rate data acquisition, which mitigates motion artifacts. Third, the dual SoSs reconstruction significantly reduces split or fringing artifacts.

The adaptive reconstruction algorithm offers both high imaging quality and fully automated dual-modal image reconstruction. ADS-UPACT represents a major technical advance and facilitates many preclinical and clinical applications, such as cancer diagnosis, surgical navigation, neuroimaging, and drug development.

Author contributions

Yachao Zhang and Lidai Wang conceived the idea and designed the study. Yachao Zhang developed the imaging system and processing algorithms. Yachao Zhang performed experiments. Yachao Zhang and Lidai Wang analyzed the results. Lidai Wang supervised the study. All authors contributed to the manuscript writing.

Declaration of Competing Interest

Lidai Wang has a financial interest in PATech Limited, which, however, did not support this work. All authors declare no competing interests.