Time Efficient Ultrasound Localization Microscopy Based on A Novel Radial Basis Function 2D Interpolation

Ultrasound localization microscopy (ULM) allows for the generation of super-resolved (SR) images of the vasculature by precisely localizing intravenously injected microbubbles. Although SR images may be useful for diagnosing and treating patients, their use in the clinical context is limited by the need for prolonged acquisition times and high frame rates. The primary goal of our study is to relax the requirement of high frame rates to obtain SR images. To this end, we propose a new time-efficient ULM (TEULM) pipeline built on a cutting-edge interpolation method. More specifically, we suggest employing Radial Basis Functions (RBFs) as interpolators to estimate the missing values in the 2-dimensional (2D) spatio-temporal structures. To evaluate this strategy, we first mimic the data acquisition at a reduced frame rate by applying a down-sampling (DS = 2, 4, 8, and 10) factor to high frame rate ULM data. Then, we up-sample the data to the original frame rate using the suggested interpolation to reconstruct the missing frames. Finally, using both the original high frame rate data and the interpolated one, we reconstruct SR images using the ULM framework steps. We evaluate the proposed TEULM using four in vivo datasets, a Rat brain (dataset A), a Rat kidney (dataset B), a Rat tumor (dataset C) and a Rat brain bolus (dataset D), interpolating at the in-phase and quadrature (IQ) level. Results demonstrate the effectiveness of TEULM in recovering vascular structures, even at a DS rate of 10 (corresponding to a frame rate of sub-100Hz). In conclusion, the proposed technique is successful in reconstructing accurate SR images while requiring frame rates of one order of magnitude lower than standard ULM.


I. INTRODUCTION
U LTRASOUND is recognized as one of the principal medical imaging modalities for visualizing subcutaneous body structures such as organs, muscles, and arteries due to its safety, cost-effectiveness, and non-invasiveness [1].However, conventional ultrasound performance is constrained by the diffraction limit, which constrains the spatial resolution to the scale of the wavelength [2], [3].While at higher frequencies the spatial resolution improves, the beam's penetration depth diminishes, thus causing a trade-off between spatial resolution and penetration depth [2].
In the last decade, Ultrasound Localization Microscopy (ULM), namely a technique for visualizing the vasculature at unprecedented resolution, has addressed the aforementioned trade-off [4].The super-resolved (SR) images obtained by ULM provide indispensable information for the understanding and diagnosis of various diseases that modify the vascular structure or blood flow, such as cancer [5], [6], [7], [8], stroke [8], [9], and atherosclerosis [6].ULM operates by individually localizing and tracking intravenously injected microbubbles (MBs), also known as ultrasound contrast agents (UCAs).The centroid of each MB can be localized with an accuracy 10 times lower than the wavelength [10], when UCAs are well separated.ULM is thus able to bypass the conventional ultrasonic diffraction limit by aggregating the centroid of single echoes produced by the UCAs [2], [11].Finally, SR images are obtained by tracking the microbubbles as they flow through the vascular system, providing essential information about the vascular bed's structure as well as its flow dynamics [12].
Despite the advantages that ULM may bring in the clinical setting, its utilization is currently limited by two main practical considerations [2], [10], [13].Firstly, low MBs concentrations and long acquisition times are necessary to record sufficiently separated MB echoes for localization and tracking [13].Secondly, a high frame rate is necessary to coherently pair the MB signals from one frame to the successive one [14], [15], [16].To date, this need for long acquisition times and high frame rates has hampered the ULM clinical translation.
Extensive research has been conducted to reduce the ULM acquisition time.One straightforward solution is to increase the concentration of MBs [17].High concentrations, on the other hand, result in MB overlapping signals, making it impossible to discriminate between different MBs' echoes and individually localize them.To tackle this issue, various research works have been conducted, e.g. based on sparse recovery methods [18], [19], deconvolution [20], [21], and velocity filtering [2].To distinguish overlapping MBs signals, Huang et al. [22] proposed establishing different subpopulations based on different flow patterns.However, these techniques require high frame rates that are generally not possible with the available clinical scanners.Other approaches take advantage of deep learning architectures to localize MBs in dense scenarios [23], [24].However, deep learning necessitates a significant amount of labeled data for model training and testing, which is difficult to obtain.
Moreover, Guo et al. [15] evaluated by a quantitative investigation the loss of ULM performance when decreasing the frame rate.The study demonstrates, by utilizing an in vivo dataset, that a decreased frame rate results in an unreliable microbubble detection and tracking, which in turn leads to a decreased spatial resolution (from 9.9 µm to 15 µm when going from 1kHz to 250Hz) as well as a decreased saturation (saturation drop above 50% from 1kHz to 100Hz).Furthermore, the frame rate has a significant impact on the precision of the velocity measurements.
To tackle the high frame rate requirement, Ackerman et al. [25] and Opacic et al. [26] implemented a modified Markov chain Monte Carlo data association to recover UCA tracks at a lower frame rate, which optimally associate MBs positions based on an underlying motion model.Another approach [27] uses a Kalman filter as a tracking backbone, and introduces MB image features on top.Although these approaches can recover UCAs tracks at a lower frame rate, their improved performance comes at a significant computational cost.Enabling lower frame rate acquisitions using TEULM would alleviate some of the data size constraints for ULM and would make the technology more pragmatic for clinical translation, where low frame rates are standard.To be more precise, adopting lower frame rate acquisitions in ULM can bring about numerous advantages in the clinical context.These benefits encompass reduced data size, compatibility with standard equipment, and cost-effectiveness.
In this study, we aim at bridging the gap between the high frame rates required by ULM (in the range of 1kHz) and the limited frame rate of the standard clinical scanners (in the sub-range of 100Hz).
To this end, we propose a new time-efficient ULM (TEULM) pipeline based on a novel 2-dimensional (2D) interpolation technique to relax the requirement on the high frame rate necessary to obtain super-resolved images.The article is structured as follows.Section II gives a detailed explanation of the interpolation technique used upstream the ULM pipeline.Next, we describe the key steps of ULM framework, as well as the TEULM pipeline.In the result section (Sec.III), we quantify the quality of the reconstructed images using TEULM, and compare these images with the SR images obtained through standard ULM.Finally, we provide the discussion and conclusions of the proposed method in Sections IV and V, respectively.

