# BSRPCA: A Simplified Blind Super-Resolved RPCA-Based Approach for Enhancing Blood Flow Estimation

Xuan-Hieu Le, Jean-Pierre Remeniéras, Denis Kouamé, Duong-Hung Pham

► To cite this version:

Xuan-Hieu Le, Jean-Pierre Remeniéras, Denis Kouamé, Duong-Hung Pham. BSRPCA: A Simplified Blind Super-Resolved RPCA-Based Approach for Enhancing Blood Flow Estimation. 2025 IEEE International Conference on Image Processing (ICIP 2025), Sep 2025, Anchorage, AL, United States. hal-05008877

HAL Id: hal-05008877

https://hal.science/hal-05008877v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# BSRPCA: A SIMPLIFIED BLIND SUPER-RESOLVED RPCA-BASED APPROACH FOR ENHANCING BLOOD FLOW ESTIMATION

Xuan-Hieu Le*, Jean-Pierre Remeniéras†, Denis Kouamé*, and Duong-Hung Pham*

*Université de Toulouse, IRIT Laboratory, CNRS UMR 5505, Toulouse, France

†Université de Tours, INSERM, Imaging Brain & Neuropsychiatry (iBraiN), UMR 1253, Tours, France

## ABSTRACT

Ultrasound Localization Microscopy (ULM) offers great advances compared to conventional ultrasound imaging, enabling the reconstruction of microvascular structures with micrometer-scale precision. Most model-based localization approaches for ULM rely on a sequential five-step framework, which complicates parameter tuning and results in intricate tissue filtering. In addition, insufficient consideration of point spread function (PSF) also affects the final results. End-to-end data-driven solutions offer promising results yet require training data and ground-truths. In this study, we introduced a novel approach, based on computational super-resolution referred to as BSRPCA - Blind Super-Resolved Robust Principal Component Analysis, to address these challenges. Our method replaces the first three steps of the ULM process with a single super-resolution step. Experimental results on both *in silico* and *in vivo* datasets have proven the effectiveness of our method compared to the other benchmarks.

**Index Terms—** Ultrasound, deconvolution, super resolution, ultrasound localization microscopy.

## 1. INTRODUCTION

Ultrasound Localization Microscopy (ULM) is an emerging field in biomedical ultrasound that holds significant potential to overcome the diffraction limit in conventional ultrasound [1, 2]. By detecting, localizing, and tracking individual intravascular microbubbles (MBs) injected into the blood-stream, ULM enables the reconstruction of deep vascular structures with micrometer-scale spatial resolution, showcasing its unique capabilities in clinical applications in both animals and humans, across different organs such as the brain, breast, kidney or liver [3].

However, ULM technique also faces a trade-off between localization accuracy and acquisition time [1]. On the one hand, lower concentrations of injected MBs yield non-overlapping and spatially sparse signals, improving accuracy and enhancing image quality but resulting in slower

acquisition times. On the other hand, higher concentrations of MBs speed up the process but may lead to overlapping signals that make it difficult to isolate MB signals, and thus negatively affecting localization accuracy and image quality. Therefore, enhancing ULM image quality in terms of contrast and resolution, in scenarios with high MB concentrations remains a crucial challenge.

Several approaches have been proposed to tackle the issue mentioned above. Traditional model-based ULM approaches usually involve a sequential framework consisting of five steps: 1) tissue filtering and preprocessing, 2) MB detection, 3) MB localization, 4) MB tracking, and 5) image rendering. However, approaches that adopt such pipelines (e.g. radial symmetry (RS)-ULM method from [4]) suffer from accumulated errors from steps. In contrast, data-driven approaches with deep learning have shown promise in reducing processing time and providing precise end-to-end reconstructions [3, 2]. Nonetheless, these methods often require a substantial amount of realistic training data, which are not available in clinical practice. Furthermore, they can produce data shifts if applied without fine-tuning and tend to be less interpretable.

To circumvent these limitations, we propose a novel approach termed BSRPCA based on the robust principal component analysis (RPCA) framework. The contributions of this paper are thereby 2 folds:

* Firstly, we introduced a model-based approach that jointly handles the three distinct steps of RS-ULM within a single super-resolution algorithm. This eliminates the need for complex localization techniques and does not require any training data. Specifically, our approach can jointly perform tissue filtering, and blind super-resolution in a unified manner.

