Compressed Sensing-Based Super-Resolution Ultrasound Imaging for Faster Acquisition and High Quality Images

Goal: Typical SRUS images are reconstructed by localizing ultrasound microbubbles (MBs) injected in a vessel using normalized 2-dimensional cross-correlation (2DCC) between MBs signals and the point spread function of the system. However, current techniques require isolated MBs in a confined area due to inaccurate localization of densely populated MBs. To overcome this limitation, we developed the ℓ1-homotopy based compressed sensing (L1H-CS) based SRUS imaging technique which localizes densely populated MBs to visualize microvasculature in vivo. Methods: To evaluate the performance of L1H-CS, we compared the performance of 2DCC, interior-point method based compressed sensing (CVX-CS), and L1H-CS algorithms. Localization efficiency was compared using axially and laterally aligned point targets (PTs) with known distances and randomly distributed PTs generated by simulation. We developed post-processing techniques including clutter reduction, noise equalization, motion compensation, and spatiotemporal noise filtering for in vivo imaging. We then validated the capabilities of L1H-CS based SRUS imaging technique with high-density MBs in a mouse tumor model, kidney, and zebrafish dorsal trunk, and brain. Results: Compared to 2DCC and CVX-CS algorithms, L1H-CS achieved faster data acquisition time and considerable improvement in SRUS image quality. Conclusions and Significance: These results demonstrate that the L1H-CS based SRUS imaging technique has the potential to examine microvasculature with reduced acquisition and reconstruction time to acquire enhanced SRUS image quality, which may be necessary to translate it into clinics.


I. Introduction
Contrast enhanced ultrasound imaging has been developed using microbubbles (MBs) to provide improved therapeutic outcomes and anatomical and functional information for the diagnosis of disease and preclinical investigations. MBs have been utilized for blood-brain barrier opening and drug/gene delivery by stable and inertial oscillations to increase the permeability of connective tissues in brain vasculature [1]- [3]. Targeting specific biomarkers using MBs with functional moieties detects tumor angiogenesis with high specificity using ultrasound imaging [4], [5]. Recently, a series of development of MB mediated ultrasound imaging techniques tries to overcome the diffraction limit of conventional ultrasound imaging while preserving penetration depth with improved detection capability for clinical and preclinical applications. Super-resolution ultrasound (SRUS) imaging technique utilizes the blinking of MBs by compounding thousands of frames to construct microvessel images to examine the structural abnormalities and dysfunction of the vessel in various types of organs and leaky vasculature in tumor mass [6]- [12]. Before constructing the SRUS image, the clutter signals of each frame are removed by applying to transmit or post-processing techniques such as pulse inversion [13], differential imaging [14], and singular value decomposition (SVD) filtering [15]. Among them, SVD filtering was demonstrated as a suitable processing step to precisely localize MBs in the range of few micrometers [16].
Normalized 2-dimensional cross-correlation (2DCC) technique is then utilized to localize MBs by comparing detected MBs signals and the point-spread function (PSF) of the ultrasound imaging system [7]. However, 2DCC algorithm does not accurately detect densely populated MBs. Reduced MBs concentration may be used to avoid inaccurate localization due to dense MBs, which increase data acquisition time. For example, a total of 75000 frames at a frame rate of 500 Hz was acquired to reconstruct the SRUS image of rat brain vasculature within 150 sec [17]. The accompanied drawback with longer acquisition time may cause misalignment between frames due to the motion of an imaging target. Recently, an isolation technique using a 3D Fourier transform was introduced to overcome this limitation by dividing the subset of MBs signal according to the flow velocities and directions at the region of interest (ROI) in a vessel [18]. However, it is challenging to determine optimal subsets due to the unpredictable flow rate in deeper microvessel in vivo, which may cause blurring and distortion of original MBs signals and inaccurate localization of MBs.
showed the capability to recover closely located structures labeled with fluorophores [21]. Furthermore, CS algorithm was demonstrated that this could be utilized for the localization of high-density MBs for the reconstruction of SRUS images in simulation and in vivo application while reducing the acquisition time of ultrasound image frames [22], [23]. However, the standard CS algorithm based on the interior point method is computationally expensive, thus, the computational time was measured as ~7.58 hours/frame with an image size of 64 × 64 pixels for SRUS imaging [23]. To solve the minimization problem efficiently for the rapid CS algorithm, alternative methods were introduced [24], [25]. This alternative method-based CS algorithm was applied for SRUS imaging of the rabbit kidney in vivo but utilized the constant noise parameter to solve the minimization problem [26]. The ℓ 1 -homotopy based CS (L1H-CS) algorithm showed the capability to search the optimized parameter naturally [27]. Fast super-resolution optical imaging of microtubules in BS-C-1 cells was performed based on the L1H-CS algorithm, which showed a ~300-fold faster computational time than the standard CS algorithm [28].
In this paper, we present an L1H-CS based SRUS imaging technique that could reduce both the computational time and the acquisition time of ultrasound image frames for the reconstruction of the SRUS image. The localization efficiency, error, and computational time of the proposed technique were evaluated by comparing it with the conventional localization technique using the two and randomly distributed point targets (PTs) generated in simulation. Furthermore, we developed clutter reduction, noise equalization, motion correction, and spatiotemporal noise filtering algorithms to perform in vivo imaging. The contrast to tissue ratio (CTR) and signal to noise ratio (SNR) were then analyzed quantitatively at ROI. Besides, we performed SRUS imaging of mouse tumor, kidney, and zebrafish dorsal trunk, and brain with the injection of high-density MBs. The number of frames and computational time required for the reconstruction of in vivo microvasculature images based on the proposed and conventional techniques were examined. Results demonstrate that the SRUS imaging technique could be substantially improved by applying L1H-CS algorithm for the localization of high-density MBs, suggesting its potential as a faster SRUS imaging technique.