A. RBF-Based 2D Interpolation
Radial Basis Functions (RBFs) have been frequently used in recent years to interpolate high-dimensional sparse data [28],  [29].The main advantage of this technique is that it delivers respectable numerical results without entailing complicated computations.This makes it simple to apply to complicated geometries in a variety of spaces.In the following, we provide the formulation of the 2D interpolation problem using RBFs.
Given the data (X j , Y j ), j = 1, 2, . . ., N , where X j and Y j ∈ R d (in our problem d = 2), we aim to find a continuous function S(X) so that S(X j ) = Y j ∀ j ∈ {1, 2, . . ., N }.In addition, the interpolation bases are chosen as {φ (X − X 1 ) , φ (X − X 2 ) , . . ., φ (X − X N )}, where X 1 , X 2 , . . ., X N are predetermined interpolation nodes.These assumptions lead to the definition of RBFs as: where ∥•∥ 2 indicates the ℓ 2 norm (i.e., the Euclidean distance) and τ represents the so-called shape parameter.Due to their radial dependence, these functions are referred to as RBFs.Some possible RBFs are listed in Table I.
Following that, we develop a linear combination of RBFs φ(r ) to move forward with the 2D interpolation: where X = (x 1 , x 2 , . . ., x d ), τ is the shape parameter, β j is the interpolation coefficient, and r = ∥X∥ 2 .If we intend to estimate f : R d → R utilizing the interpolators determined by Eq. ( 2), we need to find the interpolation coefficients β j .To do that, we use the following interpolation condition: By applying Eq. (3) at the central points of X c j for j = 1, . . ., N , we obtain the linear system: which can be expanded as follows: where X j indicates the imaging pixels, f (X j ) refers to the pixel intensity of the imaging point X j , r i, j = ∥X c i − X c j ∥ is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the Euclidean distance, and N is the number of points in the 2D grid (i.e., the imaging points in our problem).Based on the system in Eq. ( 5), to solve the interpolation problem (i.e., finding the β values), we must invert the coefficient matrix , as discussed in the next section.Now, using the coefficients calculated from the above system, if we seek to analyze the interpolator at M points X i , i = 1, 2, . . ., M, the new coefficients matrix entries must be constructed as follows: h i j = φ(∥X i −X c j ∥ 2 ) ∀ j ∈ {1, . . ., N }.Finally, the following matrix multiplication is used to evaluate the interpolant at the M points: where F est represents the estimated F (i.e., in our interpolation problem, F est indicates the interpolated data).
The invertibility of the interpolation matrix and the error analysis are described in Appendix A and B, respectively.