* Secondly, through comprehensive experimentation, we demonstrated that the proposed approach achieves comparable performance to, or even exceeding the benchmark, RS-ULM in both *in silico* and *in vivo* datasets.

The remainder of this paper is structured as follows. Section 2 provides a background on the classical model-based pipeline that is frequently utilized for this task, as well as details of the proposed methodology. Section 3 presents experimental results with their discussion, and Section 4 concludes the paper.

<footnote>This research was supported by the French National Research Agency (ANR) SONATINE project with grant number ANR-23-CE45-0002-01. The code of this paper is available at: gitlab.irit.fr/minds/theseeric/BSRPCA.</footnote>
## 2. METHODOLOGY

This section first provides the audience with the background information and problem formulation and then presents the details of the methodology. For reference, we provide in Table 1, a list of commonly used mathematical notations with brief explanations. Further explanations are given upon the first appearance of these symbols and others in the text.

<table><caption>**Table 1:** Summary of mathematical notations used.</caption><thead><tr><th>Notation</th><th>Notion</th></tr></thead><tbody><tr><td>z</td><td>the scalar</td></tr><tr><td>z̃</td><td>the vector</td></tr><tr><td>Z</td><td>the matrix</td></tr><tr><td>Z<sup>H</sup></td><td>the conjugate transpose of Z</td></tr><tr><td>I<sub>N</sub></td><td>the identity matrix with specified dimension N</td></tr><tr><td>Soft<sub>λ</sub>(w)</td><td>the soft thresholding operator of w given a threshold λ</td></tr><tr><td>F</td><td>the unitary 3D Fourier transform matrix</td></tr><tr><td>Λ</td><td>the diagonal matrix whose elements are the 3D Fourier transform of the zero-padded PSF defined as [(1<sub>dc</sub><sup>T</sup> ⊗ I<sub>N<sub>t</sub></sub>) ⊗ (1<sub>dr</sub><sup>T</sup> ⊗ I<sub>N<sub>z</sub>N<sub>x</sub></sub>)] Λ.</td></tr><tr><td>Λ<sub>1</sub><sup>T</sup><sub>w</sub></td><td>the row vector of ones of dimension of w.</td></tr><tr><td>SVD(w,r)</td><td>the partial decomposition of w subjective to the first r components</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1:** Traditional sequential pipeline for ULM processing [4]. We propose an algorithm aiming to replace the first three steps in the pipeline.</figcaption></figure>

### 2.1. Background: The Classical ULM Pipeline

As aforementioned, a traditional pipeline for ULM processing involves a sequence of 5 steps as visualized by Fig. 1.

Since the quality of ULM output is directly related to the effectiveness of filters and localization algorithms, several studies have been dedicated to improving these steps.

The conventional approaches to tissue filtering typically depend on Singular Value Decomposition (SVD) [5] or Robust Principal Component Analysis (RPCA) [6]. These methods apply the decomposition as a cluster filtering step, thereby offering simplicity and efficiency. However, they do not consider the blurring effects induced by point spread functions (PSFs), which can significantly affect the quality in terms of resolution and contrast of image reconstruction. Consequently, subsequent studies, such as DRPCA (Deconvolutive Robust Principal Component Analysis) incorporated deconvolution in the filtering step to account for the influence of the PSFs that characterize the employed imaging systems [7]. More advanced approaches, such as one proposed in [7], introduced blind deconvolution alternatives to tackle the challenge of requiring prior knowledge about the PSFs.

On the other hand, to improve localization efficiency, several localization methods have been introduced, such as interpolation with Lanczos, Spline or Cubic; Gaussian fitting; Weighted Average (WA) and Radial Symmetry (RS) [4]. Among these, Radial Symmetry has been proven as one of the most effective methods. One of the most prominent works, RS-ULM [4], utilized SVD for tissue filtering and utilized Radial Symmetry for the localization step. This approach has been shown to work effectively with *in vivo* data.

While the sequential approaches, of which RS-ULM is a representative, may show promise, they suffered from several limitations. First, due to their sequential nature, their final results are affected by errors accumulated from previous steps. Secondly, relying on several steps while utilizing such a complex localization algorithm, may complicate the overall tuning process and application of the pipeline. Moreover, we argue that using SVD in the filtering step of RS-ULM may be inaccurate, as it does not adequately account for the PSF's impact.