A. ℓ 1 -homotopy Based Compressed Sensing
Mathematically, detected MBs signals, y, using the ultrasound imaging system have a linear relationship with MBs location, x y = Φx + e, (1) where y and x consist of the one-column vectors reshaped from the 2D matrix by row-wise concatenating the matrix of the original image or the localized results in the upsampled space or pixels [ Fig. 1(a)], respectively. The matrix Φ represents the point spread function (PSF) derived by an imaging system by simulation or actual measurements. e denotes a noise vector. The i th column of matrix Φ corresponds to the acquired image when only a one-point target exists at a particular index i of x. To achieve upsampled localized results, x, at a measured frame, y, the interior point method utilizes to resolve the ℓ 1 -norm minimization problem for CS as follows [21]: minimize: x 1 subject to : x i ≥ 0 and Φx − y 2 ≤ ε, (2) where ε denotes a constant showing the balance of fidelity and sparsity in the optimization algorithm. For example, a perfect fit for an original image is estimated when the ε value is 1, but ε is generally set somewhat larger than 1 to apply uncertainty to the variance measurement. For in vivo imaging we selected ε to be 2.3 [21]. In this paper, the interior point method-based CS algorithm was implemented by using the CVX optimization package in Matlab (CVX-CS) [29].
To solve the minimization problem efficiently, an alternative method was proposed as follows [24], [25]: where λ is a positive weight. Eq. (3) was expanded to build the L1H-CS formulation using the homotopy parameter, ϵ, as follows [27], [30]: As ϵ adjusts from 0 to1, the Eq. (4) steadily transforms to the Eq. (3). The homotopy method offers a typical framework that can solve an optimization problem by continuously transforming it into a related problem with less computing cost [27]. In this paper, we used L1H-CS package provided in https://github.com/sasif/L1-homotopy [27].
To reduce the computing cost of the CVX-CS, and L1H-CS algorithms, image patches with a kernel size of 7 × 7 pixels were used [thick black square in Fig. 1(a)]. Each image patch was extended by half pixels after upsampled by the factor of 8. Adjacent image patches were overlapped by 2 pixels to localize MBs by CS algorithms. Upsampled image patches were stitched patch by patch after taking 5 × 5 pixels to form SRUS patch [thick red square in Fig. 1(a)] to prevent false localization of MBs close to edges of image patches. The image patch-based CS algorithms allowed parallel computing in Matlab to improve computing speed substantially. The computation is performed on a personal desktop, with an AMD Ryzen 7 3700X 3.6 GHz processor and 16 GB RAM.