B. Proposed Approach
1) Standard ULM Framework: As shown in Fig. 1, the traditional ULM framework consists of five basic steps: data acquisition, pre-processing and filtering, MB detection, MB localization and tracking, and accumulation.For the data acquisition step, the standard framework demands data acquired at high frame rates.Following that, the high frame rate data is pre-processed to reduce clutter and tissue noise.In this study, we employ singular value decomposition (SVD) spatio-temporal filtering to distinguish the MBs signal from the surrounding tissue [30].The SVD filter takes advantage of the fact that back-scatter tissue signals have a high spatio-temporal coherence compared to MB echoes.Therefore, by eliminating the initial k singular values, the tissue signal is filtered out from the data.
Next, MB signals are detected based on the pixels intensity in each frame.As in [17], an estimate of the number of MBs (i.e., n) in each frame is utilized instead of creating a threshold to identify the minimum pixel intensity value representative of the MB signals, which is non-trivial and may result in poor MB detection.In this paper, for each frame, the n pixels characterized by the most intense values are selected.Individual MBs are then localized using the Gaussian fitting algorithm [17], [31], [32], which fits a 2D multivariate Gaussian distribution.This approach assumes that the maximum intensity value (I ) in a grid of [N z , N x ] pixels corresponds to the MB centroid.To recover the super-resolved location, the original image is up-sampled by a resolution factor r es.Each pixel is interpolated into a r es × r es grid, allowing to localize the MBs with a precision r es times higher than the pixel size.The final image is then obtained as a grid of size [N z ×r es, N x ×r es].Here, we define r es as the empirical limit found in [17], namely 10.The localization problem can then be mathematically described as an optimization problem.After centering the subspace in the maximum value of the subsampled grid, the super-resolved z and x coordinates, (c z , c x ), of each MB is evaluated as the location of the Gaussian peak as follows: where f z and f x respectively indicate the full width at half maximum (FWHM) along z and x directions.The standard deviations σ z and σ x are empirically evaluated to be equal to 1, as in [17].Then, the bipartite tracking algorithm matches each detected MB location in a frame-to-frame fashion.Thanks to their simplicity, the bipartite graph-based pairing methods are well established in the ULM context [13], [14], [15], [16], [17], [22], [27].These algorithms solve the tracking problem through an optimal assignment, based on the MBs frame-toframe distance.In particular, the relative distance between all the UCAs in two successive frames is evaluated.The MBs are then considered to be the same, or matched, in such a way that the total distance between all the MBs in the two frames is minimized.Since the assignment depends on the frameto-frame distance, it is easy to understand the importance of a high frame rate to accurately pair MB signals [15], [16].This paper implements a modified Khun-Munkres bipartite tracking algorithm, as in [33].Generally, the method works on the assumption that each MB in one frame is paired with one MB in the subsequent frame.However, due to the uncertainty originated by the arbitrary estimation of the number of MBs in each frame in the localization step, the assumption may result in a sub-optimal solution.Furthermore, a sub-optimal solution might be caused by unpredictable UCAs events, such as MB appearance and/or disappearance.
As such, the Khun-Munkres algorithm is modified to cope with partial assignment.In particular, a MB in a frame f which has depleted all the possible assignment options in the next frame f + 1 is discarded [14].Moreover, we establish a threshold for the minimum length of a track in order to further strengthen the robustness of the tracking performance, i.e., a track is discarded if it is lower than the threshold because it is considered as noise.
Finally, an SR intensity map is obtained by accumulating the tracks, previously smoothed and interpolated to fill in the missing points, over all frames.The velocity map, which defines a MB's velocity by its displacement over time, can also be retrieved and superimposed onto the SR intensity map.
2) TEULM Framework: TEULM introduces, in the preprocessing stage, the RBF-based 2D interpolation technique to relax the frame rate required by standard ULM.The interpolation method aids in recovering the lost data, which is the result of adopting a reduced frame rate acquisition.More specifically, we propose sparse (low frame rate acquisition) combined with 2D interpolation to reconstruct the unknown information corresponding to the frames that are not recorded in the low frame rate acquisition (i.e., for a given downsampling (DS) rate, interpolating DS -1 middle frames to achieve high frame rate IQ data needed for high-quality SR imaging).Following that, we adopt the standard ULM framework, described in Section II-B.1.We pre-process the interpolated data by applying SVD filer first, keeping Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II DATASETS SUMMARY IN THIS STUDY
the k-values used for filtering as the original ones, and successively a Butterworth filter.Then TEULM's framework comprises of detection, localization and tracking of the MBs.The final images are then obtained through accumulation (see Fig. 1).ULM and TEULM parameters are summarized in Table III.In the case of using ULM at the original frame rate and TEULM, we use the same parameters used in [17].However, for the ULM at down-sampling rates > 1 (non-interpolated data), we adapt the maximum linking distance and minimum track length by a factor that is proportionate to the level of down-sampling.For example, when down-sampling by a factor of 2 (DS=2), the maximum linking distance is adjusted to be twice the original linking distance, while the minimum track length is set to half of its original value.This adjustment is made for a fair comparison, which accounts for the fewer frames used yielding in larger gaps between microbubble positions in consecutive frames.Similarly, due to the reduced number of frames, a shorter minimum track length is necessary.