Reflecting on these facts, we propose BSRPCA (Blind Super-Resolved Robust Principal Component Analysis) to replace and jointly handle the first three steps in the classical ULM pipeline within a single algorithm, based on resolving a super-resolution problem. By doing so, we could alleviate excessive tuning efforts, eliminate the need for a complex localization algorithm, and reduce accumulated errors.

### 2.2. Problem formulation

We consider the initial ultrasound in-phase and quadrature (IQ) signal, denoted as $S$, which is modeled as the sum of three components: the tissue signal $T$, the super-resolved MB signal $x$ observed through $DHx$, where $D$, $H$ are respectively
the decimation operator and the system PSF, and the noise N, as expressed in the following equation:

$$S = T + DHx + N, \quad (1)$$

where $x \in \mathbb{R}^{N_z dz N_x dx \times N_t dt}$ represents the super-resolved target image to be estimated. The terms $dx$ and $dz$ indicate the decimation rates in the spatial directions $x$ and $z$, respectively, while $N_x$, $N_z$, and $N_t$ are the original dimensions of the observed corrupted data. In this work, we assume no super-resolution along the temporal dimension, and therefore $dt$ is always set to 1. We recall that $D \in \mathbb{R}^{N_t \times N_h}$ and $H \in \mathbb{R}^{N_h \times N_h}$ are respectively the decimation operator (e.g. downsampling) and the blurring PSF. $N_h$ and $N_t$ represent the super-resolved spatial size and the original spatial size, for example, one can consider $N_h = d \times N_t$ with $d = dx \times dz$. Finally, $N$ is assumed to be an additive noise, classically used for ultrasound radiofrequency images.

Building on the previous approach that utilized RPCA framework as outlined in [7, 8], we consider the tissues to be low-rank due to their inherent minimal variations over time, and the target $x$ as sparse, reflecting the sparse nature of the MBs. This enables effective separation of tissue from signals [9]. The main goal of this paper is to develop an optimization strategy to estimate the true clean, noise-free, tissue-free, and super-resolved image $x$ from the corrupted data $y$. In other words, our optimization problem is formulated as:

$$[\hat{x}, \hat{T}] = \arg \min_{x, T} \|S - DHx - T\|_F^2 + \lambda \|x\|_1 + \mu \|T\|_* \quad (2)$$

where $\lambda$ and $\mu$ are hyperparameters defining sparsity and low-rankness of $T$.

## 2.3. Proposed BSRPCA Method

To optimize (2), we use the ADMM framework and transformed the problem into the optimization of the summation of two function blocks. In this scenario, convergence is typically guaranteed, as demonstrated in [10]. We begin by introducing two auxiliary variables, $Z$ and $y$, and define some additional variables as follows:

$$M_1 = [DH \ I_{N_h}]; M_2 = [I_{N_h} \ 0]; M_3 = [0 \ I_{N_h}],$$

$$\tilde{x} = [x \ y]^T; \tilde{Z} = [Z \ T]^T,$$

$$f(\tilde{x}) = \frac{1}{2} \|S - M_1 \tilde{x}\|_F^2,$$

$$g(\tilde{Z}) = \lambda \|M_2 \tilde{Z}\|_1 + \rho \|M_3 \tilde{Z}\|_* \quad (3)$$

where $I_{N_h}$ denotes an identity matrix that has the same dimension as $DH$. With these definitions, our optimization problem (2) becomes:

$$\arg \min_{\tilde{x}, \tilde{Z}} f(\tilde{x}) + g(\tilde{Z}) \quad \text{such that} \quad \tilde{x} = \tilde{Z}, \quad (4)$$

This equation has the augmented Lagrangian expression:

$$L(\tilde{x}, \tilde{Z}, \gamma) = f(\tilde{x}) + g(\tilde{Z}) + \frac{\rho}{2} \left\| \tilde{x} - \tilde{Z} + \frac{\gamma}{\rho} \right\|_2^2 + \text{cst} \quad (5)$$

where $\text{cst}$ denotes a constant, $\rho$ is the positive convergent parameter and $\frac{\gamma}{\rho}$ is the scale dual variable. Using the ADMM framework [11, 12], the following 3 iterative steps are introduced to address the above equation:

$$\tilde{x}^{k+1} = \arg \min_{\tilde{x}} f(\tilde{x}) + \frac{\rho}{2} \left\| \tilde{x} - \tilde{Z}^{k+1} + \frac{\gamma^k}{\rho} \right\|_2^2 \quad (6)$$