B. Ultrasound Imaging Configuration
For the acquisition of the in-phase and quadrature (IQ) data, we utilized an ultrasound imaging research platform (Vantage 256, Verasonics Inc., Kirkland, WA, USA) with a high-frequency linear array transducer (L35-16vX, Verasonics Inc., Kirkland, WA, USA). The transducer was incorporated with a 3-axis motorized stage (ILS150CC, Newport Corp., Irvine, CA, USA). A five-angle (from −7° to 7° with a step size of 3.5°) plane-wave imaging with an effective frame rate of 500 Hz and transmit frequency of 25 MHz was used. The excitation peak-to-peak voltage of 10 V with a one-cycle sinusoidal signal was generated for each transmission angle. A total of 500 compounded IQ data sets were acquired for 1 sec with pixel resolutions of 27.5 and 34.5 μm in axial and lateral directions, respectively. Furthermore, PSF at the elevation focus of ~ 8 mm was acquired using the Verasonics simulator. The PSF acquired in the simulation was utilized to evaluate localization techniques using two-PTs and randomly distributed PTs. For in vivo imaging, the PSF of MBs was generated as follows: where (z, x) and (c z , c x ) are the coordinate system and the center position of generated PSF.
The σ a and σ l represent the average full width at half maximum (FWHM) of MBs in axial and lateral directions, respectively. In this work, the average FWHM for each direction was measured by manually selecting 10 MBs at different frames. The average FWHM of PSF in axial and lateral direction was measured to be 39.6 and 48.7 μm, respectively.