C. Data
To assess the proposed strategy, we use four in vivo datasets.We denote the in vivo datasets as dataset A for the Rat Brain, dataset B for the Rat Kidney, dataset C for the Rat tumor and dataset D for the Rat Brain bolus.
The datasets, presented in [17] and summarized in Table II, are recorded with a 15 MHz center frequency probe (Vermon, France).We opted for the use of in vivo datasets due to their capacity to encompass a wide range of hemodynamic structures across different organs (such as the brain and kidney) and various injection modalities.

D. In Vivo Evaluation Metrics
We use the following metrics to evaluate the TEULMgenerated SR image quality and compare it with that obtained by standard ULM.
1) Root Mean Square Error (RMSE): The similarity between the original image, reconstructed at 1KHz (I 1 ), and the images interpolated at different DS rates (I DS ), with DS = 2, . . ., 10, is assessed pixel-by-pixel using the root mean square error: where N pi x denotes the total number of pixels in the reconstructed image.Both density maps, which depict the vascular anatomy, and velocity maps are evaluated.The Root Mean Square Error (R M S E) is employed to gauge the resemblance in terms of vascular structure and intensity among the reconstructed images acquired through TEULM and ULM at various levels of DS rates.Indeed, the R M S E conducts a pixel-level assessment of the error, taking into account the image reconstructed from the initial high frame-rate data (1KHz for datasets A, B, and D, and 500Hz for dataset C) as the reference standard.
2) DICE Score: The D I C E score evaluates the similarity between two images by measuring the overlapping area of two images over the area of their union [34]: The D I C E score quantitatively assesses the resemblance of the reconstructed binary images generated by ULM and TEULM across various DS rates.
3) Percentage of Preserved Tracks: Since the number of detected tracks has a significant impact on the saturation and on the quality of the reconstructed images, we analyze the statistics related to the number of tracks at various DS rates.The ratio of the number of tracks found at the various DS rates (DS = 2, . . ., 10) to those found at DS = 1 is used to calculate the percentage of preserved tracks.The percentage of preserved tracks is evaluated for density maps only.In ULM, a density map refers to a spatial representation of vessel density within a region of interest.
4) Saturation: We further assess the TEULM performance at different DS rates through saturation (i.e., indicating the number of pixels crossed by MBs) [15], [16].The evaluation of the saturation is meaningful since the reconstructed images are the result of the accumulation of the pixels crossed by the MBs.In this study, the saturation curve is calculated as the percentage of filled pixels over the total area of the image.The saturation is evaluated only for density maps.
5) Fourier Ring Correlation (FRC): Fourier Ring Correlation is a commonly used method in the ULM context to estimate the spatial resolution of the super-resolved images [35].The F RC evaluates the spatial frequency correlation along isospatial frequency rings (R) between two sub-images in the frequency domain, formed by random splitting the tracks: The resolution is the inverse value of the F RC curve, when the F RC drops below the threshold value.In this study, we select the 1/2 bit as thresholding method, which corresponds to the highest spatial frequency necessary to collect the information to fill half a bit [36], [37].To achieve an unbiased and objective quantification, we conduct global (i.e., whole image) F RC resolution estimates.