$$\tilde{Z}^{k+1} = \arg \min_{\tilde{Z}} g(\tilde{Z}) + \frac{\rho}{2} \left\| \tilde{x}^k - \tilde{Z} + \frac{\gamma^k}{\rho} \right\|_2^2 \quad (7)$$

$$\gamma^{k+1} = \gamma^k + \rho (\tilde{x}^{k+1} - \tilde{Z}^{k+1}) \quad (8)$$

Using the definitions in (3), the resolution of each step is:

$$y^{k+1} = \arg \min_{y} \frac{1}{2} \|S - DHx^k - y\|_F^2 + \frac{\rho}{2} \left\| y - T^k + \frac{\gamma_2^k}{\rho} \right\|_2^2 \quad (9)$$

$$x^{k+1} = \arg \min_{x} \frac{1}{2} \|S - DHx - y^{k+1}\|_F^2 + \frac{\rho}{2} \left\| x - Z^k + \frac{\gamma_1^k}{\rho} \right\|_2^2 \quad (10)$$

$$T^{k+1} = \arg \min_{T} \mu \|T\|_* + \frac{\rho}{2} \left\| y^{k+1} - T + \frac{\gamma_2^k}{\rho} \right\|_2^2 \quad (11)$$

$$Z^{k+1} = \arg \min_{Z} \lambda \|Z\|_1 + \frac{\rho}{2} \left\| x^{k+1} - Z + \frac{\gamma_1^k}{\rho} \right\|_2^2 \quad (12)$$

$$\gamma_2^{k+1} = \gamma_2^k + \rho (y^{k+1} - T^{k+1}) \quad (13)$$

$$\gamma_1^{k+1} = \gamma_1^k + \rho (x^{k+1} - Z^{k+1}) \quad (14)$$

Unlike the other steps, which have closed-form solutions, solving (10) requires more effort, as it involves tackling a single-image super-resolution problem ($l_2 - l_2$). Specifically, this involves reconstructing a high-resolution image $x$ from its blurred, downsampled, and noisy version, $S - y^{k+1}$, using a prior/regularization term $\|x - (Z^k - \frac{\gamma_i^k}{\rho})\|_2$. Inspired by previous work [13], we propose handling the decimation and blurring operators simultaneously by leveraging their specific properties in the frequency domain [14, 15]. This approach leads to an efficient analytical solution for (10). As proposed in [13], we also defined the prior of the target as an interpolated approximation derived from the observation, e.g. by replacing the prior term with $\|x - D^H(S - T)\|_2$ for better convergence. Moreover, previous studies have highlighted that solving step (11) is computationally intensive due to the need for a full singular value decomposition (SVD) of a large matrix. To improve efficiency and accelerate the algorithm, we adopted techniques similar to those in [7], replacing the full SVD with a partial decomposition limited to the first $r_f$ components. The pseudo-code for the proposed BSRPCA algorithm is provided in Algorithm 1.

## 3. RESULTS AND DISCUSSIONS

In this section, we first introduce the datasets and metrics used. Following this, we present a qualitative and quantita-
<caption>**Algorithm 1:** Pseudo algorithm of BSRPCA</caption>

Initialize: **x**, **Z**, **y**, **T**, $\gamma_1$, $\gamma_2$.

While not convergent do:

Step 1: $\mathbf{y}^{k+1} = \frac{1}{1+\rho}(\mathbf{S} - \mathbf{D}\mathbf{H}\mathbf{x}^k + \rho\mathbf{T}^k - \gamma_2^k)$.

Step 2: $\mathbf{x}^{k+1} = \mathbf{k} - \mathbf{F}^H \mathbf{\Delta}^H (dr\mathbf{I}_{N_l} + \mathbf{\Delta}\mathbf{\Delta}^H)^{-1} \mathbf{\Delta}\mathbf{F}\mathbf{k}$,
where $\mathbf{k} = \mathbf{F}^H \mathbf{\Delta}^H \mathbf{F}\mathbf{D}^H (\mathbf{S} - \mathbf{y}^{k+1}) + \rho (\mathbf{Z}^k - \frac{1}{\rho}\gamma_1^k)$.

Step 3: $\mathbf{T}^{k+1} = \mathbf{U}\mathbf{S}\mathbf{o}\mathbf{f}^\mu \frac{\rho}{\rho} [\Delta]\mathbf{V}^T$,
with $\mathbf{U}\Delta\mathbf{V}^T = \mathbf{S}\mathbf{V}\mathbf{D}(\mathbf{y}^{k+1} + \frac{1}{\rho}\gamma_2^k, r_f)$.

