Real-Time Passive Acoustic Mapping Using Sparse Matrix Multiplication

Abstract

Passive acoustic mapping enables the spatiotemporal monitoring of cavitation with circulating microbubbles during focused ultrasound (FUS)-mediated blood-brain barrier opening. However, the computational load for processing large data sets of cavitation maps or more complex algorithms limit the visualization in real-time for treatment monitoring and adjustment. In this study, we implemented a graphical processing unit (GPU)-accelerated sparse matrix-based beamforming and time exposure acoustics in a neuronavigation-guided ultrasound system for real-time spatiotemporal monitoring of cavitation. The system performance was tested in silico through benchmarking, in vitro using nonhuman primate (NHP) and human skull specimens, and demonstrated in vivo in NHPs. We demonstrated the stability of the cavitation map for integration times longer than 62.5 [Formula: see text]. A compromise between real-time displaying and cavitation map quality obtained from beamformed RF data sets with a size of 2000 ×128 ×30 (axial [Formula: see text]) was achieved for an integration time of [Formula: see text], which required a computational time of 0.27 s (frame rate of 3.7 Hz) and could be displayed in real-time between pulses at PRF = 2 Hz. Our benchmarking tests show that the GPU sparse-matrix algorithm processed the RF data set at a computational rate of [Formula: see text]/pixel/sample, which enables adjusting the frame rate and the integration time as needed. The neuronavigation system with real-time implementation of cavitation mapping facilitated the localization of the cavitation activity and helped to identify distortions due to FUS phase aberration. The in vivo test of the method demonstrated the feasibility of GPU-accelerated sparse matrix computing in a close to a clinical condition, where focus distortions exemplify problems during treatment. These experimental conditions show the need for spatiotemporal monitoring of cavitation with real-time capability that enables the operator to correct or halt the sonication in case substantial aberrations are observed.

Fig. 1.

Sparse matrix construction. (a) 3D matrix of distance from the pixel to transducer element in sample units calculated for an imaging array with 128 elements, and a reconstructed image of 2000×128=256000 pixels. (b) 3D matrix of re-indexed distance to follow data output by Verasonics Vantage. (c) Standard basis vectors used for the DAS iterative calculation. (d) Matrix resulting from DAS operation on a given standard basis vector. (e) Sparse matrix values allocation following 2D representation.

Fig. 2.

TEA-PAM real-time algorithm using sparse matrix operation. (a) DAS beamforming algorithm using GPU-accelerated sparse matrix operation. (b) TEA operation. (c) Cavitation maps.

Fig. 3.

Experimental setups for (a) the in vitro skull and phantom and (b) the in vivo BBB opening in NHP. In the in vitro experiment, the FUS transducer was placed on top of the phantom and orthogonal to the PCD array. The skull was placed between the phantom and the PCD array for assessing the skull effects on the cavitation mapping. In the in vivo experiment, the FUS transducer was targeted to the region-of-interest based on the neuronavigation coordinates while the PCD array was placed against the temporal bone window toward the FUS focus.

Fig. 4.

Effect of integration time on computational time and cavitation mapping characteristics. (a) The computational time T_c increased linearly with the integration time T_i, which limited the integration time to a maximum of 1.44 μs for PRF= 2 Hz (0.5 s of pulse repetition period). (b) Ti also affected the mapping quality with lower values providing maps with discrete cavitation spots out of focus and values higher than 62.5 μs reaching a steady state of cavitation map. This representative case was performed at 450 kPa. Computational times T_c refers to the reconstruction of large acoustic maps (2,000 axial pixels × 128 lateral pixels).

Fig. 5.

Cavitation mapping sensitivity through primate skull. Cavitation maps using 62.5-μs integration time at variable pressure levels were acquired (a) without a skull, (b) with NHP skull, and (c) with human skull between the PCD array and the phantom.

Fig. 6.

Cavitation activity recorded during BBB opening in NHP 1. (a) BBB opening (in color) induced by sonication at 450 kPa revealed in contrast-enhanced T1-weighted MR image. (b) Frequency spectra obtained from the beamformed signal at the location of maximum image intensity. (c) Cavitation activity evaluated with a single-element PCD transducer (cavitation dose). (d) Normalized power detected with a PCD array positioned at the temporal window and a single-element PCD transducer co-aligned with the FUS transducer. Power was defined as the sum of squared beamformed data using 1.44-μs integration time for the PCD array and the entire signal duration for the single-element PCD transducer. (e) Spatial distribution of the cavitation activity shown in cavitation maps (−6 dB) for different samples.

Fig. 7.

Cavitation activity recorded during BBB opening in NHP 2. (a) BBB opening (in color) induced by sonication at 450 kPa revealed in contrast-enhanced T1-weighted MR image. (b) Frequency spectra obtained from the beamformed signal at the location of maximum image intensity. (c) Cavitation dose detected using a single-element PCD transducer indicating only stable cavitation throughout the sonication duration. (d) Normalized power detected with a PCD array positioned at the temporal window and a single-element PCD transducer co-aligned with the FUS transducer using 1.44-μs integration time. (e) Reconstructed cavitation maps (−6 dB) for different samples.

Fig. 8.

Benchmarking for sparse matrix operation. (a) Computational time of CPU and GPU implementations for the sparse matrix multiplication and standard DAS using the NHP in vivo data set for integrated sample N_T varying from 10 to 1500 (Z: 2000 samples, N: 128 elements, N_x: 128 lateral pixels, N_z: 2000 axial pixels). (b) Sample processing rate at the same conditions. (c) Computational time with N_T = 30 for different sizes of FOV, where the lateral map size N_x remained constant equal to 128 pixels and the axial map size varied: N_z = 50, 100, 200, 500, 1000, and 2000 pixels. The benchmarking was performed offline in a Dell Precision T7910 workstation (dual processor Intel Xeon CPU E5–2650 v4 @ 2.20GHz, 128 GB of RAM) equipped with a GPU (NVIDIA Quadro P6000, 24 GB memory, 3840 cores, driver: 392.56) running MS Windows 10 Pro 64-bits and Matlab 2017b.