III. RESULTS
In this section, we present the results obtained through TEULM using the in vivo datasets described in section II-C.
We provide the in vivo results of standard ULM on the original high frame rate data (namely, 1kHz for dataset A, B and D and 500Hz for dataset C), see Fig. 2.These images are considered as a reference.Figures 3, 4, 5, and 6 depict the ULM and TEULM reconstructed density and velocity maps for datasets A, B, C, and D, respectively, across different DS rates.It's worth emphasizing that for dataset C, we provide the results until down-sampling 6 (which implies a frame rate of 83Hz) which demonstrates the capability of our proposed TEULM approach to the frame rates close to the current clinical system's frame rate.In addition, it should be noted that conventional ULM fails to reconstruct any tracks using a higher down-sampling rates (DS ≥ 7).The results show that for DS ≥ 2, TEULM performs better than standard ULM in recovering the local and global information (i.e., ULM even fails to recover the global information at higher DS rates).As expected, a decrease in the frame rate (from left to right) results in a proportional degradation in the image quality.Qualitatively, the vessel definition appears to be less continuous and more blurry at higher DS rates compared to the original images, see Fig. 2. We quantitatively assess this loss of definition in TEULM's generated density maps through the spatial resolution obtained using F RC on the entire image (Table IV), for all the in vivo datasets.The Fourier Ring Correlation evaluates the spatial resolution as the correlation between 2 sub-images formed by randomly splitting the tracks forming the density maps.In addition, we evaluate the consistency of the density maps generated by TEULM and ULM with the ones obtained using the high frame rate data using the DICE score, see Table V.For all the in vivo datasets, we further evaluate the efficiency of the proposed technique by means of two of the previously introduced metrics: the percentage of preserved tracks, and the saturation.The metrics are shown in Fig. 8 considering downsampling 2, 4, 8 and 10.Fig. 8 clearly shows that TEULM preserves the majority of the tracks, as well as a constant saturation level compared to ULM, for which both the metrics decrease exponentially with increasing down-sampling rates.When under-sampling, we can generate multiple datasets from the high frame rate data.To assess the consistency across different datasets, we generate the permutation starting from frame 6 at down-sampling of 10 and compare the results with respect to the dataset generated starting from frame 1.The choice of the permutation allows us to evaluate the performance of the proposed approach for the most asynchronous couple of data that we could generate.Particularly, we assess the performance of TEULM through the statistical results (mean ± standard deviation) of the D I C E score, the saturation percentage, and the percentage of preserved tracks.The results of the DICE score are presented in Table V.The percentage of preserved tracks for dataset A, B, C, and D is 0.98 ± 7 × 10 −4 , 0.77 ± 4 × 10 −2 , 0.37 ± 1 × 10 −2 , and 0.83 ± 2 × 10 −2 , respectively.In addition, the saturation percentage reaches the value of 0.96 ± 2 × 10 −2 , 0.68 ± 1 × 10 −2 , 0.77 ± 3 × 10 −2 , and 0.9 ± 1 × 10 −2 for dataset A, B, C, and D, respectively.The results certify the consistency of TEULM across permutations.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.In addition, we analyze the influence of the frame rate on the velocity maps obtained with TEULM and ULM.To this end, we calculate the flow direction and the velocity magnitude maps for the in vivo datasets.For the rendering of the flow direction of the blood flow, we use the red and blue colors to encode the different directions (see second row of Figs. 3, 4, 5, and 6, for dataset A, B, C and D respectively).About the mean velocity magnitude encoding (see the third row of Figs. 3, 4, 5, and 6, for dataset A, B, C and D respectively), the slowest velocities are encoded in blue, whereas the high velocities are encoded in red.The results of the velocity profiles, provided in row two and three of the Figures 3 (dataset A), 4 (dataset B), 5 (dataset C), and 6 (dataset D), clearly show that there is a lower limit to the frame rate reduction under which it is not possible to recover the global and local velocity information of the vascular structure.
Moreover, we quantify the error for both the intensity and velocity maps of the sparse ULM and TEULM with respect to the original ULM (i.e., obtained with a frame rate equal to 1kHz for dataset A, B and D and frame rate equal to 500Hz for dataset C) utilizing the R M S E metric (Figure 7).The results allow us to evaluate the trend of these two metrics with respect to different DS rates, up to 50Hz. 1  Finally, we have conducted an evaluation of the mean velocity profiles for all the in vivo datasets (see Fig. 9).These profiles represent the mean velocity values (in mm/s) against the number of tracks associated with each mean velocity across all the tracks within their respective datasets.Fig. 9 shows that the mean velocities are misestimated at higher down-sampling, confirming that recovering the dynamic profile is a challenging task for the RBFs interpolators.