Step 4: $\mathbf{Z}^{k+1} = \mathbf{S}\mathbf{o}\mathbf{f}^\lambda \frac{1}{\rho} [\mathbf{x}^{k+1} + \frac{1}{\rho}\gamma_1^k]$.

Step 5: $\gamma_2^{k+1} = \gamma_2^k + \rho (\mathbf{y}^{k+1} - \mathbf{T}^{k+1})$.

Step 6: $\gamma_1^{k+1} = \gamma_1^k + \rho (\mathbf{x}^{k+1} - \mathbf{Z}^{k+1})$.

tive comparison of the results between our proposed methods and a standard benchmark, RS-ULM on different datasets. All experiments were conducted on a computer equipped with an Intel(R) Xeon(R) W-1290 CPU and 128GB of RAM. For BSRPCA's tuning, $\lambda$ is set empirically, and $\rho$ is fixed.

## 3.1. Datasets and Metrics

To validate the contributions of our algorithm, we evaluate its effectiveness using 2 main datasets: the augmented *in silico* dataset, and the *in vivo* WKY rat brain dataset.

The *in silico* dataset, originally adopted from [4], consists of 20 blocks, each containing 1,000 frames. As the original version only contained MBs without tissue patterns and had low MB concentrations, which are unrealistic under practical conditions, we augmented the dataset to generate another version that better reflects reality. Specifically, we first increased the overlapping MB concentration in the dataset by summing groups of four blocks to allow more MBs in every frame. Each group consists of 5 blocks that are separated by equal periodic time intervals, for example, the first group includes blocks 1, 6, 11, and 16. Subsequently, we projected the MB locations into a super-resolved space at a scale of 10 times to create the target. Following this, we convolved the target with a Gaussian PSF and downsampled the results 10 times before the tissue and white Gaussian noise addition to generate the input observations. For benchmarking purposes, we reported the performance of methods utilizing several metrics adopted from [4]: Jaccard Accuracy (%), Lateral Error, Axial Error, Structure Similarity Index Measure (SSIM), and Normalized Root Mean Square Error (nMSE).

The *in vivo* dataset, corresponding to an angiography of a WKY 200 grams rat brain using Doppler contrast- enhanced US (with bolus injections at regular intervals), was acquired at UMR 1253, iBrain, Universite de Tours in France. Experiments on the WKY rat were conducted with approval from the regulatory authorities at the University of Tours. The dataset includes 40,000 frames across 80 blocks. For our ex-

periments, we excluded the bottom part of the images which contained irrelevant information. For PSF estimation, we applied the same assumptions made in [7], that is to consider the PSF to be spatially and temporally invariant. We utilized the method proposed in [7] with block size of 19x29 to estimate PSF from the temporal mean of all filtered frames. Since the ground-truths are not available for this *in vivo* dataset, we employed 2 image quality metrics commonly utilized in *in vivo* studies to quantitatively compare the methods: Contrast Ratio [16], Saturation [4] and Fourier Ring Correlation [17].

<figure><img src="image_3.png" /><figcaption>(a) The target rendering on the augmented *in silico* dataset.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) RS-ULM rendering on the augmented *in silico* dataset.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) BSRPCA rendering on the augmented *in silico* dataset.</figcaption></figure>

<caption>**Fig. 2:** Rendering results of both methods on the *in silico* dataset.</caption>

## 3.2. Results and Discussion

Fig. 2 visualizes the qualitative results of our proposed method and the RS-ULM, Table 2 compares both approaches on different metrics on the augmented silico dataset.
<figure><img src="image_6.png" /><figcaption>(a) RS-ULM rendering on the *in vivo* dataset.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) BSRPCA rendering on the *in vivo* dataset.</figcaption></figure>

**Fig. 3:** Rendering results of both methods on the *in vivo* dataset. A careful examination of the regions highlighted by the green rectangles reveals the distinctions between the methods.