C. Post-Processing Procedures for in Vivo Super-Resolution Ultrasound Imaging
The post-processing procedures for SRUS imaging in vivo are shown in Fig. 1(b). To remove the clutter signal, SVD based spatiotemporal filtering was applied to the IQ data. The SVD algorithm was briefly summarized as follows. First, the stacked spatial IQ data are converted into a 2D space-time matrix called Casorati matrix S. The SVD of S was described as follows [15]: where the column of U and V matrices are the associated spatial singular and temporal vectors, and * represents the conjugation transpose. Λ is a diagonal matrix that consists of singular values. In this paper, the SVD filtering was implemented using the "svd.m" in Matlab. To remove the clutter signals, the lower (T l ) value thresholds were set by the examination of the turning point on the singular value orders [31], [32]. Then, the filtered signal was reconstructed by: After removing the clutter signals, IQ data sets were converted to the intensity data by the square root of the sum of squares. Then, the noise equalization algorithm was applied to compensate for the ramp-shaped background noise which might occur due to the varying gain along with the depth [33]. The background noise was derived by the reconstructed signal using only the lowest singular value. Then, the noise equalized MBs images were extracted by dividing the background noise profile. The rigid motion of noise equalized MBs images were corrected. The rigid motion in axial and lateral direction throughout the frames at the ROI was estimated by applying the subwavelength motion correction algorithm with the first frame of clutter signals as a reference [34]. To further reduce the noise while preserving the MBs signal, we applied nonlocal means (NLM) filter to the spatiotemporal signals [7].
MBs were then localized by applying the CVX-CS, and L1H-CS algorithms to denoised MBs signals as described in the previous section. We also implemented 2DCC based localization technique to compare the performance of CS-based localization technique. Denoised MBs and PSF signals were spatially interpolated 8× under the cubic spline interpolation. 2DCC between the interpolated MBs and PSF signals was conducted. The centroid of MBs was found using the regional peak points finding algorithm ("imregionalmax.m" in Matlab). These processes were iterated for every frame. Each MB localized in each frame was tracked by pairing with the nearest MB in the following frames using the Hungarian algorithm (https://github.com/tinevez/simpletracker) [35], [36]. Finally, the tracked MBs were superimposed to reconstruct the final SRUS image. Furthermore, the Power Doppler (PWD) image was produced using the denoised MBs signal as follows [15]: To evaluate denoising techniques including SVD, SVD + noise equalization, SVD + noise equalization + NLM filtering, we measured CTR and SNR of denoised images as follows: SNR(dB) = 20log 10 V MBs std(noise) , (10) where V MBs and V Tissue are the signals of MBs and tissues, respectively. The std(noise) represents the standard deviation of the noise signal. Besides, to compare the localization techniques quantitatively, we defined the vessel density (VD) in the tumor as follows [37].
V D(%) = number of vessel pixels using n frames number of pixels in tumor × 100, Tumor areas were selected manually by examining the B-mode and SRUS images of mouse tumors.

D. Point Target Preparation Using Simulation
To evaluate the performance of localization using 2DCC, CVX-CS, and L1H-CS algorithms, two points aligned along with axial and lateral directions and randomly distributed PTs were generated in simulation. First, two-PTs were constructed at each inter-point distance from 13.75 to 110 μm with a step size of 13.75 μm in axial and lateral directions, respectively. The wavelength in this study was 61.6 μm. We localized these two points using 2DCC, CVX-CS, and L1H-CS algorithms to examine the minimum distance that can be localized in axial and lateral directions. Localized PTs were paired with the nearest PTs generated in simulation using Hungarian pairing algorithm with a limitation of a maximum pairing distance of 50 μm. Then, the localization efficiency and errors were identified by counting the localized PTs and calculating the distance deviated from the locations of simulated PTs. These were iterated 5 times for each density. Then, the average localization efficiency and error with the standard deviation of 2DCC, CVX-CS, and L1H-CS algorithms were measured (n = 5). Furthermore, we compared the computational time with an image size of 62 × 62 pixels at each density of PTs. The computational time was measured by using 'tic' and 'toc' MATLAB function.

E. In Vivo Imaging of Mouse Tumor Model and Normal Adult Zebrafish
The

A. Evaluation of the Localization Techniques Using Two-Point Targets in Simulation
The performance of the localization capability was analyzed by comparing 2DCC, CVX-CS, and L1H-CS based localization techniques. Two-PTs separated by a known distance in axial and lateral directions were used. Figs. 2(a-f) show the minimum distance between two-PTs that can be localized by 2DCC, CVX-CS, and L1H-CS algorithms, respectively.
These were overlaid with B-mode images. Blue, green, and red asterisks indicate the locations of localized points by 2DCC, CVX-CS, and L1H-CS, respectively. Black circles indicate the location of simulated points. Figs. 2(g-h) present the number of localized points depending on the distance between two points by 2DCC, CVX-CS, and L1H-CS. Because only two-PTs were used, the maximum number of localizations was 2 2DCC localizes two simulated points correctly when the distance between two points becomes 82.5 μm in the axial direction and 96.3 μm in the lateral direction. CVX-CS can localize two points when they are 41.3 μm apart from each other in axial and lateral directions. L1H-CS algorithm can localize two points when the separation is 27.5 μm in both directions.

B. Evaluation of the Localization Techniques Using Randomly Distributed Point Targets in Simulation
We compared the performance of the localization techniques using randomly distributed PTs depending on point density. Figs. 3(a-c) Fig. 3(d)]. CS-based algorithms followed the black dashed line closely until the target density increased to 55.62 PTs/mm 2 as indicated by the black arrow in Fig. 3