IV. DISCUSSION
The primary goal of this study is to present an innovative approach to enable ULM at low frame rates, which will facilitate ULM clinical translation.In this regard, we propose a modified ULM framework, called TEULM, that can cope with low frame rate data by employing the suggested 2D interpolation technique.
The proposed technique has been tested on four pre-clinical in vivo data (dataset A, B, C and D).For the in vivo datasets, 1 Intensity maps available in the link.
we demonstrate the high efficiency of TEULM in recovering the information lost when using a low frame rate (Figures 3,  4, 5 and 6, for datasets A, B, C and D respectively).More specifically, while both TEULM and ULM face an image quality loss at lower frame rates, as we can see in Figures 3,  4, 5 and 6, the extent of the loss varies significantly between the two techniques.In other words, the vessel's appearance for ULM degrades quickly.For instance, at a frame rate lower than 250Hz, the vascular structure is almost unrecognizable for all the datasets, which shows the ULM's low performance at lower frame rates.In contrast, when using TEULM, the overall vascular structure remains unchanged.Besides, an increase in the DS rate results in a loss of vascular definition, manifested by more blurry areas, mostly in the fine vasculature (See an example shown in the magnified Region of Interest (ROI) in Fig. 4 and highlighted by the green box.This loss of definition is quantitatively assessed by the spatial resolution estimates using the F RC (Table IV).Despite the loss in spatial resolution, the saturation and the percentage of preserved tracks are invariant in the reconstructed images using the suggested TEULM method (see Fig. 8).Further, for dataset B, we observe the appearance of vascular structures in TEULM at DS = 10, i.e., 100Hz, (Fig. 4 highlighted by the yellow arrow) that are absent in the original frame rate, i.e., 1K H z (see Fig 2).It should be noted that, in the absence of the ground truth, it is difficult to judge whether the appearing vascular structures are artifacts or real vessels.Nevertheless, we observe that for the dataset B the global bloodstream structure in the proposed TEULM is not substantially impacted by either the loss of definition or the emergence of vascular artifacts, resulting in a D I C E score of 0.81 at DS = 10.In Table V, we provided the D I C E score values for the in vivo datasets, which allows to quantitatively evaluate TEULM's performance in recovering the vessel's structure at a lower frame rate.The minimum D I C E score values (at DS = 10) obtained using TEULM are 79%, 81%, 68% and 82% for dataset A,B,C and D, respectively.TEULM's high D I C E scores show the effectiveness of the suggested technique in reconstructing the missing information.On the other hand, ULM demonstrates poor D I C E similarity indices (19%, 14% and 52% for dataset A, B and D respectively), which could be justified by the failure of the tracking algorithm at lower frame rates.Thus, for ULM, in regions characterized by complex structures (e.g.crossing vessels or junctions) the algorithm falsely pairs MBs, resulting in unreliable vessel depiction.This assumption would also explain why dataset B and C has a much lower similarity compared to dataset A and D. It seems in fact that the vasculature of the kidney and tumor (dataset B and C) is much more tortuous and characterized by short vessels, which are easily lost at lower frame rates.
The reason behind ULM's poor performance at low frame rates can be attributed to the tracking algorithm, which significantly relies on the frame rate as it employs a frameby-frame matching approach.In other words, with a decrease in the frame rate, MBs may be discarded due to the inability to find the correct match in successive frames.Consequently, several tracks are not found or rejected because of an insufficient amount of known information, leading to  a decrease in saturation and vessel definition.In fact, the number of preserved tracks in ULM for both datasets decreases exponentially as we increase the DS rate.As can be seen in Figs. 8, which presents the saturation curve and the percentage of recovered tracks at different frame rates for datasets A, B, C and D, in comparison to the original image.Additionally, at the lower frame rates, the saturation curves clearly show a decline in the ULM's capability to reconstruct the vasculature (i.e., at a DS rate of 5, the saturation on all the in vivo datasets decreases by more than the 50%).
TEULM, on the other hand, maintains a high percentage of recovered tracks.The value remains consistently around 1 for Fig. 4. ULM-and TEULM-generated density maps (first row), flow direction map (second row), and velocity magnitude map (third row) for Dataset B at different frame rates of 500Hz (DS = 2), 250Hz (DS = 4) and 100Hz (DS = 10).The color bar of the figure is the same as original ones in Fig. 2. Fig. 5. ULM-and TEULM-generated density maps (first row), flow direction map (second row), and velocity magnitude map (third row) for Dataset C at different frame rates of 250Hz (DS = 2), 125Hz (DS = 4) and 83Hz (DS = 6).We have presented images up to a DS rate of 6 because, at higher DS rates, ULM is unable to recover any tracks.The color bar of the figure is the same as original ones in Fig. 2.
dataset A and approximately 0.85 for dataset D across all DS rates.However, we do notice a reduction in the percentage of preserved tracks at DS = 6 for dataset B and DS = 4 for dataset C. Nevertheless, TEULM's ability to preserve the number of tracks is reflected by the saturation of the density maps.The saturation levels for datasets A, B, and D remain consistently near 95%, while dataset C exhibits a saturation level of around 80%, when compared to the SR image generated using the original high frame rate data.The number of preserved tracks and saturation percentage indicate that datasets B and C are more challenging for TEULM.There might be a physiological reason for this.In fact, dataset B relative to the Rat Kidney and dataset C representative of the Rat Tumor are characterized by shorter and broken tracks, which are more difficult to recover when using a higher DS rates.
TEULM and ULM super-resolution density and velocity maps are further evaluated through the R M S E at different DS = 2, 4, 8, 10 rates (shown as a heatmap in Fig. 7).In both TEULM and ULM, the error, which is color-coded, exhibits an upward trend as the down-sampling rate increases.When examining TEULM's density maps R M S E across all the in vivo datasets, it is noteworthy that the error gradient gradually rises.Remarkably, it appears that this gradient has  not yet reached a plateau.Specifically, the R M S E remains below 50% for all datasets when employing TEULM, whereas it surpasses 80% when using ULM.This observation implies the potential for further increasing the DS value.In contrast, the error gradient for ULM stabilizes rapidly, particularly for dataset C (it's important to mention that ULM cannot recover any tracks with a DS greater than 6).This phenomenon underscores ULM's challenge in accurately reconstructing the morphology of vessels in a clinical setting marked by limited frame rate acquisition capabilities.
We investigated the efficacy of TEULM in reconstructing flow dynamics in terms of axial velocity and its magnitude.The maps of the axial velocity and velocity magnitude for the in vivo datasets are given in rows two and three of Figs. 3 (dataset A), 4 (dataset B), 5 (dataset C) and 6 (dataset D).The results confirm that recovering the flow dynamics is a more challenging task even when using TEULM, due to the loss of spatio-temporal information at lower frame rates.At DS rates greater than 4, we observe high deviations in the flow direction information even in larger vessels.TEULM's poor performance in recovering the velocity information is further highlighted by the R M S E evaluated on the velocity magnitude and flow direction maps for the in vivo datasets presented in Fig. 7.In particular, the error reaches above the 80% already at DS = 4 for all the datasets when using TEULM, and above 95% when using ULM.
Furthermore, TEULM misestimates the velocity magnitude, which tends to become more homogeneous as the DS rate increases for all the datasets, as seen in Figs. 3, 4, 5 and 6, for datasets A, B, C and D respectively.The loss of the magnitude information may be explained by the reduced amount of information available at lower frame rates, which yields to a loss of short micro-vessels, as well as a loss of longer tracks, which become shorter.To further analyze the velocity misestimation phenomena, we evaluate the mean velocity across all the tracks within their respective datasets in Fig. 9).As can be seen from the figure, there is a noticeable decrease in mean velocities for the down-Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.sampling rate of 10 for all the in vivo datasets.The results clearly indicate that accurately estimating velocities becomes a formidable challenge at lower frame rates, particularly when the microbubbles are moving at high speeds.Consequently, this leads to inaccuracies in velocity estimation.Moreover, the tracking algorithms heavily rely on both spatial and temporal information to match microbubbles across frames.In our current approach, the interpolation technique utilizes spatial and temporal information along a 2D plane to reconstruct the missing frames, see Fig. 1.Hence, we believe that implementing a 3D interpolation technique has the potential to enhance the accuracy of velocity estimation.This approach may better capture the intricate dynamics of microbubble movement and improve the overall performance of the tracking process.
Another consideration concerns the clinical utility of TEULM.The proposed TEULM represents a step towards transitioning ULM into clinical practice, addressing the current constraint of high frame rate requirements.However, despite the striking similarity in global information between the SR images reconstructed at DS = 10 using TEULM and those obtained through high frame rate ULM, we observe a substantial loss in the definition of vessels, particularly in the intricate microvascular structure.This loss may adversely affect both the quantity and quality of clinical information within the SR images.Furthermore, we note a marked decline in the retrieval of dynamic information at DS values exceeding 4, which imposes a limitation on TEULM's applicability in clinical scenarios.Therefore, if these SR images are to be used in a clinical context, we assert the need for further enhancements in the representation of vessels within the microvasculature and in the recovery of dynamic information at sub-100Hz frame rates before considering the extension of utilization with higher DS factors.
There are some limitations that are not explored in the study.First, in this paper, we only tested TEULM on 2D datasets.As future work, we envision extending the proposed approach to a 3D dataset, since 3D TEULM would represent an additional and fundamental step towards facilitating the clinical translation of this technology.
Second, the impact of motion on TEULM performance has not been explored in this paper.However, motion may lead to errors in reconstructing the signal path and target localization when using higher DS rates.Therefore, addressing motion effects through the implementation of motion compensation techniques could be explored as a potential avenue for future research.In addition, we aim at implementing a 3D interpolation technique.We anticipate that this approach will enhance the performance of TEULM in capturing dynamic profiles, even at lower frame rates.