<table><thead><tr><th></th><th>RS-ULM</th><th>BSRPCA</th></tr></thead><tbody><tr><td>Jaccard Accuracy (%) ↑</td><td>4.82</td><td><strong>20.40</strong></td></tr><tr><td>nMSE ↓</td><td>0.058</td><td><strong>0.056</strong></td></tr><tr><td>SSIM (%) ↑</td><td>0.894</td><td><strong>0.903</strong></td></tr><tr><td>Lateral Error ↓</td><td>0.22</td><td><strong>0.16</strong></td></tr><tr><td>Axial Error ↓</td><td><strong>0.16</strong></td><td><strong>0.16</strong></td></tr></tbody></table>

**Table 2:** Reported metrics of approaches on the *in silico* dataset. Bold number indicate better performance. Directed arrows indicate direction of better performance in metrics.

<table><thead><tr><th></th><th>RS-ULM</th><th>BSRPCA</th></tr></thead><tbody><tr><td>Contrast Ratio (dB) ↑</td><td>26.32</td><td><strong>30.55</strong></td></tr><tr><td>Saturation (%) ↑</td><td>16.69</td><td><strong>21.86</strong></td></tr><tr><td>Fourier Ring Correlation ↓</td><td>0.301</td><td><strong>0.289</strong></td></tr></tbody></table>

**Table 3:** Reported metrics of approaches on the *in vivo* dataset. Numbers in bold indicate better performance.

Inspection of Fig. 2 and Table 2, one may remark that the proposed BSRPCA produced more accurate results. Specifically, BSRPCA demonstrated better performance in reconstructing the horizontal curved structures, implying its enhancement in scenarios of high overlapping MB concentration. Moreover, RS-ULM exhibits unusual oscillation behaviors on the small curved tubes possibly due to its lack of a convolution mechanism. In contrast, BSRPCA tends to maintain the overall structures better without losing details, as evidenced by its better SSIM score and the maintenance of four straight lines in the corners. The numerical results presented in Table 2 align with these qualitative findings, showing that BSRPCA outperformed RS-ULM in all quantitative metrics.

On the *in vivo* dataset, we report the qualitative results for both methods in Fig. 3. Table 3 presents the quantitative results

measured from the *in vivo* dataset using both methods. A careful and deep glance into Fig. 3 and Table 3 reveals the superiority of BSRPCA in the reconstruction of microstructures which are significant degraded when using RS-ULM in the same areas. While both methods produced similar imaging contrasts, our proposed BSRPCA offered enhanced contrast for small structures. The quantitative results in Table 3 are consistent with the qualitative outcomes indicating that BSRPCA produced slightly better results compared to RS-ULM.

## 4. CONCLUSION

In this paper, we introduced a novel approach named BSRPCA that can jointly perform tissue filtering and blind super-resolution within a single algorithm. The proposed technique eliminates the needs for complicated localization algorithms while offering a simplified procedure. Both quantitative and qualitative results from simulated and *in vivo* datasets demonstrated the effectiveness of the proposed method compared to the classical RS-ULM pipeline. However, BSRPCA comes with increased computational costs as a trade-off. Future work will therefore focus on optimizing its execution time. Additionally, the study assumes that the PSF is spatiotemporally invariant; thus, exploring alternative methods for PSF estimation (e.g., [18]) will be an important direction for future research. Lastly, investigating techniques such as deep unfolding, as demonstrated in [19], for automatic hyperparameter tuning will be another focus in future work.

## 5. ACKNOWLEDGMENTS

The authors thank Coretin Alix and Dr. Ayache Bouakaz (INSERM iBraiN, Tours, France) for their help in acquiring the *in vivo* dataset used in this study.
## 6. REFERENCES

[1] Yirang Shin, Matthew R Lowerison, Yike Wang, Xi Chen, Qi You, Zhijie Dong, Mark A Anastasio, and Pengfei Song, “Context-aware deep learning enables high-efficacy localization of high concentration microbubbles for super-resolution ultrasound localization microscopy,” *Nat. Commun.*, vol. 15, no. 1, pp. 2932, Apr. 2024.

[2] Claudia Errico, Juliette Pierre, Sophie Pezet, Yann Desailly, Zsolt Lenkei, Olivier Couture, and Mickael Tanter, “Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging,” *Nature*, vol. 527, no. 7579, pp. 499–502, Nov. 2015.

[3] Brice Rauby, Paul Xing, Maxime Gasse, and Jean Provost, “Deep learning in ultrasound localization microscopy: Applications and perspectives,” *IEEE Trans. Ultrason. Ferroelectr. Freq. Control.*, vol. 71, no. 12: Breaking the Resolution, pp. 1765–1784, Dec. 2024.