C. Evaluation of Post-Processing Techniques Using Images of in Vivo Mouse Tumor Model
After evaluating the localization technique using simulated PTs, we first examined the effects of combinations of post-processing techniques including SVD, noise equalization, and NLM filtering to improve CTR and SNR of MBs signals. Figs. 4(a-d) show images before and after applying SVD, noise equalization, and NLM filtering to the data set acquired by ultrasound imaging of in vivo mouse tumor model. To calculate CTR, MBs and tissue signals were collected from the B-mode, SVD, SVD + noise equalization, and SVD + noise equalization + NLM filtered images. ROIs were selected at the same depth indicated as the orange boxes shown in Fig. 4(b) for all images. CTR of B-mode, SVD, SVD + noise equalization, and SVD + noise equalization + NLM filtering were measured to be -0. 23, 12.84, 13.83, and 15.1 dB, respectively, as shown in Fig. 4(e). We measured SNR of SVD, SVD + noise equalization, and SVD + noise equalization + NLM filtering along the scan line indicated by the orange dotted line in Fig. 4(b). Fig. 4(f) shows the intensity profiles along the same scan line shown in Fig. 4(b) images after each step as shown in Figs. 4(b-d). SNR improves from 21.02 to 27.35 dB as noise equalization and NLM filtering were applied as shown in Fig. 4(g).
Using denoised data sets, we reconstructed PWD and SRUS images using 2DCC, CVX-CS, and L1H-CS algorithms as shown in Figs. 5(a-d). As indicated by the white arrows in Fig.  5(a), smaller vessels were visualized efficiently by SRUS imaging [ Fig. 5(a-d)]. We further examined the FWHM of vessels along the orange dotted line. As the black dotted line indicates in Fig. 5(e) We further examined the effect of the number of frames by measuring the VD in the tumor. The SRUS image reconstructed by the L1H-CS algorithm using the 500 frames was overlaid with the B-mode image as shown in Fig. 6(a). The tumor area was indicated by the white dotted circles shown in Fig. 6(a). The maximum VD inside the tumor area of the SRUS image was reconstructed by L1H-CS [ Fig. 6(a)] with 500 frames were measured to be 13.76%. We also measure the maximum VD within the same region of images reconstructed by 2DCC and CVX-CS, which are 4.08% and 9.82%, respectively as shown in Fig. 6(b).
To calculate the computational cost to acquire the similar quality of vessel images between three algorithms, we determined the number of frames to be used at VD of 4.08% as shown in Fig. 6(b). 500 frames, 125 frames, and 75 frames were used for the reconstruction of SRUS images using 2DCC, CVX-CS, and L1H-CS techniques as shown in Figs. 6(c-e) to achieve VD of 4.08%. Computational times were measured to be approximately 4 min (2DCC), 59 min (CVX-CS), and 3 min (L1H-CS), respectively, as shown in Fig. 6(f).

D. In Vivo Mouse and Zebrafish Super-Resolution Ultrasound Imaging Using ℓ 1 -homotopy Based Compressed Sensing Algorithm
Figs. 7(a-d) represent reconstructed SRUS image overlaid with B-mode images of mouse tumor and kidney and zebrafish dorsal trunk and brain, respectively. We measured the vessel size and distance of adjacent vessels at the ROI along the yellow lines in Figs 7(a-d). In SRUS images of mouse tumor, the vessel sizes were measured to be 25.68 μm in profile 1, and the distances between adjacent vessels were measured to be 109.9 and 109.01 μm in profiles 2 and 3, respectively [ Fig. 7(e)]. The vessel sizes were measured as 32.2, and 22.2 μm in profiles 1, and 2, and the distance of adjacent vessels was estimated as 149.9 in profile 3 of the SRUS image of mouse kidney [ Fig. 7(f)]. Besides, the vessel sizes were measured as 55.7, and 85.85 μm in profile 1, 64.8 μm in profile 2, and 61.8 μm in profile 3, respectively, in the SRUS image of zebrafish dorsal trunk [ Fig. 7(g)]. The vessels size of 63.23 and 64.2 μm in profiles 1, and 2 were detected in the SRUS image of the zebrafish brain [ Fig. 7(h)]. The smallest vessel size is 22.2 μm by SRUS imaging. Considering the half wavelength of 25 MHz ultrasound as 32 μm, the diffraction limit was overcome by L1H-CS based SRUS imaging.