V. CONCLUSION
This study presents TEULM, a technique to alleviate the high frame rate data acquisition (in the range of 1kHz), generally necessary to frame-to-frame pair the microbubbles in the ULM framework.Relaxing such constraint would promote the ULM clinical transition, as the US systems used in hospital can reach a frame rate of a magnitude lower.Results show that TEULM allows the generation of SR density maps with data acquired at frame rates as low as 50 Hz (an order of magnitude lower than the standard ULM).On the other hand, even if still capable to relax the frame rate by factor 4, TEULM required slightly higher frame rates (250 Hz) for the generation of consistent velocity maps with respect to the high frame rate's generated maps.the invertibility of the interpolation matrix for GA and MQ functions [38] since these two are the RBFs that are used in this work.
Definition 1: The following conditions should be satisfied to have a positive definite function φ on the domain R d .Firstly, the function φ must be even.Secondly, the following condition should be fulfilled for all the points X 1 , X 2 , . . ., X N and any vector v = [v 1 , v 2 , . . ., v N ] T : Remark: The function φ is strictly positive definite if and only if the null vector is the only vector v that converts the inequality in Eq. ( 11) into an equivalence.Since all of the eigenvalues of the interpolation matrix are positive with a strictly positive definite φ, is positive definite and thus invertible.
Remark 2: Despite the invertibility of matrix , a high condition number of this matrix may still render the interpolation problem ill-posed.As a consequence, we must choose the appropriate τ parameter in order reduce the error.The parameter τ can be set using a variety of methods, such as the Hardy [41] and the Frank [42] approaches.It should be noted that the τ value is typically chosen at random in the literature, based on the interpolation matrix condition number.