[4] Baptiste Heiles, Arthur Chavignon, Vincent Hingot, Pauline Lopez, Eliott Teston, and Olivier Couture, “Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy,” *Nat. Biomed. Eng.*, vol. 6, no. 5, pp. 605–616, May 2022.

[5] Charlie Demené, Thomas Deffieux, Mathieu Pernot, Bruno-Félix Osmanski, Valérie Biran, Jean-Luc Gennisson, Lim-Anna Sieu, Antoine Bergel, Stéphanie Franchui, Jean-Michel Correas, Ivan Cohen, Olivier Baud, and Mickael Tanter, “Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity,” *IEEE Trans. Med. Imaging*, vol. 34, no. 11, pp. 2271–2285, Nov. 2015.

[6] Ricardo Antonio Maronna, “Robust m-estimators of multivariate location and scatter,” *The Annals of Statistics*, vol. 4, no. 1, pp. 51–67, 1976.

[7] Duong-Hung Pham, Adrian Basarab, Jean-Pierre Remenieras, Paul Rodriguez, and Denis Kouame, “Fast high resolution blood flow estimation and clutter rejection via an alternating optimization problem,” in *2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI)*. Apr. 2021, IEEE.

[8] Kailiang Xu, Xingyi Guo, Yihui Sui, Vincent Hingot, Olivier Couture, Dean Ta, and Weiqi Wang, “Robust pca-based clutter filtering method for super-resolution ultrasound localization microscopy,” in *2021 IEEE International Ultrasonics Symposium (IUS)*. Sept. 2021, p. 1–4, IEEE.

[9] Yann Desailly, Anne-Marie Tissier, Jean-Michel Correas, Frédéric Wintzenrieth, Mickaël Tanter, and Olivier Couture, “Contrast enhanced ultrasound by real-time spatiotemporal filtering of ultrafast images,” *Physics in Medicine and Biology*, vol. 62, no. 1, pp. 31–42, Dec. 2016.

[10] Tianyi Lin, Shiqian Ma, and Shuzhong Zhang, “On the global linear convergence of the ADMM with Multi-Block variables,” *SIAM J. Optim.*, vol. 25, no. 3, pp. 1478–1497, Jan. 2015.

[11] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” *Found. Trends Mach. Learn.*, vol. 3, no. 1, pp. 1–122, 2011.

[12] Tianyi Lin, Shiqian Ma, and Shuzhong Zhang, “On the global linear convergence of the admm with multiblock variables,” *SIAM Journal on Optimization*, vol. 25, no. 3, pp. 1478–1497, 2015.

[13] Ningning Zhao, Qi Wei, Adrian Basarab, Nicolas Dobiogeon, Denis Kouame, and Jean-Yves Tourneret, “Fast single image super-resolution using a new analytical solution for 12-12 problems,” *IEEE Trans. Image Process.*, vol. 25, no. 8, pp. 3683–3697, Aug. 2016.

[14] Per Christian Hansen, James G. Nagy, and Dianne P. O Leary, *Deblurring Images*, SIAM J. Appl. Math, Jan. 2006.

[15] William W. Hager, “Updating the inverse of a matrix,” *SIAM Review*, vol. 31, no. 2, pp. 221–239, June 1989.

[16] Alfonso Rodriguez-Molares, Ole Marius Hoel Rindal, Jan D'hooge, Svein-Erik Måsøy, Andreas Austeng, Muyinatu A. Lediju Bell, and Hans Torp, “The generalized contrast-to-noise ratio: A formal definition for lesion detectability,” *IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, vol. 67, no. 4, pp. 745–759, 2020.

[17] V. Hingot, A. Chavignon, B. Heiles, and O. Couture, “Measuring image resolution in ultrasound localization microscopy,” *IEEE Transactions on Medical Imaging*, vol. 40, no. 12, pp. 3812–3819, 2021.

[18] Loïc Denis, Eric Thiébaut, Ferréol Soulez, Jean-Marie Becker, and Rahul Mourya, “Fast approximations of shift-variant blur,” *Int. J. Comput. Vis.*, vol. 115, no. 3, pp. 253–278, Apr. 2015.

[19] C Bertocchi, E Chouzenoux, M-C Corbineau, J-C Pesquet, and M Prato, “Deep unfolding of a proximal interior point method for image restoration,” *Inverse Problems*, vol. 36, no. 3, pp. 034005, Feb. 2020.