IV. Discussions
Recently, CVX-CS based SRUS imaging technique was introduced with a single plane wave transmission [22]. The validation was performed by using simulated high density MBs. In our paper, we significantly reduced the computational time by using L1H and evaluated the performance of L1H-CS with in vivo experiments. SRUS images of the vasculature of rabbit kidney was reconstructed by using a fast iterative shrinkage thresholding algorithm (FISTA) based compressed sensing technique [26], [39]. The computational cost of FISTA is similar to L1H per iteration to optimize Equation (4); however, FISTA requires more iterations to reach approximated solutions [40]. Therefore, we developed L1H-CS based SRUS imaging technique to improve computational cost and localization efficiency of densely populated MBs. SRUS imaging based on a modified sub-pixel convolutional neural network (mSPCN ULM) was proposed to image vasculature of rat kidneys [23]. In this experiments, mSPCN ULM took approximately 23 seconds per one frame using GPU acceleration, while L1H-CS takes approximately 2.4 seconds per one frame. While quantitative comparison between mSPCN-ULM and L1H-CS based SRUS imaging technique may not be adequate due to different computing resources and setup for data acquisition and processing, these results support a potential of L1H-CS to improve current SRUS imaging.
The localization results of the two-PTs and randomly distributed-PTs demonstrated that CS-based algorithms improved localization efficiency compared to 2DCC algorithm. L1H CS algorithm exhibited slightly better performance than the CVX-CS algorithm to localize the PTs. In particular, CVX-CS localized two-PTs when the axial and lateral distances between two points approach 41.25 μm. L1H-CS algorithm localized two points when the axial and lateral distances approach 27.5 μm [ Fig. 2]. Optimization of ϵ in Equation (4) improved localization efficiency of L1H-CS algorithm, while ε in equation (2) was fixed for CVX-CS. Randomly distributed PTs with a density of 55.62 PTs/mm2 were localized to be 50.74 ± 0.62 PTs/mm2 by L1H-CS, which demonstrated 10% localization accuracy compared to CVX-CS. One possible explanation of improved localization efficiency of L1H-CS compared to CVX-CS is that L1H-CS algorithm searches optimized regularization parameter ϵ, while CVX-CS fixed the constant ε as 2.3 [Eqs. (2) and (4)] [21]- [23]. Besides, L1H-CS required a much lower computing cost than the CVX-CS algorithm [ Fig. 3(f) and Fig. 5(f)]. This is because the computational complexity of L1H-CS is no longer affected by the inverse matrix, which is a major factor that demands the high computational cost in CVX-CS algorithm [28].
Total computational time of L1H-CS is higher than that of 2DCC as shown in Fig. 5 L1H-CS localization technique proposed here can detect densely populated PTs or MBs efficiently as demonstrated in Fig. 2 and Fig. 3. Thus, we performed the SRUS imaging of mouse and zebrafish by injecting the high-density MBs (0.75 -2.0 × 10 9 MBs/mL) than the standard density of MBs (1.5 -5.6 × 10 8 MBs/mL). In the in vivo experiment, VD in tumors was measured as 4.08% based on microvessel image reconstructed by 2DCC, while L1H-CS could reconstruct microvessel of 13.76% in tumor regions using 500 frames [ Fig.  6(b)]. These resultssuggestthattheL1H-CSalgorithmallowsSRUSimaging with the injection of high-density MBs, suggesting that the acquisition time of ultrasound image frames can be reduced, resulting in a faster imaging session would be implemented.
Based on the results demonstrated in this study, L1H-CS can be used to decrease imaging acquisition time to acquire a similar quality of SRUS images [ Fig. 6] or improve the image quality of SRUS with significantly improved computational cost [ Fig. 5]. L1H-CS only requires 75 frames to reconstruct the similar quality, validated by VD, projected acquisition time is only 15% of that of the typical 2DCC method. To acquire high quality SRUS images, validated by VD, 500 frames are used with 5 times higher computational cost. However, the computational cost of L1H-CS was decreased by 5 times compared to CVX-CS algorithm. This exhibits the potential to be used in clinics.
PWD image could not represent microvessels indicated by white arrows as shown in

V. Conclusion
In conclusion, we described the L1H-CS based SRUS imaging technique that could reduce computational time and acquisition time or significantly improve SRUS image quality with some expenses of post-processing cost. In the simulation, we confirmed that a minimum distinguishable distance of PTs was 27.5 μm, and efficient localization of PTs compared to a conventional localization technique. By in vivo mouse imaging, the minimum detectable blood vessel size in mouse kidney was measured to be 22 μm, which confirms the sub diffraction limit capability of L1H-CS. Based on results both in simulation and in vivo experiments, the proposed algorithm significantly improves SRUS image quality and data acquisition time. (a) For compressed sensing-based reconstruction of SRUS image, upsampled image patches (thick black square) were used by overlapping 2 pixels between adjacent image patches to remove false localization of point targets near boundaries and to correctly build SRUS patches (red thick square). By stitching SRUS patches side by side, a complete SRUS image was acquired from B-mode image (black dashed square). (b) The schematic picture describes overall post-processing steps to reconstruct SRUS images and the administration of microbubbles to mouse tumor models and zebrafish.