B. Error Analysis
The rate of convergence (RoC) of the RBF technique for adequately smooth functions is exponential (e.g., the RoC for MQ is O(ℑ N ), where 0 < ℑ < 1).However, for other methods, including the finite differences and polynomial spline bases, the RoC is in the order of O(N −m ) for a given constant m.
The MQ technique exponentially converges to proper smooth functions for a particular shape parameter value (τ ).In this case, h, i.e., the distance between X and X j , is reduced to its minimum upper bound: where indicates a series of centers and represents the interpolation space.Therefore, we can formulate the interpolation error as follows: In the other points (i.e., the middle points), the error is expressed as: The same findings are obtained for the GA basis in [43].

Fig. 2 .
Fig. 2. Reconstructed density maps (first row), flow direction map (second row) and velocity magnitude maps (third row) for the in vivo datasets using the original high frame rate data (1kHz for dataset A, B and D and 500Hz for dataset C).The colorbar of the density maps shows the number of counts of MBs centroids found in each pixel of the image.The colorbar of the Flow direction map indicates values in MBs/(pix.sec).The Velocity magnitude color-map encodes different mean velocities.Red corresponds to the maximum velocity (V max ) of 70 mm/s, 40 mm/s, 48 mm/s and 63 mm/s for dataset A, B, C and D, respectively).

Fig. 3 .
Fig. 3. ULM-and TEULM-generated density maps (first row), flow direction map (second row), and velocity magnitude map (third row) for Dataset A at different frame rates of 500Hz (DS = 2), 250Hz (DS = 4) and 100Hz (DS = 10).The color bar of the figure is the same as original ones in Fig. 2.

Fig. 6 .
Fig. 6.ULM-and TEULM-generated density maps (first row), flow direction map (second row), and velocity magnitude map (third row) for Dataset D at different frame rates of 500Hz (DS = 2), 250Hz (DS = 4) and 100Hz (DS = 10).The color bar of the figure is the same as original ones in Fig. 2.

Fig. 7 .
Fig. 7. Heatmaps representing the normalized Root Mean Square Error (RMSE) of the reconstructed density and velocity maps at DS = 2, 4, 8, 10 for TEULM and ULM for the in vivo datasets (A, B, C and D).The RMSE is normalized with respect to the maximum value for each dataset.

Fig. 8 .
Fig.8.Variation of the percentage of preserved tracks (shown as bar plots) and the percentage of saturation (shown as line plot) versus the changes in DS rate (for the in vivo datasets).The results are normalized with respect to the tracks and the saturation of the reconstructed images using the high frame rate datasets, corresponding to 1KHz for Dataset A, B and D, and 500Hz for Dataset C. For dataset C, since ULM fails to recover any track for DS >7, the saturation and the number of preserved tracks in equal to 0.

Fig. 9 .
Fig.9.Mean velocities (mm/s) evaluated for ULM at the original frame rate and TEULM at down-sampling rate 10, for all the in vivo datasets.

TABLE I RBF
BASIS USED IN THIS STUDY (ϵ IS A CONSTANT, r IS THE RADIUS OF THE RBF, AND Ln(.) INDICATES THE NATURAL LOGARITHM)

TABLE III ULM
AND TEULM PARAMETERS USED IN THIS WORK.* AND * * INDICATE THE PARAMETERS WHICH WERE ADJUSTED FOR THE NON-INTERPOLATED DATA (ULM AT DS > 1).IN PARTICULAR, * AND * * RESPECTIVELY MULTIPLY AND DIVIDE BY A FACTOR EQUAL TO THE DS RATE USED

TABLE IV SPATIAL
RESOLUTION (IN µM) MEASURED BY THE FRC CURVE AND THE 1/2 BIT THRESHOLDING METHOD FOR THE HIGH FRAME RATE ULM (1KHZ FOR DATASETS A, B AND D AND 500HZ FOR DATASET C) AND TEULM AT DIFFERENT DS RATES FOR THE in Vivo DATASETS

TABLE V DICE
SCORE FOR ALL THE in Vivo DATASETS OF ULM AND TEULM DENSITY MAPS.⋆ INDICATES THAT FOR TEULM AT DS=10, WE PRESENT THE MEAN AND THE STANDARD DEVIATION OF THE DICE SCORE WITH RESPECT TO TWO PERMUTED DATASETS