# Ultrasound matrix imaging for 3D transcranial in vivo localization microscopy

Flavien Bureau, Louise Denis, Antoine Coudert, Mathias Fink, Olivier Couture, Alexandre Aubry

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

## To cite this version:

Flavien Bureau, Louise Denis, Antoine Coudert, Mathias Fink, Olivier Couture, et al.. Ultrasound matrix imaging for 3D transcranial in vivo localization microscopy. Science Advances, 2025, 11 (11), pp.eadt9778. 10.1126/sciadv.adt9778. hal-04744243v2

HAL Id: hal-04744243

https://hal.science/hal-04744243v2

Submitted on 3 Oct 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_2.png" /></figure>
Ultrasound matrix imaging for 3D transcranial in-vivo
localization microscopy

Flavien Bureau<fnref n="1" />,†, Louise Denis<fnref n="2" />,†, Antoine Coudert<fnref n="2" />, Mathias Fink<fnref n="1" />,
Olivier Couture<fnref n="2" />, and Alexandre Aubry<fnref n="1" />,*

**Summary Sentence:** Ultrasound matrix localization microscopy provides a three-dimensional, artifacts-free and super-resolved image of brain vessels.

Transcranial ultrasound imaging is usually limited by skull-induced attenua-
tion and high-order aberrations. By using contrast agents such as microbub-
bles in combination with ultrafast imaging, not only can the signal-to-noise
ratio be improved, but super-resolution images down to the micrometer scale
of the brain vessels can be obtained. However, ultrasound localization mi-
croscopy (ULM) remains impacted by wave-front distortions that limit the
microbubble detection rate and hamper their localization. In this work, we
show how ultrasound matrix imaging (UMI), which relies on the prior record-
ing of the reflection matrix, can provide a solution to those fundamental issues.

As an experimental proof-of-concept, an in-vivo reconstruction of deep brain microvessels is performed on three anesthetized sheeps. The compensation of wave distortions is shown to drastically enhance the contrast and resolution of ULM. This experimental study thus opens up promising perspectives for a

<footnote n="1">¹ PSL University, ESPCI Paris, CNRS, Institut Langevin, Paris, France</footnote>

<footnote n="2">²Sorbonne Université, CNRS, INSERM, Laboratoire d'Imagerie Biomedicale, Paris, France</footnote>

<footnote>†These authors contributed equally to this work.</footnote>

<footnote>*To whom correspondence should be addressed; E-mail: alexandre.aubry@espci.fr.</footnote>
transcranial and non-ionizing observation of human cerebral microvascular pathologies, such as stroke.

INTRODUCTION

Transcranial imaging is essential to understand the complex vascular mechanisms underlying pathologies. For cerebrovascular events, such as strokes, X-rays and computed tomography (CT) are currently used to observe hemorrhagic events (?) and magnetic resonance imaging (MRI) to diagnose and date ischemic stroke (?). Although these methods allow visualization of many neurovascular diseases, MRI or CT-scans require large and expensive equipment that is only primarily available in health care centers.

Alternatively, transcranial ultrasound is used in Doppler mode to monitor reperfusion after a stroke (?). This technique allows portable real-time imaging at low cost. However, it remains a major challenge as the complexity of the skull layer with its unpredictable porosity causes severe attenuation, strong aberrations and multiple scattering, leading to a drastic reduction in resolution and contrast of the ultrasound image.

Microbubbles (MB) are often used to improve the visualization of blood flow in ultrasound imaging (?, ?, ?). In the clinic, they are injected intravenously and are excellent contrast agents due to their high impedance mismatch compared to soft tissues. The use of ultrafast ultrasound scanners, i.e. with high frame rates, has also enabled organ perfusion dynamics to be studied with greater precision (?, ?).

In addition to enhancing the blood signal, microbubbles can also be used to generate super-resolved images. This is the principle of ultrasound localization microscopy (ULM) (?, ?, ?, ?), which is based on a similar idea developed a few years earlier in fluorescence microscopy (?, ?). By separating the echoes of individual injectable microbubbles, their position can be localized with a sub-wavelength precision and tracked over time. This allows an accuracy that no longer
depends on the diffraction limit but on the accuracy of detection of their center. As a consequence, ULM results in a significant improvement in resolution by a factor of ten (?) compared to standard ultrasound images. ULM has already proved its worth in many organs in depth, both pre-clinically (?, ?, ?, ?, ?, ?, ?, ?) and clinically (?, ?, ?, ?, ?, ?, ?, ?). The real-time information contained in the displacement of microbubbles enables complex vascular networks to be mapped at the micrometer scale, even going as far as the capillary-level reconstruction of functional units, such as glomeruli in the kidney (sULM) (? , ?, ?). ULM imaging has even been successfully used in transcranial clinical imaging of the adult human brain (?), where aberration correction is critical (? , ?).

Recently, 3D ULM ( ?, ?, ?, ?, ?, ?, ?, ?, ?) has overcome many 2D limitations, such as the
projection of 3D structures in 2D, suboptimal tracking of microbubbles, and out-of-plane probe
movements that prevent vascular reconstruction. Nevertheless, aberrations caused by the skull
remain a problem for ULM, especially in 3D. The sensitivity of bubble detection is lower, which
affects the contrast of the super-resolved image, but also the accuracy of bubble localization,
where the high sidelobes of a degraded PSF can cause artifacts such as vessel doubling ( ?, ?).

To solve the fundamental problem of aberrations, adaptive focusing was introduced years
ago in astronomy to compensate for phase distortions undergone by starlight when it passes
through the atmosphere ( ?, ?, ?, ?). This idea was then transposed to ultrasound imaging by
tailoring a focusing law that either optimizes an image metrics ( ?, ?) or maximizes the spatial
correlations of back-scattered echoes ( ?, ?, ?, ?). In the specific case of ULM, several methods
have been proposed to correct the aberrations caused by the skull. A first strategy (?) is based
on deep learning algorithms trained with numerical simulations. Although it gives impressive
in-vivo results for small animals, it is hampered by the specificity of the training set, so its
transfer to other experimental configurations (probe, organ, SNR) may require tedious training
procedures. The second method (?) uses microbubbles as ultrasound guide stars but suffers
several drawbacks. On the one hand, microbubbles should be detected and isolated enough, which can be extremely difficult to assess because of the strongly distorted focal spot across the skull. On the other hand, the number of detected micro-bubbles should be large enough for self-averaging in each isoplanatic patch. Last but not least, the demonstration was limited to 2D imaging using a 1D phased array probe, while skull heterogeneities are 3D distributed and require a bi-dimensional control of the ultrasonic wave-field to properly compensate for such high-order aberrations (?).

Interestingly, inspired by previous seminal works (?, ?, ?), a matrix approach of ultrasound imaging has been recently developed to overcome aberrations without any guide star (? , ?). Based on the recording of the reflection associated with the ultrasonic probe and the linearity of the wave equation, it allows the simulation of wave focusing in post-processing and an iteration of this process towards local adaptive focusing laws in ultrasonic speckle (?). An optimized contrast and close-to-ideal resolution is then obtained for each voxel of the ultrasound image. Ultrasound Matrix Imaging (UMI) has been successfully applied to transcranial 3D imaging of a head phantom (?).

The goal of this paper is to demonstrate the interest of UMI for in-vivo 3D brain imaging.
In particular, we will show how: (i) UMI can be used to quantify aberrations and multiple
scattering in in-vivo transcranial imaging; (ii) UMI can exploit static speckle from brain tissue to
tailor complex 2D focusing laws in a strong aberration and multiple scattering regime; (iii) ULM
can leverage these focusing laws despite motion to increase the detection rate of microbubbles
and improve their localization; (iv) UMI can be fruitfully combined to ULM in order to provide
a three-dimensional, artifacts-free, super-resolved and contrasted image of brain vessels.

To that aim, an experimental proof-of-concept will consist in brain imaging on three anes-
thetized sheeps. The first step of the experimental procedure consists in the acquisition of the
reflection matrix and the estimation of local focusing laws by means of UMI. The second step
consists in a microbubble injection and the beamforming of ultra-fast images using the focusing
laws provided by UMI. The last part consists in the usual ULM post-processing that consists
in localizing and tracking micro-bubbles. Final ULM image contrast will be shown to be dras-
tically improved by the better detection rate of smaller micrometer-sized vessels provided by
UMI. The comparison with MRI will show how UMI allows the removal of artifacts from which
the initial ULM image suffered.

RESULTS

<figure><img src="image_3.png" /><figcaption>Figure 1: **Sheep experiment : A multi-sequence acquisition.** (A) Angiographic MRI. (B) Transcranial power Doppler for positioning the probe. (C) Acquisition of the reflection matrix to estimate local aberration laws with 3D-UMI. (D) Ultrafast imaging for super-resolved ULM images. The images shown in this figure are for illustrative purposes only and do not provide quantitative results.</figcaption></figure>

**Reflection matrix in a virtual source basis**

In this paper, we demonstrate transcranial brain imaging by overcoming the major challenges
of: (i) skull bone aberration with 3D UMI; (ii) diffraction-limited resolution with the super-
resolution feature provided by ULM. As a proof-of-concept, the same experimental protocol
(Fig. ??) has been applied to three anesthetized sheeps (Methods). A 32 × 32 array of trans-
ducers is placed on the shaved head of each sheep directly over the crest of the frontal bone (Fig. ??B). The probe is driven in 1-2 frequency bandwidth in order to limit the attenuation by skull bone (?). The first ultrasound sequence consists in recording a high-dimension reflection matrix $R_{uv}(\tau)$ (Fig. ??C) that contains the scattered wave-field $R(u_{out}, v_{in}, \tau)$ recorded by each transducer $u_{out}$ of the probe as a function of the echo time $\tau$ for a set of 324 incident diverging waves. This illumination basis is associated with a set of virtual transducers $v_{in}$ placed behind the probe (see Methods and Fig. ??). It maximizes the signal-to-noise ratio across the targeted field-of-view (?). The number of incident waves is dictated by our choice of the maximal angular range spanning from -30° to +30° (Methods).

<figure><img src="image_4.png" /><figcaption>Figure 2: **Principle of ultrasound matrix imaging.** (A) Focused reflection matrix contains the time-gated response between virtual source ($\rho_{in}$) and detector ($\rho_{out}$) located at the same depth $z$. (B) An output projection of the focused reflection matrix in the transducer basis ($u_{out}$) provides the reflected wave-fronts induced by each virtual source ($\rho_{in}$) at depth $z$. (C) Those wave-fronts are realigned to extract the wave-front distortions induced by the mismatch between the real speed-of-sound distribution and the wave velocity model. Seen from the focused basis, this operation leads to an angular de-scan of each virtual source at the same position, leading to the synthesis of a virtual guide star. (D) Exploiting the correlations between each distorted wave-front, an iterative phase reversal algorithm extract an aberration phase law that can be used to compensation for wave distortions induce by the skull and ideally retrieve a diffraction-limited focal spot across the field-of-view.</figcaption></figure>
Quantifying skull-induced aberrations and scattering

In contrast with standard beamforming that relies on the principle of confocal imaging, the
cornerstone of UMI consists in decoupling the input and output focal spots in post-processing.
Mathematically, such a beamforming procedure can be written as follows:

$$
R(\rho_{\text{out}}, \rho_{\text{in}}, z) = \sum_{\mathbf{v}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} R(\mathbf{u}_{\text{out}}, \mathbf{v}_{\text{in}}, \tau_{\text{in}}(\mathbf{v}_{\text{in}}, \rho_{\text{in}}, z) + \tau_{\text{out}}(\mathbf{u}_{\text{out}}, \rho_{\text{out}}, z)), \quad (1)
$$

where **r**<sub>in</sub> = (ρ<sub>in</sub>, z) and **r**<sub>out</sub> = (ρ<sub>out</sub>, z) are the coordinates of the input and output focal points,
respectively. τ<sub>in</sub> and τ<sub>out</sub> are the expected travel times for the wave to go from the source/detector
to the input/output focusing point (Methods, Eq. ??). The result of Eq. ?? is a focused reflection
matrix **R**<sub>ρρ</sub>(z) whose coefficients R(ρ<sub>out</sub>, ρ<sub>in</sub>, z) probe the cross-talk between a virtual source
and detector at lateral positions ρ<sub>in</sub> and ρ<sub>out</sub>, respectively, both located at the same depth z
inside the medium (Fig. ??A). This matrix enables a local quantification of the focusing quality
and multiple scattering by providing a local point spread function in reflection (RPSF) (?, ?, ?)
(Methods, Eq. ??). Figure ??B display maps of RPSFs at depth z = 35 mm for each sheep. In
the single scattering regime and in absence of aberrations, the RPSF would exhibit a confocal
diffraction-limited peak. On the contrary, the various RPSFs here display a distorted focal spot
induced by wave velocity heterogeneities on top of an incoherent background due to multiple
scattering events induced by the skull. Strikingly, whatever the imaging depth, the RPSF looks
more distorted and less contrasted for sheep n°6 than for sheep n°4 and 5 (Fig. ??B).

On the one hand, the aberration level can be quantified by investigating the spatial extent
δρ of the RPSF (?, ?) (Fig. ??). The depth evolution of δρ is shown in Fig. ??C. A general
trend is a degradation of the resolution with depth as predicted by diffraction theory but also a
more important aberration level in sheep n°6. This difference of behavior can be understood by
looking at the micro-CT of each skull (Fig. ??A). While the skulls of sheep n°4 and 5 show a
regular thickness ($L \sim 4.2 \pm 1.1$ mm and $L \sim 5.3 \pm 1$ mm, respectively), skull n°6 displays a more irregular shape: $L \sim 5.7 \pm 1.7$ mm (?). This difference of morphology between the three skulls could explain the higher degree of aberrations observed in sheep n°6.

On the other hand, the multiple scattering rate can be estimated from the incoherent background of the RPSF (Methods, Eq. ??). It appears that the multiple scattering rate is far from being negligible even after the beamforming process since it reaches a value of 30 % in the $z = 35$ mm-region for sheep n°6 (Fig. ??D). As the aberration level, the multiple scattering rate is larger for this sheep over almost the whole depth range. This feature can be understood by looking at the internal structure of the bone. Indeed, the diploë volume ratio $v$ is much higher in sheep n°6 ($v = 57\%$) than in sheeps n°4 and 5 ($v = 14\%$ and $v = 29\%$, respectively) (?). This high degree of heterogeneity seems to account for the higher multiple scattering rate observed in sheep n°6.

## Overcoming wave distortions

The complex speed-of-sound distribution highlighted by Fig. ??B hampers ultrasound brain imaging. The confocal image of the brain, that can be extracted from the diagonal coefficients of $R_{\rho\rho}(z)$ (Methods, Eq. ??), is blurred by the skull heterogeneities (Fig. ??A). Fortunately, an optimized contrast and diffraction-limited resolution can be recovered using the distortion matrix concept (?, ?, ?). This matrix contains the wave distortions exhibited by each reflected wave-front in the transducer basis (Fig. ??C). The choice of this basis is dictated by the fact that most aberrations are induced by the skull heterogeneities, i.e in the vicinity of the probe. An iterative phase reversal (IPR) process is then applied to the distortion matrix to provide an estimation of an aberration phase law, $\phi(\mathbf{u}_{\text{out}}, z)$, at each depth (Methods, Fig. ??D and Movie S1). Examples of aberration phase laws are displayed in Fig. ??C. They show a complex feature,
<figure><img src="image_5.png" /><figcaption>Figure 3: **Influence of cranial heterogeneities on aberration and multiple scattering.** (A) Micro-CT of each sheep skull. Left: sheep n°4, middle: sheep n°5, right: sheep 6. (B) Corresponding RPSF maps at depth z = 35 mm. (C and D) RPSF extension and multiple scattering rate as a function of depth, extracted by UMI, for sheeps n°4 (green), 5 (purple line) and 6 (orange line).</figcaption></figure>

as quantified by the Strehl ratio *S* that provides the ratio between the energy at focus with and without aberrations (?) (Methods, Eq. ??). For instance, the aberration law measured at depth *z* = 50 mm exhibits a Strehl ratio *S* = 0.03. This extremely low value illustrates the detrimental impact of skull heterogeneities on the focusing process. A more detailed quantification of the aberration magnitude is provided in Section ?? by evaluating the modulus of the aberration transfer function. As expected from the previous RPSF analysis (Fig. ??), the aberration magnitude is found to be 30 % larger in sheep n°6 compared to sheep n°4 and 5.
The aberration phase laws can be leveraged to compensate for the wave distortions exhibited by the focused reflection matrix (Methods). The resulting confocal image is displayed in Fig. ??E. It exhibits a clear improvement compared to the initial confocal image (Fig. ??A), with a contrast gain of 8 dB for the deepest bright scatterers ($z = 50$ and $60$ mm). The efficiency of the correction process can be assessed by comparing the maps of RPSFs before (Fig. ??B) and after (Fig. ??D) aberration compensation. The spatial resolution is improved by a factor between 2 and 3 and the multiple scattering background is decreased by more than 10 dB beyond $z = 50$ mm (Fig. ??). Nevertheless, the obtained RPSFs remain imperfect in Fig. ??D. High-order aberrations subsist because of their anisoplanicity. Last but not least, even though the confocal image shows a clear improvement, the diffraction-limited resolution $\delta\rho_0$ remains limited by the probe size $D$: $\delta\rho_0 \sim \lambda z/D$.

## Matrix imaging for ultrasound localization microscopy

To overcome diffraction, ULM is then performed by injecting boluses of microbubbles and tracking them using an ultra-fast imaging sequence (Methods). To that aim, a hybrid transmit basis is used with three cylindrical waves emitted by the entire probe and two spherical diverging waves transmitted successively by each probe panel (?). This sparse emission basis, referred to as $s_{in}$, yields a volume rate of 209 Hz (Methods). For each incident wave-field, the raw data is stored in a set of new reflection matrices $R'_{us}(\tau, t_m) = [R'(u_{out}, s_{in}, \tau, t_m)]$ recorded at different times $t_m$. A set of three-dimensional images is then obtained by beamforming the reflection matrix rephased by the conjugate of the aberration law to compensate for skull-induced
<figure><img src="image_6.png" /><figcaption>Figure 4: Ultrasound matrix imaging in the sheep brain. (A) Confocal volume extracted from the diagonal of the focused R-matrices (Methods, Eq. ??). (B, C, and D) Maps of initial RSPFs, aberration phase laws, and corrected RPSFs, respectively, at three different depths $z = 32$ mm (top), $z = 50$ mm (middle) and $z = 68$ mm (bottom). Each RPSF is displayed in a scan range $\Delta\rho = \rho_{out} - \rho_{in}$ (see Methods) that varies from -15 to +15 mm in both x and y directions. (E) Confocal image after aberration correction. The results shown here correspond to ultrasound data acquired on sheep n°6.</figcaption></figure>

aberrations:

$$ \mathcal{I}'(\mathbf{r}, t_m) = \sum_{\mathbf{s}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} \exp[-j\phi(\mathbf{u}_{\text{out}}, z)] \times R'(\mathbf{u}_{\text{out}}, \mathbf{s}_{\text{in}}, \tau_{\text{in}}(\mathbf{s}_{\text{in}}, \rho, z) + \tau_{\text{out}}(\mathbf{u}_{\text{out}}, \rho, z), t_m) \quad (2) $$

After filtering out the static component in the image series (Methods), the dynamic component of the brain can be highlighted. A two-dimensional cross-section of the resulting image is displayed in Fig. ??B at a given time $t_m$. This image is compared to its raw counterpart (Fig. ??A) obtained without aberration compensation. Thanks to UMI, isolated microbubbles show less distorted focal spots, as shown by two examples in the purple and red squares in Figs. ??A,B. Just above (~50 mm), a halo of microbubbles in the blue square seem to indicate the pres-
ence of a vessel but the individual detection of each microbubble is complicated by their high concentration.

<figure><img src="image_7.png" /><figcaption>Figure 5: **Enhancing micro-bubble localization accuracy with UMI.** (A and B) Sagittal middle slice of the ultrasound image without and with aberration correction at a given time, respectively. The position of the detected microbubbles is highlighted by red dots. The blue and orange boxes show areas where false alarms are reduced by reducing the side lobes of the PSF. The purple region shows a region where aberration compensation allows the detection of a microbubble. The red square highlights a zone where a microbubble is reassigned at a new location. (C) Number of localized microbubbles as a function of the acquisition time *t* with (green) and without (pink) UMI. (D) Increase of the microbubble detection rate with UMI for different track lengths. The data shown here correspond to sheep n° 6 (acquisition 4).</figcaption></figure>

A deconvolution of the dynamic image with a theoretical PSF and the application of a detection threshold provides a localization of bubbles in the field-of-view (Methods). The position of the detected bubbles are highlighted with red points in Fig. ??A and B. While several microbubbles are wrongly assigned to the high side-lobes of each distorted focal spot in Fig. ??A, UMI drastically reduces this number of false detection in Fig. ??B by restoring a satisfying PSF, in particular in the colored squares. Despite this reduction of false alarms, the number of localized
microbubbles is increased. Figure ??C illustrates this assertion by showing, for sheep n°6, the higher number of detected bubbles over time with UMI. The different peaks corresponds to each injection of microbubbles.

To further reduce the false alarm rate and better appreciate the effect of UMI, microbubbles can be tracked over time (Methods). In absence of aberration correction, the focal spot is distorted and bubbles cannot be detected on each frame during their journey across the field-of-view. The bubble trajectory is then subdivided into shorter tracks. As a consequence, the track length is a relevant observable to demonstrate the benefit of the aberration correction process (?). Fig. ??D shows the increase of detected tracks provided by UMI as a function of their length. UMI significantly increases the number of detected tracks by a factor ranging from 40 to 250% according to the track length. Table ?? shows the same observable for different minimum path lengths. The general trend is the same for each acquisition, thereby showing the robustness of UMI and its benefit in each imaging configuration. This result is crucial for ULM images as the longest tracks are the most reliable to delineate small vessels.

<table><caption>Table 1: **Quantification of UMI benefits on ULM tracks.** Percentage of additional tracks provided by UMI for various minimal tracks length expressed in frames.</caption><thead><tr><th rowspan="2">Sheep number</th><th rowspan="2">Acq. number</th><th colspan="5">Additional tracks detected with UMI correction [%]</th></tr><tr><th>≥ 2 frames</th><th>>3 frames</th><th>>5 frames</th><th>>10 frames</th><th>>30 frames</th></tr></thead><tbody><tr><td>S4</td><td>1</td><td>24.0</td><td>41.7</td><td>49.1</td><td>63.3</td><td>121.7</td></tr><tr><td rowspan="2">S5</td><td>1</td><td>14.1</td><td>26.6</td><td>35.7</td><td>54.5</td><td>76.6</td></tr><tr><td>2</td><td>7.5</td><td>13.7</td><td>16.2</td><td>28.5</td><td>42.8</td></tr><tr><td rowspan="5">S6</td><td>1</td><td>15.4</td><td>16.3</td><td>19.2</td><td>23.5</td><td>43.5</td></tr><tr><td>2</td><td>8.7</td><td>13.3</td><td>15.3</td><td>19.1</td><td>36.1</td></tr><tr><td>3</td><td>6.7</td><td>11.5</td><td>13.2</td><td>17.3</td><td>18.3</td></tr><tr><td>4</td><td>6.9</td><td>22.3</td><td>32.5</td><td>47.2</td><td>66.2</td></tr><tr><td>5</td><td>8.7</td><td>27.9</td><td>32.5</td><td>53.5</td><td>77.1</td></tr><tr><td colspan="2">Mean</td><td>11.5</td><td>21.7</td><td>26.7</td><td>38.4</td><td>60.3</td></tr><tr><td colspan="2">Std</td><td>6.0</td><td>9.0</td><td>11.2</td><td>16.1</td><td>28.5</td></tr></tbody></table>

ULM density maps are built by counting the number of tracks crossing each voxel of the field of view. Such volumes are shown from two different angles in Fig. ???. The typical structure
of the cerebral vessels can already be recognized before aberration correction (Figs. ??A and C). This structure is known as the Willis circle (?) and can be seen as the main cross-road of the main vessels in the brain.

<figure><img src="image_8.png" /><figcaption>Figure 6: Correcting ULM images with UMI. 3D images shown as maximum intensity projections from two different points of view. ULM volumes (A, C) before and (B, D) after aberration compensation. The polygon of Willis is depicted by a dashed white line and is located around $z = 40 – 45$ mm. The same dynamic range is displayed for all images. The data shown here correspond to sheep n° 6, acquisition 4.</figcaption></figure>

However, the skull heterogeneities also induce artifacts such as the duplication of one artery that can be noticed at the bottom of the conventional ULM image (Figs. ??A and C). This type of artifact has already been observed in transcranial imaging in both mouse (?) and human (?). Strikingly, the focusing laws extracted by UMI (Fig. ??C) allow the removal of such artifacts as highlighted by the corrected ULM images in Figs. ??B and D. More generally, the higher detection rate and correct repositioning of micro-bubbles lead to a better contrast for ULM as highlighted by the vessel network on top of the Willis circle ($z = 25–35$ mm) that appears much brighter (+6 dB) and better resolved. A similar resolution and contrast improvement can be observed for the other 3D ULM images shown in Fig. ?? and corresponding to other
acquisitions on the three sheeps.

To appreciate more quantitatively the benefit of UMI, the ULM density maps are superim-
posed on MRA, which is considered as the gold standard for observing large vascular structures
in the brain. After UMI, the ULM density maps match more closely with the vessels revealed by
MRA. In the coronal slice (Fig. ??A,C), UMI allows a more complete and brighter reconstruc-
tion of the vascular network (Fig. ??C) that matches in depth with the Willis polygon revealed
by MRA, especially in the anterior cerebral arteries that deviate from the polygon (green and
blue arrows). A zoom on the right anterior cerebral artery is displayed in Fig. ?? and highlights
the drastic gain in contrast provided by ULM compared to MRA around the Willis circle. This
benefit of UMI can be quantified by a loss of image entropy (?) increased by a factor 5 with
UMI.

In the sagittal section (Fig. ??B,D), ULM provides a better adequation between MRA and
ULM images of venous system (yellow arrows). This observation can be quantified by a com-
puting structural similarity index between MRA and ULM images. This similarity metrics is
increased by a factor 2 with UMI (Fig. ??). Besides the repositioning of large vessels, the in-
crease in the total SNR of the ULM volume and the appearance of small vessels provided by
UMI are noteworthy compared to the initial ULM image but also with respect to MRA that only
provides an image of the largest vessels (Fig. ??).

Finally, note that the volume rate (209 Hz) of our system limits the ability of the ULM
pipeline to track microbubbles traveling in the larger cerebral vessels and to provide an appro-
priate flow velocity map (?). In the future, a 3D ultrasound scanner with a higher volume rate
could alleviate this issue and provide additional physiological information.
<figure><img src="image_9.png" /><figcaption>Figure 7: **Confronting ULM density maps to MRA images.** (A,B) Coronal slices and sagittal slices, respectively: The ULM sections are shown as maximum intensity projections of 5 mm (coronal) and 1 cm (sagittal) thickness. Arbitrary colormaps are used with fixed dynamic range for the ULM (hot colormap) and MRA (B&W scale) images to facilitate their comparison. (C,D) Same images as in panels A and B with ULM images corrected by UMI. Green and blue arrows show a better reconstruction of the posterior cerebral arteries in deviation from the Willis polygon. Yellow arrows show the correction in the venous system. The data shown here correspond to sheep number 6 (acquisition 4).</figcaption></figure>

DISCUSSION

Compared to a previous study on a static brain phantom (?), we have here demonstrated the benefit of UMI for transcranial in-vivo imaging. Besides improving the contrast and resolution of standard ultrasound brain images, we also showed that it can be fruitfully coupled to other imaging modalities such as ULM. While the gain in contrast provided by UMI leads to increase the overall number of detected microbubbles, it also drastically increased the length of the microbubble tracks, which is affected by SNR, point-source separability and localiza-
tion precision. Because UMI reduces the side lobes of the imaging PSF, the localization of microbubbles is actually much sharper. Doing so, UMI provides a more detailed ULM image and the removal of common artifacts such as duplication of vessels. The success of UMI has been confirmed by confronting ULM images with gold standard MRA.

More fundamentally, UMI also showed some correlation between the aberration level and the skull thickness variations. It also highlighted the impact of diploë volume fraction on the multiple scattering level. The striking result of UMI on sheep n°6 is promising for applications of UMI on human brain imaging, since this skull heterogeneity is in line with what is expected for human brain imaging across the temporal window (?, ?).

UMI approach for aberration compensation is perfectly complementary with the alternative route proposed recently by Robin *et al*. (?) that exploited bubbles as ultrasound guide stars for trans-cranial adaptive focusing. To some extent, UMI is more robust since it directly synthesizes virtual guide stars from speckle (?) and can work whatever the concentration of bubbles as well as in stronger aberration and multiple scattering conditions. Typically, we were not able to apply the approach proposed in Ref. (?) due to the low quality of the ULM data. Indeed, the latter method requires the presence of bubbles whose associated focal spot displays a correlation degree $\gamma$ larger than 0.6 with the ideal diffraction-limited PSF. Here, the normalized Strehl ratio, whose value is a maximal bound for $\gamma$, ranges from 0.1 to 0.4 (Section ??). Such an aberration level makes difficult the use of microbubbles as guide stars in an adaptive focusing scheme. In Section ??, we also show how UMI provides a flexible benchmarking framework to show the benefits of the iterative phase reversal process (?) compared with conventional aberration correction methods in speckle such as iterative time reversal (?) or SVD beamforming (?, ?).

Nevertheless, in the current realization, UMI also showed some limits. First, it requires the prior recording of a high-dimension reflection matrix. Second, in the present case, it only provided the depth-dependence of aberration phase laws and was not able to capture the lateral
variations of aberrations. Last but not least, it did not compensate for the multiple reverbera-
tions inside the skull that prevented us from brain imaging right below the skull. Those limits
are explained by the fact that high-order aberrations and reverberations are associated with ex-
tremely small isoplanatic patches. This feature limits the efficiency of the spatial averaging
process required for the synthesis of a coherent guide star in matrix imaging (Fig. ??).

In future works, these limits will be addressed by exploiting the dynamics of blood vessels.
Indeed, rapid decorrelation of speckle may offer an opportunity as it provides numerous speckle
realizations for a single voxel. Therefore, the aberration phase laws could in principle be ex-
tracted at a higher spatial resolution, since UMI would no longer have to rely on the isoplanicity
assumption (?, ?). In other words, time averaging will replace spatial averaging, and the aber-
ration phase laws can be extracted directly from the ultrafast sequence [Fig. ??D] and not from
a prior static sequence [Fig. ??C]. This access to a large number of disorder realizations can
also be leveraged for tailoring complex spatio-temporal focusing laws capable of harnessing
multiple reverberations.

In the future, we therefore expect that aberration and scattering compensation provided by
UMI could help ULM in the detection and classification of cerebrovascular accidents in humans.
An important goal will be to distinguish between ischemic and hemorrhagic strokes in the early
phase, as has been done in the rat brain (?). Therefore, one perspective of this work is to conduct
a pilot study with patients who have recently suffered a stroke.

In this paper, we have shown the great advantages and relative simplicity of the UMI and
ULM combination. 3D transcranial in-vivo imaging of sheep brain has been performed and
shown to be in excellent agreement with gold standard MRA. Although we have proven this
claim for the specific case of ultrasound data acquired with a multiplexed matrix probe, these
results are much more general. Indeed, UMI can be applied to any other ultrasound probe
(e.g. sparse or RCA array), any other modality (e.g. conventional Doppler or shear wave
elastography) and in either 2D or 3D configurations.

Beyond the specific case of ultrasound, this study paves the way towards the combination of matrix imaging (?, ?) with super-resolution localization techniques (?) in all fields of wave physics, ranging from optical microscopy with fluorescent molecules (?) to seismology for dynamic imaging of glaciers (?, ?).

# MATERIALS AND METHODS

## Experimental procedure

All experiments were performed between October 2022 and January 2023 in accordance with ARRIVE guidelines, the European Directives and the French Legislation on Animal Experimentation (#32738) and approved by the local ethics committee (CENOMEXA, No54). A total of 6 female domestic sheep were included in the study (Ovis aries weighing 35-40 kg and aged 15-16 months). Three of them were used to optimize the experimental protocol, the other three are presented here (Table ??).

The sheep were sedated with ketamine (12 mg/kg) and xylazine (1 mg/kg) for placement of the venous catheter. Propofol (100 µg/kg) was then administered and endotracheal intubation was performed. During the surgical and imaging procedures, anesthesia was maintained by administration of sevoflurane (1.5-2.0 %) in air. Cardiac activity, oximetry and arterial pressure (between 80 and 120 mm Hg) were monitored throughout the procedure. Respiratory rate and tidal volume were adjusted to maintain physiologic limits, and body temperature was maintained at 38 °C with a warming blanket.

At the end of the ultrasound examination, i.e. after 5 hours, the anesthesia was terminated so that the animal could wake up and be returned to its box. Three imaging sessions were performed within one to two weeks. At the end of the last imaging session, the sheep was euthanized with a barbiturate overdose (pentobarbital 0.2g/kg) after 2 minutes of 8% sevoflurane
inhalation. A controlled arterial pressure of less than 20 mmHg and no cardiac activity were
used to assess death.

After the last imaging session and the death of the animal, the part of the skull below the
probe was removed for cleaning and preserved in saline and bleach. This subject was then fixed
in a bath of formalin-free fixative F13 (Morphisto, Germany) for one day, then dried and placed
under vacuum for micro-CT imaging (Skyscan 1176, Bruker, Kontich, Belgium, voxel size 18
µm) (?).

Ultrasound sequences

A low-frequency matrix array probe (Vermon, Tours, France), as described in Table ??, is placed
on the shaved head of the sheep directly over the crest of the frontal bone. 3D ultrasound
sequences were then recorded with a Vantage 256 4-to-1 multiplex system (Verasonics UTA
1024-MUX adapter, Kirkland, USA). The 1024 elements were connected to the 256 channels
of the ultrasound scanner using a 4-to-1 multiplexer (Verasonics UTA 1024-MUX Adapter).
The probe is therefore divided into four panels, the multiplexer only allowing emission with
one panel and receiving with another (? , ?). For the reflection matrix acquisition, a set of 16
emissions had to be performed for any incident wave-field in order to record the 4 × 4 responses
between each panel at emission and reception. Those recordings are then reassembled in post-
processing to get the full reflected-wave-front associated with each synthetic incident wave-
field. For the ULM sequence, this constraint has been partially mitigated since it is possible
to emit with the four connected panels at the same time, provided that the emission waveform
and delays were identical along the elevation axis (? , ?). As detailed further, the use of incident
cylindrical wave-fronts was therefore used to reach a sufficient frame rate for ULM. The suc-
cessive ultrasound imaging sequences of the entire sheep experiment are summarized in Fig. ??.
<table><caption>Table 2: Matrix array datasheet.</caption><tr><td>Number of transducers</td><td>32 × 32 = 1024 (with 5 dead elements)</td></tr><tr><td>Geometry (y-axis)</td><td>3 inactive rows between each block of 32x8 elements</td></tr><tr><td>Pitch</td><td>δu = 0.55 mm (≈ λ/2 at c = 1540 m/s)</td></tr><tr><td>Aperture</td><td>Δu = (Δu<sub>x</sub>) = (17.6 mm)<br>Δu<sub>y</sub> = (19.3 mm)</td></tr><tr><td>Central frequency</td><td>f<sub>c</sub> = 1.56 MHz</td></tr><tr><td>Bandwidth</td><td>60% → [1 − 2] MHz</td></tr><tr><td>Transducer directivity</td><td>θ<sub>max</sub> = 64° at c = 1540 m/s</td></tr></table>

Initial real-time power Doppler imaging is used to correctly position the probe before acquiring longer ultrasound data. This sequence is performed with the injection of microbubble contrast agent (Sonovue, Bracco, Italy) (?) to obtain SNR-enhanced images of the cerebral vasculature. Once a sufficiently good signal is obtained, the probe is fixed, and we wait until the microbubbles has disappeared, i.e. 10 minutes.

A UMI sequence is then performed. It consists in the acquisition of the reflection matrix using a set of 324 synthetic spherical waves whose associated virtual sources are located along a sphere radius $R = 39.5$ mm: $R_{uv}(t) = [R(u_{out}, v_{in}, t)]$, where the vector $v_{in}$ indicates the position of the corresponding virtual sources (see Supplementary ??) and $u_{out}$ refers to the position of the transducers. The insonification of the medium is made with diverging waves associated with virtual sources located above the probe (?) such that:

$$v = \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix} = -\mathcal{R} \begin{pmatrix} \sin \theta_x \\ \sin \theta_x \\ \sqrt{1 - \sin^2 \theta_x - \sin^2 \theta_y} \end{pmatrix}, \quad (3)$$

with $\theta = (\theta_x, \theta_y)$ the angle associated with each virtual source. The angular range spans from $-\sin \theta_m$ to $+\sin \theta_m$, with $\theta_m = 30^\circ$. The angular pitch $\delta(\sin \theta)$ is governed by the probe extension: $\delta(\sin \theta) = \lambda_c/\Delta u$, with $\lambda_c = c_0/f_c$, the central wave length. This angular range and angular pitch dictate the number of incident diverging waves: $N = (2 \sin \theta_m/\delta(\sin \theta))^2 \sim 324$.
Finally, a hybrid ultrasound localization microscopy (ULM) sequence is performed. To limit the number of incident waves, only the azimuth and elevation shifts ($\theta_x, \theta_y$) of (0,0); (-5, 0); (5, 0); (0, -5); (0, 5) degrees are considered. As stated above, the 4-to-1 multiplexer imposes a sequence of 16 insonifications to record the response of the medium for each virtual source. Such a sequence is prohibitive for high-frame rate imaging. However, this condition can be relaxed by only recording the responses between neighboring panels, thereby restricting the number of insonifications to 10 for each virtual source (?). Moreover, the four panels can be connected together for the emission if the delays are similar in the elevation direction. For such an incident cylindrical wave-front, only 4 successive transmissions are required for the 4 reception apertures. However, full cylindrical waves could only be induced when the angle $\theta_x$ is 0 in the azimuthal direction. In the y-direction, waves still needed to be split in multiple emission-reception. Consequently, an hybrid sequence was chosen to optimize the frame rate while maintaining a sufficient angular diversity for the acquisition:

* Three cylindrical waves were emitted sequentially from the entire probe and had to be repeated four times to record the reflected wave-front on each panel.

* Two spherical waves (virtual source) were decomposed by panels, corresponding to 10 emissions.

The full ULM sequence thus corresponds to 32 emissions to achieve a volume rate of 209 Hz (?,
?). Since only a few insonifications are used, the ULM acquisition consists of a partial reflection
matrix that is acquired at multiple times $t_m$ so that the reflected wavefronts are stored in a
second matrix $R' = [R'(u_{\text{out}}, i_{\text{in}}, \tau, t_m)]$. The sparse feature of the recorded reflection matrices
implies an ultrasound image of lower quality compared to a full measurement of the $R$-matrix.
Both the contrast and, to a lessen extent, the resolution of the beamformed image are degraded.
However, those detrimental effects are compensated, at least partially, by the echogenicity of the
bubbles for the contrast and by the localization feature that allows the resolution to be no longer limited by diffraction but by the signal-to-noise ratio (?). During this acquisition, microbubble injection boluses (Sonovue, Bracco, Italy) were injected, i.e. from 1.5 to 2 ml/min, during 2.5 to 6.5 minutes.

A total of 1 to 6 acquisitions were performed per sheep. The experimental conditions are summarized in Table ?? . The raw ULM data used in the paper corresponds to the ULM data presented in Ref. (?). It should also be noted that, for some acquisition, the UMI sequence is performed after the ULM sequence and not before, but this does not change the principle of the method.

Focused reflection matrix

A general procedure to build each focused reflection matrix, noted R$_{ρρ}$ ( z ), is to use a delay- and-sum (DAS) algorithm, which consists of applying appropriate time delays to create a double synthetic focusing procedure at both input and output (Bureau 2023) (?). Mathematically, such a beamforming procedure is described by Eq. ?? . The input and output times-of-flight, τ$_{in}$ and τ$_{out}$ , are given by:

$$
\tau_{\text{in/out}}(\mathbf{w}, \mathbf{r}) = \frac{|\mathbf{w} - \mathbf{r}|}{c_0} = \frac{\sqrt{(x - w_x)^2 + (y - w_y)^2 + (z - w_z)^2}}{c_0}. \quad (4)
$$

with w = v$_{in}$ at input and w = u$_{out}$ at the output.

The confocal volume shown in Fig. ?? e can be extracted from the diagonal elements of the reflection matrix ( ρ$_{in}$ = ρ$_{out}$ ), so that :

$$
\mathcal{I}(\rho, z) = |R(\rho, \rho, z)|^2 \qquad (5)
$$
Reflection point spread function

The off-diagonal points in R<sub>ρρ</sub>(z) can be exploited for a quantification of the focusing quality at any pixel of the ultrasound image. To that aim, the spreading of energy is investigated along each antdiagonal of R<sub>ρρ</sub>(z). Mathematically, this can be done by means of a change of variable in order to express the focused R-matrix in a common mid-point basis (?):

$$
R(\Delta\rho, \rho_m, z) = R(\rho_{\text{out}}, \rho_{\text{in}}, z) \qquad (6)
$$

with $\rho_m = (\rho_{in} + \rho_{out})/2$, the mid-point between $\rho_{in}$ and $\rho_{out}$, and $\Delta\rho = \rho_{out} - \rho_{in}$, their relative position. A local average of the back-scattered intensity is then performed needed to smooth out the speckle reflectivity and provide an estimator of the reflection point spread function (RPSF):

$$
RPSF(\Delta\rho, \mathbf{r}_p) = \langle |R(\Delta\rho, \rho_m, z)|^2 \mathcal{P}(\rho_m - \rho_p, z - z_p) \rangle_{\rho_m, z} \quad (7)
$$

with $\mathbf{r}_p = (\rho_p, z_p)$. The symbol $\langle \cdots \rangle$ denotes a spatial average over the variables in the subscript. $\mathcal{P}(\rho_m - \rho_p, z - z_p) = 1$ for $|\rho_m - \rho_p| < p_\rho/2$ and $|z_m - z_p| < p_z/2$, and zero otherwise.

The dimensions of $\mathcal{P}$ used for Figs. ??c,d and ??e-h are $\mathbf{p} = (p_\rho, p_z) = (4.3, 3)$ mm.

As shown by Figs. ??c,d and Figs. ??e-h, the RPSF displays the following typical shape:
A distorted confocal spot associated with singly-scatterered echoes on top of a flat background
induced by multiple scattering events. On the one hand, the extension $\delta\rho$ of the confocal peak
is directly related to the local resolution of the ultrasound image. $\delta\rho$ can be obtained by consid-
ering the de-scanned area $A_{-3dB}$ over which the normalized RPSF is larger than -3dB (Supple-
mentary Fig. ??), such that

$$
\delta\rho = \sqrt{A_{-3dB}/\pi}. \tag{8}
$$

On the other hand, the multiple scattering rate $\alpha_M$ can be estimated by averaging the RPSF
beyond the confocal peak, i.e for $4\delta\rho_0 < ||\Delta\rho|| < 6\delta\rho_0$:

$$ \alpha_M(\mathbf{r}_p) = \frac{\langle RPSF(\Delta\rho, \mathbf{r}_p) \rangle_{4\delta\rho_0 < ||\Delta\rho|| < 6\delta\rho_0}}{RPSF(\Delta\rho = \mathbf{0}, \mathbf{r}_p)} \quad (9) $$

## Aberration phase law

An aberrations law is extracted using the distortion matrix framework (?, ?) following these different steps:

- (i) We first project $\mathbf{R}_{\rho\rho}(z)$ onto the transducer basis ($\mathbf{u}$) either at input or output using a homogeneous propagation model:

$$ \mathbf{R}_{\mathbf{u}\rho}(z) = \mathbf{G}_0(z) \times \mathbf{R}_{\rho\rho}(z), \quad (10) $$

and

$$ \mathbf{R}_{\rho u}(z) = \mathbf{R}_{\rho\rho}(z) \times \mathbf{G}_0(z)^{\top} \quad (11) $$

with $\mathbf{G}_0 = [G_0(\mathbf{u}, \rho, z)]$, the transmission matrix between the probe and the set of voxels at depth $z$ for an ideal homogeneous medium:

$$ T_0(\mathbf{u}, \rho, z) = \frac{z \exp \left( j k_0 \sqrt{\|\mathbf{u} - \rho\|^2 + z^2} \right)}{4\pi (\|\mathbf{u} - \rho\|^2 + z^2)} \quad (12) $$

As sketched in Fig. ??B, this dual reflection matrix contains the wave-fronts associated with virtual sources corresponding to the focusing points ($\rho_{in/out}, z$). Those wave-fronts correspond to the sum of a geometric parabolic component related to the position of each focusing point ($\rho_{in/out}, z$) and of a distorted component resulting from the mismatch between the propagation model and the real speed-of-sound distribution. Interestingly,
this distorted component displays strong correlations because each reflected wave-front goes through the same part of the skull, the probe being placed at the skull surface. This is the so-called angular (or tilt-tilt) memory effect (?) that we will exploit in the following.

- (ii) To do so, we isolate the distorted part of each wavefront by subtracting the geometric component that would be ideally obtained in absence of aberrations (i.e. constant speed of sound $c_0$). The result is the input and output distortion matrices which are expressed as follows:

$$D_{u\rho}(z) = R_{u\rho}(z) \circ \bar{T}_0 \quad (13)$$

and,

$$D_{\rho u}(z) = R_{\rho u}(z) \circ \bar{T}_0^\top \quad (14)$$

where the symbol $\circ$ stands for the Hadamard (element wise) product and with, $\bar{T}_0$, the normalized transmission matrix such that

$$\bar{T}_0(\mathbf{u}, \rho, z) = T_0(\mathbf{u}, \rho, z) / |T_0(\mathbf{u}, \rho, z)|.$$

In the focused basis, the operation described in Eqs. ?? and ?? amounts to a de-scan of each virtual source ($\rho_{in/out}, z$) at the origin $(0, z)$ of each focal plane, as sketched in Fig. ??C. In the transducer basis, it consists in a realignment of associated wave-fronts, thereby highlighting the angular correlations between wave distortions.

- (iii) We then exploit this angular memory effect by computing the correlation matrices of wave-front distortions,

$$C_{in}(z) = D_{u\rho}(z) \times D_{u\rho}^{\dagger}(z), \quad (15)$$

and

$$C_{out}(z) = D_{\rho u}^{T}(z) \times D_{\rho u}^{*}(z). \quad (16)$$
(iv) A numerical iterative phase reversal (IPR) algorithm (?) is performed to extract an estimation of aberration laws $\Phi_{in/out}(z) = [\phi_{in/out}(\mathbf{u}_{in/out}, z)]$ by solving the following equation:

$$ \Phi_{in/out}(z) = \arg \{C_{in/out} \times \exp(i\Phi_{in/out}(z))\}. \quad (17) $$

The IPR process, sketched in Fig. ??D, allows the synthesis of a coherent guide star from the set of virtual sources ($\rho_{in/out}, z$). The associated wave-front is the resulting aberration phase law $\Phi_{out/in}(z)$. The complexity of each aberration law can be quantified by the Strehl ratio which is defined as follows (?) (see also Section ??):

$$ S(z) = |\langle\exp[i\phi(\mathbf{u}, z)]\rangle_{\mathbf{u}}|^2. \quad (18) $$

The aberration phase laws yield an estimator of the input/output transmission matrix between the probe and each focal plane:

$$ T_{in/out}(z) = \exp(i\Phi_{out}) \circ T_0 \quad (19) $$

(vi) Aberrations are compensated by applying the phase conjugate of the transmission matrices to provide an updated focused reflection matrix:

$$ \mathbf{R}_{\rho\rho}(z) = \mathbf{T}_{\mathrm{in}}^{\dagger}(z) \times \mathbf{R}_{\rho\mathbf{u}}(z), \quad (20) $$

and

$$ \mathbf{R}_{\rho\rho}(z) = \mathbf{R}_{\rho\mathbf{u}}(z) \times \mathbf{T}_{\mathrm{out}}^{*}(z), \quad (21) $$

The whole process is repeated two times at input and output to improve gradually the focusing process. The whole UMI process is summarized by the flow chart displayed in Fig. ???. The confocal volume (Eq. ??) displayed in Fig. ??f and the map of RPSFs (Eq. ??) shown in Fig. ??d are extracted from the diagonal coefficients of this corrected focused reflection matrix.
Ultrasound localization microscopy

From the reflection matrix R' recorded using the hybrid sequence, a ULM image can be built
by performing the following post-processing steps.

- (i) High volume rate ultrasound images are constructed using a conventional delay-and-sum algorithm into a 300×300×600 mm³ volume with 0.5×0.5×0.5 mm³ voxels (?):

$$
\mathcal{I}'(\rho, z, t_m) = \sum_{\mathbf{s}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} R'(\mathbf{u}_{\text{out}}, \mathbf{s}_{\text{in}}, \tau_{\text{in}}(\mathbf{s}_{\text{in}}, \rho, z) + \tau_{\text{out}}(\mathbf{u}_{\text{out}}, \rho, z), t_m)
$$

- (ii) A high-pass time filter is applied to the set of images by means of a two frames sliding average window:

$$
\mathcal{I}''(\rho, z, t_m) = \mathcal{I}'(\rho, z, t_m) - \langle \mathcal{I}'(\rho, z, t_n) \rangle_{n \in [m-1;m+1]}
$$

Albeit crude, this simple operation allows the removal of most of the slowly-varying
component of the ultrasound image (i.e tissues) while enhancing the dynamic part of the
image (moving microbubbles).

- (iii) The microbubbles were then detected in each volume by localizing regional maxima on each image $\mathcal{I}''$. The center of each PSF was sub-localized using the radial symmetry method (?, ?, ?), i.e. by assuming a PSF Gaussian shape and by applying a gradient algorithm. Microbubbles with a signal-to-noise ratio (SNR) below 10.5 dB were discarded. This SNR was calculated as the intensity at the center of the microbubbles divided by the mean intensity of the neighboring voxels over a 5-mm-wide patch. The detected bubbles and their localization are superimposed to one cross-section section of the ultrasound image in Figs. ??a,b.

- (iv) Tracking algorithm has been performed using the Hungarian method (?) to build the
histogram displayed in Fig. ??d. A maximum linking distance between two microbubbles
is fixed to 1.5 mm. A one frame gap is allowed and no restriction on track duration is
imposed.

(v) A 3D ULM density map $\mathcal{U}(\mathbf{r})$ is obtained by accumulating and projecting all tracks of microbubbles onto a three-dimensional grid whose spatial sampling is $\lambda/10$. The result is displayed in Fig. ??a,c.

The same process is applied to the ultrasound images $\mathcal{I}''$ corrected by UMI (Eq. ??). The resulting ULM maps are displayed in Figs. ??b,d.

## Computational Insights

The computation time of the UMI process (?) is limited by the beamforming of the focused R-matrix (Eq. ??) that takes here 30 min (on GPU with CUDA language) while the aberration correction process only lasts for a few minutes. All the post-processing was realized with Matlab (R2024a) on a working station with 2 processors @2.20GHz, 128Go of RAM, and a GPU with 48 Go of dedicated memory. This computation time for UMI remains acceptable compared to the ULM process that takes 1h30min for 9000 frames. It should be noted that the UMI process (beamforming and correction) can be further speed up by reducing the number of focusing points required to extract the aberration laws.

## Angiographic MRI sequences

Angiographic MRI sequences (MRA) were performed with a 3T MRI scanner (GE SIGNA PREMIER) in a thermoregulated room, i.e. at approximately 20°C (axial propeller, TR=5539 ms, TE=165.744 ms, slice thickness 2 mm, 31 slices, in-place resolution 0.375×0.375). The superimposition of MRA and ULM images in Fig. ?? has been performed manually using Amira
software (Thermofischer, v2019.4).

## SUPPLEMENTARY MATERIALS

Section ???. Virtual source basis.

Section ???. Flow chart of the UMI process.

Section ???. Magnitude of aberration correction.

Section ???. Resolution enhancement.

Section ???. Multiple scattering rate.

Section ???. Comparison with other aberration correction methods.

Section ???. Quantitative comparison between MRA and ULM images.

Section ???. ULM images for each dataset.

Fig. ???. Virtual source basis.

Fig. ???. Flow chart of the UMI process.

Fig. ???. Quantifying the magnitude of aberration correction.

Fig. ???. Quantifying the ultrasound image resolution.

Fig. ???. Single and multiple scattering rates.

Fig. ???. Comparison with other aberration correction methods.

Fig. ???. Quantitative comparison between ULM and MRA cross-sections in the sagittal plane.

Fig. ???. Quantitative comparison between ULM and MRA cross-sections in the coronal plane.

Fig. ???. Impact of UMI on different ULM acquisitions.

Tab. ???. Experimental conditions for each acquisition.

Movie S1. Result of UMI at each depth

Movie S2. ULM images for different acquisitions before and after UMI correction.

References (?, ?).
# REFERENCES AND NOTES

## References

1. M. Bahrami, M. Keyhanifard, M. Afzali, Spontaneous intracerebral hemorrhage, initial computed tomography (ct) scan findings, clinical manifestations and possible risk factors. *Am. J. Nucl. Med. Mol. Imaging* **22**, 106–112 (2022).

2. L. Scheldeman, et al., Different mismatch concepts for magnetic resonance imaging-guided thrombolysis in unknown onset stroke. *Annals of Neurology* **87**, 931–938 (2020).

3. J. D. Kirsch, M. Mathur, M. H. Johnson, G. Gowthaman, L. M. Scoutt, Advances in tran-scranial Doppler US: Imaging ahead. *Radiographics* **33**, E1–E14 (2013).

4. M. Versluis, E. Stride, G. Lajoinie, B. Dollet, T. Segers, Ultrasound Contrast Agent Mod-eling: A Review. *Ultrasound Med. Biol.* **46**, 2117–2144 (2020).

5. A. Tarighatnia, M. R. Fouladi, N. D. Nader, A. Aghanejad, H. Ghadiri, Recent trends of contrast agents in ultrasound imaging: A review of the classifications and applications. *Mater. Adv.* **3**, 3726–3741 (2022).

6. H.-L. Liu, C.-H. Fan, C.-Y. Ting, C.-K. Yeh, Combining Microbubbles and Ultrasound for Drug Delivery to Brain Tumors: Current Progress and Overview. *Theranostics* **4**, 432–444 (2014).

7. O. Couture, S. Bannouf, G. Montaldo, J.-F. Aubry, M. Fink, M. Tanter, Ultrafast Imaging of Ultrasound Contrast Agents. *Ultrasound Med. Biol.* **35**, 1908–1916 (2009).

8. O. Couture, M. Fink, M. Tanter, Ultrasound contrast plane wave imaging. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **59**, 2676-2683 (2012).
9. C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture, M. Tanter, Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging. *Nature* **527**, 499–502 (2015).

10. O. Couture, V. Hingot, B. Heiles, P. Muleki-Seya, M. Tanter, Ultrasound Localization Microscopy and Super-Resolution: A State of the Art. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **65**, 1304–1320 (2018).

11. K. Christensen-Jeffries, O. Couture, P. A. Dayton, Y. C. Eldar, K. Hynynen, F. Kiessling, M. O'Reilly, G. F. Pinton, G. Schmitz, M.-X. Tang, M. Tanter, R. J. G. van Sloun, Super-resolution Ultrasound Imaging. *Ultrasound Med. Biol.* **46**, 865–891 (2020).

12. P. Song, J. M. Rubin, M. R. Lowerison, Super-resolution ultrasound microvascular imaging: Is it ready for clinical use? *Z. Med. Phys.* **33**, 309–323 (2023).

13. E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, H. F. Hess, Imaging intracellular fluorescent proteins at nanometer resolution. *Science* **313**, 1642–1645 (2006).

14. M. J. Rust, M. Bates, X. Zhuang, Stochastic optical reconstruction microscopy (STORM) provides sub-diffraction-limit image resolution. *Nat. Methods* **3**, 793–795 (2006).

15. Y. Desailly, J. Pierre, O. Couture, M. Tanter, Resolution limits of ultrafast ultrasound localization microscopy. *Phys. Med. Biol.* **60**, 8723–8740 (2015).

16. V. Hingot, C. Errico, M. Tanter, O. Couture, Subwavelength motion-correction for ultrafast ultrasound localization microscopy. *Ultrasonics* **77**, 17–21 (2017).
17. M. R. Lowerison, N. V. C. Sekaran, W. Zhang, Z. Dong, X. Chen, D. A. Llano, P. Song, Aging-related cerebral microvascular changes visualized using ultrasound localization mi-croscopy in the living mouse. *Sci. Rep.* **12**, 619 (2022).

18. J. Claron, V. Hingot, I. Rivals, L. Rahal, O. Couture, T. Deffieux, M. Tanter, S. Pezet, Large-scale functional ultrasound imaging of the spinal cord reveals in-depth spatiotemporal responses of spinal nociceptive circuits in both normal and inflammatory states. *Pain* **162**, 1047–1059 (2021).

19. O. Demeulenaere, Z. Sandoval, P. Mateo, A. Dizeux, O. Villemain, R. Gallet, B. Ghaleh, T. Deffieux, C. Deméné, M. Tanter, C. Papadacci, M. Pernot, Coronary Flow Assessment Using 3-Dimensional Ultrafast Ultrasound Localization Microscopy. *JACC Cardiovasc. Imaging* **15**, 1193–1208 (2022). Place: United States.

20. K. B. Hansen, C. A. Villagómez-Hoyos, J. C. Brasen, K. Diamantis, V. Sboros, C. M. Sørensen, J. A. Jensen, 2016 IEEE International Ultrasonics Symposium (IUS) (2016), pp. 1–4. ISSN: 1948-5727.

21. I. Taghavi, S. B. Andersen, C. A. V. Hoyos, M. B. Nielsen, C. M. Sørensen, J. A. Jensen, In Vivo Motion Correction in Super-Resolution Imaging of Rat Kidneys. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **68**, 3082–3093 (2021).

22. J. Zhu, et al., 3D Super-Resolution US Imaging of Rabbit Lymph Node Vasculature in Vivo by Using Microbubbles. *Radiology* **291**, 642–650 (2019).

23. F. Lin, S. E. Shelton, D. Espíndola, J. D. Rojas, G. Pinton, P. A. Dayton, 3-D Ultrasound Localization Microscopy for Identifying Microvascular Morphology Features of Tumor Angiogenesis at a Resolution Beyond the Diffraction Limit of Conventional Ultrasound. *Theranostics* **7**, 196–204 (2017).
24. T. Opacic, S. Dencks, B. Theek, M. Piepenbrock, D. Ackermann, A. Rix, T. Lammers, E. Stickeler, S. Delorme, G. Schmitz, F. Kiessling, Motion model ultrasound localization microscopy for preclinical and clinical multiparametric tumor characterization. *Nat. Commun.* **9**, 1527 (2018).

25. C. Huang, W. Zhang, P. Gong, U.-W. Lok, S. Tang, T. Yin, X. Zhang, L. Zhu, M. Sang, P. Song, R. Zheng, S. Chen, Super-resolution ultrasound localization microscopy based on a high frame-rate clinical ultrasound scanner: an in-human feasibility study. *Phys. Med. Biol.* **66** (2021).

26. C. Demené, J. Robin, A. Dizeux, B. Heiles, M. Pernot, M. Tanter, F. Perren, Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients. *Nat. Biomed. Eng.* **5**, 219–228 (2021).

27. L. Denis, S. Bodard, V. Hingot, A. Chavignon, J. Battaglia, G. Renault, F. Lager, A. Aissani, O. Hélénon, J.-M. Correas, O. Couture, Sensing ultrasound localization microscopy for the visualization of glomeruli in living rats and humans. *eBioMedicine* **91**, 104578 (2023).

28. S. Bodard, L. Denis, V. Hingot, A. Chavignon, O. Hélénon, D. Anglicheau, O. Couture, J.-M. Correas, Ultrasound localization microscopy of the human kidney allograft on a clinical ultrasound scanner. *Kidney Int.* **103**, 930–935 (2023).

29. S. Bodard, L. Denis, G. Chabouh, J. Battaglia, D. Anglicheau, O. Hélénon, J.-M. Correas, O. Couture, Visualization of Renal Glomeruli in Human Native Kidneys With Sensing Ultrasound Localization Microscopy. *Investig. Radiol.* **59**, 561–568 (2024).

30. F. Knieling, et al., Ultrasound super-resolution imaging of neonatal cerebral vascular reorganization during neurovascular interventions. *Research Square* (2023). Research Square.
31. J. Yan, B. Huang, J. Tonko, M. Toulemonde, J. Hansen-Shearer, Q. Tan, K. Riemer, K. Ntagiantas, R. A. Chowdhury, P. D. Lambiase, R. Senior, M.-X. Tang, Transthoracic ultrasound localization microscopy of myocardial vasculature in patients. *Nat. Biomed. Eng.* **8**, 689–700 (2024).

32. L. Denis, E. Meseguer, A. Gaudemer, G. Jaklh, S. Bodard, G. Chabouh, D. Hervé, E. Vicaut, P. Amarenco, O. Couture, Transcranial ultrasound localization microscopy in moyamoya patients using a clinical ultrasound system (2024).

33. D. E. Soulioti, D. Espindola, P. A. Dayton, G. F. Pinton, Super-Resolution Imaging Through the Human Skull. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **67**, 25–36 (2020).

34. J. Robin, C. Demené, B. Heiles, V. Blanvillain, L. Puke, F. Perren, M. Tanter, In vivo adaptive focusing for clinical contrast-enhanced transcranial ultrasound imaging in human. *Phys. Med. Biol.* **68**, 025019 (2023).

35. B. Heiles, M. Correia, V. Hingot, M. Pernot, J. Provost, M. Tanter, O. Couture, Ultrafast 3D Ultrasound Localization Microscopy Using a 32×32 Matrix Array. *IEEE Trans. Med. Imag.* **38**, 2005–2015 (2019).

36. A. Chavignon, V. Hingot, C. Orset, D. Vivien, O. Couture, 3D transcranial ultrasound localization microscopy for discrimination between ischemic and hemorrhagic stroke in early phase. *Sci. Rep.* **12**, 1–11 (2022).

37. O. Demeulenaere, A. Bertolo, S. Pezet, N. Ialy-Radio, B. Osmanski, C. Papadacci, M. Tanter, T. Deffieux, M. Pernot, In vivo whole brain microvascular imaging in mice using transcranial 3D Ultrasound Localization Microscopy. *EBioMedicine* **79**, 103995 (2022).

38. J. R. McCall, R. DeRuiter, M. Ross, F. Santibanez, S. D. Hingtgen, G. F. Pinton, P. A. Dayton, Longitudinal 3-D Visualization of Microvascular Disruption and Perfusion Changes
in Mice During the Evolution of Glioblastoma Using Super-Resolution Ultrasound. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **70**, 1401–1416 (2023).

39. C. Bourquin, *et al.*, Quantitative pulsatility measurements using 3D dynamic ultrasound localization microscopy. *Phys. Med. Biol.* **69**, 045017 (2024).

40. S. Lei, C. Zhang, B. Zhu, Z. Gao, Q. Zhang, J. Liu, Y. Li, H. Zheng, T. Ma, In vivo ocular microvasculature imaging in rabbits with 3D ultrasound localization microscopy. *Ultrasonics* **133**, 107022 (2023).

41. L. Wei, *et al.*, High-Frame-Rate Volumetric Porcine Renal Vasculature Imaging. *Ultra-sound Med. Biol.* **49**, 2476–2482 (2023).

42. A. Coudert, L. Denis, A. Chavignon, S. Bodard, M. Naveau, P. P. Sistiaga, R. Saulnier, C. Orset, D. Vivien, C. Chappard, O. Couture, 3D Transcranial ultrasound localization microscopy reveals major arteries in the sheep brain. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* (2024).

43. P. Xing, V. Perrot, A. U. Dominguez-Vargas, S. Quessy, N. Dancause, J. Provost, Towards Transcranial 3D Ultrasound Localization Microscopy of the Nonhuman Primate Brain (2024). ArXiv:2404.03547.

44. G. Chabouh, L. Denis, S. Bodard, F. Lager, G. Renault, A. Chavignon, O. Couture, Whole organ volumetric sensing Ultrasound Localization Microscopy for characterization of kidney structure. *IEEE Trans. Med. Imag.* **PP** (2024).

45. P. Xing, J. Porée, B. Rauby, A. Malescot, E. Martineau, V. Perrot, R. L. Rungta, J. Provost, Phase Aberration Correction for in vivo Ultrasound Localization Microscopy Using a Spatiotemporal Complex-Valued Neural Network. *IEEE Trans. Med. Imag.* **43**, 662-673 (2024).
46. H. W. Babcock, The possibility of compensating astronomical seeing. *Publications of the Astronomical Society of the Pacific* **65**, 229–236 (1953).

47. A. Labeyrie, Attainment of Diffraction Limited Resolution in Large Telescopes by Fourier Analysing Speckle Patterns in Star Images. *Astronomy and Astrophysics* **6**, 85 (1970).

48. R. A. Muller, A. Buffington, Real-time correction of atmospherically degraded telescope images through image sharpening. *J. Opt. Soc. Am.* **64**, 1200–1210 (1974).

49. F. Roddier, ed., *Adaptive Optics in Astronomy* (Cambridge University Press, Cambridge, 1999).

50. M. Hirama, O. Ikeda, T. Sato, Adaptive ultrasonic array imaging system through an inhomogeneous layer. *J. Acoust. Soc. Am.* **71**, 100–109 (1982).

51. L. Nock, G. E. Trahey, S. W. Smith, Phase aberration correction in medical ultrasound using speckle brightness as a quality factor. *J. Acoust. Soc. Am.* **85**, 1819–1833 (1989).

52. S. Flax, M. O'Donnell, Phase-aberration correction using signals from point reflectors and diffuse scatterers: basic principles. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **35**, 758–767 (1988).

53. R. Mallart, M. Fink, Adaptive focusing in scattering media through sound-speed inhomogeneities: The van Cittert Zernike approach and focusing criterion. *J. Acoust. Soc. Am.* **96**, 3721–3732 (1994).

54. G. Montaldo, M. Tanter, M. Fink, Time Reversal of Speckle Noise. *Phys. Rev. Lett.* **106**, 054301 (2011).

55. S.-E. Måsøy, T. Varslot, B. Angelsen, Iteration of transmit-beam aberration correction in medical ultrasound imaging. *J. Acoust. Soc. Am.* **117**, 450–461 (2005).
56. M. O'Donnell, S. Flax, Phase-aberration correction using signals from point reflectors and diffuse scatterers: measurements. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **35**, 768–774 (1988).

57. F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, A. Aubry, Three-dimensional ultrasound matrix imaging. *Nat. Commun.* **14**, 6793 (2023).

58. T. Varslot, H. Krogstad, E. Mo, B. A. Angelsen, Eigenfunction analysis of stochastic backscatter for characterization of acoustic aberration in medical ultrasound imaging. *J. Acoust. Soc. Am.* **115**, 3068–3076 (2004).

59. J.-L. Robert, M. Fink, Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging. *J. Acoust. Soc. Am.* **123**, 866–877 (2008).

60. A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, A. Aubry, Distortion matrix concept for deep optical imaging in scattering media. *Sci. Adv.* **6**, aay7170 (2020).

61. W. Lambert, L. A. Cobus, M. Couade, M. Fink, A. Aubry, Reflection Matrix Approach for Quantitative Imaging of Scattering Media. *Phys. Rev. X* **10**, 021048 (2020).

62. W. Lambert, L. A. Cobus, T. Frappart, M. Fink, A. Aubry, Distortion matrix approach for ultrasound imaging of random scattering media. *Proc. Nat. Acad. Sci. USA* **117**, 14645–14656 (2020).

63. W. Lambert, J. Robin, L. A. Cobus, M. Fink, A. Aubry, Ultrasound Matrix Imaging–Part I: The focused reflection matrix, the *f*-factor and the role of multiple scattering. *IEEE Trans. Med. Imag.* **41**, 3907–3920 (2022).
64. F. J. Fry, J. E. Barger, Acoustical properties of the human skull. *J. Acoust. Soc. Am.* **63**, 1576–1590 (1978).

65. J. Provost, C. Papadacci, J. E. Arango, M. Imbault, M. Fink, J.-L. Gennisson, M. Tanter, M. Pernot, 3D ultrafast ultrasound imaging in vivo. *Phys. Med. Biol.* **59**, L1–L13 (2014).

66. W. Lambert, L. A. Cobus, M. Robin, Justine and Fink, A. Aubry, Ultrasound Matrix Imaging–Part II: The distortion matrix for aberration correction over multiple isoplanatic patches. *IEEE Trans. Med. Imag.* **41**, 3921–3938 (2022).

67. V. N. Mahajan, Strehl ratio for primary aberrations: some analytical results for circular and annular pupils. *J. Opt. Soc. Am.* **72**, 1258 (1982).

68. A. Coudert, A. Chavignon, L. Denis, O. Couture, Volumetric ultrasound localization mi-
croscopy with diverging cylindrical waves. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* pp. 1–1 (2024).

69. E. Kalsoum, X. Leclerc, A. Drizenko, J.-P. Pruvo, *Encyclopedia of the Neurological Sciences* (Elsevier, 2014), pp. 803–805.

70. D.-Y. Tsai, Y. Lee, E. Matsuyama, Information entropy measure for evaluation of image quality. *J. Digit. Imaging* **21**, 338–347 (2007).

71. S. Boruah, G. R. Paskoff, B. S. Shender, D. L. Subit, R. S. Salzar, J. R. Crandall, Variation of bone layer thicknesses and trabecular volume fraction in the adult male human calvarium. *Bone* **77**, 120–134 (2015).

72. K. Adanty, K. N. Rabey, M. R. Doschak, K. B. Bhagavathula, J. D. Hogan, D. L. Romanyk, S. Adeeb, S. Ouellet, T. A. Plaisted, S. S. Satapathy, C. R. Dennison, Cortical and trabecular morphometric properties of the human calvarium. *Bone* **148**, 115931 (2021).
73. H. Bendjador, T. Deffieux, M. Tanter, The SVD Beamformer: Physical Principles and Application to Ultrafast Adaptive Ultrasound. *IEEE Transactions on Medical Imaging* **39**, 3100–3112 (2020).

74. H. Bendjador, S. Décombas-Deschamps, M. Dioguardi Burgio, R. Sartoris, B. Van Beers, V. Vilgrain, T. Deffieux, M. Tanter, The svd beamformer with diverging waves: a proof-of-concept for fast aberration correction **66**, 18LT01 (2021).

75. D. Zhao, L. N. Bohs, G. E. Trahey, Phase aberration correction using echo signals from moving targets I: Description and theory. *Ultrasonic Imaging* **14**, 97–110 (1992).

76. B.-F. Osmanski, G. Montaldo, M. Tanter, M. Fink, Aberration correction by time reversal of moving speckle noise. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **59**, 1575–1583 (2012).

77. P. Balondrade, V. Barolle, N. Guigui, E. Auriant, N. Rougier, C. Boccara, M. Fink, A. Aubry, Multi-spectral reflection matrix for ultrafast 3d label-free microscopy. *Nat. Photonics* **18**, 1097–1104 (2024).

78. E. Giraudat, A. Burtin, A. Le Ber, M. Fink, J.-C. Komorowski, A. Aubry, Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise. *Commun. Earth & Environ.* **5**, 509 (2024).

79. S. Park, Y. Jo, M. Kang, J. H. Hong, S. Ko, S. Kim, S. Park, H. C. Park, S.-H. Shim, W. Choi, Label-free adaptive optics single-molecule localization microscopy for whole zebrafish. *Nat. Commun.* **14**, 4185 (2023).

80. U. Nanni, F. Gimbert, P. Roux, A. Lecointre, Observing the subglacial hydrology network and its dynamics with a dense seismic array. *Proc. Nat. Acad. Sci.* **118**, e2023757118 (2021).
81. U. Nanni, P. Roux, F. Gimbert, A. Lecointre, Dynamic imaging of glacier structures at high-resolution using source localization with a dense seismic array. *Geophys. Res. Lett.* **49**, e2021GL095996 (2022).

82. M. Schneider, Characteristics of SonoVue™. *Echocardiography* **16**, 743–746 (1999).

83. G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, I. M. Vellekoop, Generalized optical memory effect. *Optica* **4**, 886 (2017).

84. B. Heiles, A. Chavignon, A. Bergel, V. Hingot, H. Serroune, S. Pezet, M. Pernot, M. Tanter, O. Couture, Volumetric ultrasound localization microscopy of the whole brain microvasculature. in *IEEE Open Journal of Ultrasonics, Ferroelectrics, and Frequency Control* **2**, 261–282 (2022).

85. B. Heiles, A. Chavignon, V. Hingot, P. Lopez, E. Teston, O. Couture, Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy. *Nat. Biomed. Eng* **6**, 605–616 (2022).

86. R. Parthasarathy, Rapid, accurate particle tracking by calculation of radial symmetry centers. *Nat. Methods* **9**, 724–726 (2012).

87. J.-Y. Tinevez, simpletracker [code], https://github.com/tinevez/simpletracker (2019).

88. M. Tanter, J.-L. Thomas, M. Fink, Time reversal and the inverse filter. *J. Acoust. Soc. Am.* **108**, 223–234 (2000).

89. L. Denis, E. Meseguer, A. Gaudemer, G. Jaklh, S. Bodard, G. Chabouh, D. Hervé, E. Vičaut, P. Amarenco, O. Couture, Transcranial ultrasound localization microscopy in moyamoya patients using a clinical ultrasound system. *Theranostics* **15**, 4074–4083 (2025).
90. F. Bureau, L. Denis, A. Coudert, M. Fink, O. Couture, A. Aubry, Ultrasound matrix imaging for 3D transcranial in-vivo localization microscopy [data] (Zenodo). Zenodo.

**Acknowledgements:** The authors wish to thank the laboratory Cyceron, Institut Blood and Brain (Normandie Université, INSERM) where the sheep experiments were initially performed. In particular, we would like to thank Mikael Naveau, Palma Pro Sistiaga, Romaric Saulnier, Cyrille Orset and Denis Vivien. We would also like to thank Christine Chappard for the micro-CT imaging.

**Funding:** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement nos. 772786 and 610110, ResolveStroke and REMINISCENCE projects, OC and AA, respectively).

**Author contributions:**

OC contributed to the conception of the ULM study.

AA initiated and supervised the UMI part of the project.

AC, FB and LD contributed to the conception of the ultrasound sequences.

LD and AC contributed to data acquisition.

FB performed the UMI analysis.

LD and AC performed the ULM analysis.

LD, FB, AC, OC and AA contributed to the interpretation of the experimental data.

FB, LD, AC and AA prepared the figures.

FB, LD and AA contributed to drafting the manuscript.

FB, LD, AC, MF, OC and AA discussed the results and contributed to finalizing the manuscript.

**Competing interests:** AA and MF are inventors on a patent related to UMI held by CNRS (no.
US11346819B2, published May 2022). AC and OC are inventors on a patent related also held by CNRS (no. US20240041423A1, published August 2024). OC is co-founder and shareholder of the startup ResolveStroke. All authors declare that they have no other competing interests.

**Data and materials availability:** The ultrasound data generated in this study are available at Zenodo (?) (https://zenodo.org/records/15148477).
# Supplementary Material

## S1 Virtual source basis

The distribution of virtual sources considered for the acquisition of the reflection matrix (UMI acquisition) is displayed in Fig. ??a. As described in the Methods Section, the virtual sources are located on a sphere centered on the probe midpoint. The time delays applied to each transducer to create the virtual source associated with the white spot in Fig. ??A are shown in Fig. ??B.

<figure><img src="image_10.png" /><figcaption>Figure S1: Virtual source basis. (A,B,C) The R-matrix is acquired with a set of spherical diverging wave associated with virtual sources whose spatial distribution is displayed with blue dots. Each panel corresponds to a different virtual source whose position is indicated by a white disk. (D,E,F) Time delay applied to transducers to generate the diverging wave associated with the virtual source highlighted by a white disk in panels A, B and C, respectively.</figcaption></figure>

## S2 Flow chart of the UMI process
<figure><img src="image_11.png" /><figcaption>Figure S2: Flow chart of the UMI process.</figcaption></figure>
## S3 Magnitude of aberration correction

The impact of skull aberrations can be estimated by the Strehl ratio which directly quantifies the decrease of the PSF. Under a random phase screen model, this Strehl ratio can be shown to scale as $\exp(-\sigma_{\phi}^2)$, with $\sigma_{\phi}$, the standard deviation of the aberration phase law. However, in the present case, the skull transfer function $A = [A(\mathbf{u}, z)]$ is not of constant amplitude. It can actually be estimated through the following matrix product:

$$A(z) = C \times \exp(i\Phi(z)). \qquad (S1)$$

Examples of phase and modulus of the estimated transfer function $A(z)$ are displayed in Figs. ??B and C, respectively. While the phase of $A(z)$ is by definition equal to the aberration phase law $\Phi(z)$, its modulus describes the wave amplitude across the probe aperture. As shown by Fig. ??, this amplitude is clearly not constant over the probe aperture. The evaluation of aberration correction magnitude should therefore include this amplitude variation of the aberration transfer function. To do so, a weighted Strehl ratio can be defined in order to account for both the amplitude and phase of the aberration transfer function:

$$\bar{S}(z) = \frac{|\langle A(\mathbf{u}, z)\phi(\mathbf{u}, z)\rangle_{\mathbf{u}}|^2}{\langle |A(\mathbf{u}, z)\phi(\mathbf{u}, z)|^2 \rangle_{\mathbf{u}}}. \qquad (S2)$$

This normalized Strehl ratio is an important parameter for ULM since it corresponds the correlation degree $\gamma$ between the imaging PSF and its theoretical diffraction-limited value.

The depth-evolution of $\bar{S}(z)$ is displayed for the three sheeps in Fig. ??A. The weighted Strehl ratio is lower for sheep n°6 ($\langle\bar{S}\rangle \sim 0.21$) than for sheeps n°4 and 5 ($\langle\bar{S}\rangle \sim 0.27$ and 0.28, respectively). This lower Strehl ratio indicates a stronger aberration magnitude for sheep n°6. This observation is in agreement with the stronger distortions exhibited by the initial RPSF for sheep n°6 (Fig. ??B) and is explained by the complex structure of its skull (Fig. ??A).
<figure><img src="image_12.png" /><figcaption>Figure S3: **Quantifying the magnitude of aberration correction.** (A) Depth-evolution of the weighted Strehl ratio (Eq. ??) associated with the aberration transfer functions extracted for sheep n°4 (green), n°5 (purple) and n°6 (orange). (B,C,D) Aberration phase laws extracted at depth z = 35 mm for sheeps n° 4, 5 and 6, respectively (scale bar: 5 mm). (E,F,G) Modulus of the aberration transfer function (Eq. ??) estimated at depth z = 35 mm for sheeps n° 4, 5 and 6, respectively.</figcaption></figure>

## S4 Resolution enhancement

The method used for the resolution estimation is illustrated in Fig. ?? by considering an example of RPSF for a specific data set (sheep 6, acquisition n°5) and at a specific depth (z = 36 mm). The RPSF map before and after the correction is shown in Figs. ??A,D. An example of RPSF is displayed in Fig. ??B,C. The -3 dB-area, noted $A$, is evaluated for each RPSF. An estimate of the resolution length $\delta\rho$ is then deduced using Eq. ?? (Fig. ??E). The depth evolution of the resolution $\delta\rho$ is shown for each data set in Fig. ??F. The benefit of UMI in terms of resolution is clear whatever the acquisition. The gain in resolution displayed in Fig. ??G increases from a factor 2 at shallow depths to 3 at large depths.
<figure><img src="image_13.png" /><figcaption>Figure S4: Quantifying the ultrasound image resolution. (A,D) RPSF map at *z* = 36 mm before and after correction, respectively. (B,C) Example of one RPSF at a specific location with the -3 dB-area surrounded by a dashed contour, before and after correction, respectively. (E) Radial average of the RPSFs displayed in panels C (red line) and D (green line). (F) Resolution *δρ* as a function of depth before (red) and after UMI correction (green). The transparent curves refer to the seven individual data sets, while the opaque curves represent their average. (G) Resolution enhancement provided by UMI. The gray lines refer to each individual data set, the black line shows their average value.</figcaption></figure>
S5 Multiple scattering rate

Figure ??A displays the depth evolution of the single and multiple scattering rates (Methods) before and after UMI for the acquisition 4 on sheep n°6. Corresponding RPSF maps and con- focal volumes are displayed in Fig. ???. The expression of the multiple scattering rate α<sub>M</sub> is provided in the Methods (Eq. ??) and the single scattering rate α<sub>S</sub> is deduced from α<sub>M</sub> as fol- lows: α<sub>M</sub> = 1 − α<sub>S</sub>. The reduction of the multiple scattering rate provided by UMI is displayed in Fig. ??B.

<figure><img src="image_14.png" /><figcaption>Figure S5: **Single and multiple scattering rates.** (A) Depth evolution of the single scattering (blue) and multiple scattering (red) rates without (dashed line) and with UMI (continuous line). (B) Ratio between the multiple scattering rate measured after and before UMI. The y-scale is in dB. The ultrasound data shown here corresponds to sheep n°6, acquisition 4.</figcaption></figure>
## S6 Comparison with other aberration correction methods

The UMI framework offers the opportunity to compare our approach with other aberration correction methods. In particular, our approach shares some similarities with the iterative time-reversal (ITR) approach of speckle noise developed by Montaldo et al. (?). While matrix imaging is inspired by this work and others studies in the field of ultrasound imaging (?), it features several distinct and important advances. The first is its primary building block: the broadband focused reflection matrix that precisely selects the echoes originating from a single scattering event at each depth. This time-gating operation is decisive in terms of signal-to-noise ratio since it drastically reduces the detrimental contribution of out-of-focus and multiply scattered echoes. Second, the IPR process used in the matrix imaging process is more efficient than an iterative time reversal process. Time reversal is an adaptive filter (?) that is far from being optimal for imaging since it acts as a low pass spatial frequency filter.

This detrimental effect can be highlighted by implementing the ITR approach in the UMI pipeline. To do so, the IPR process (Eq. ??) should be replaced by ITR, such that:

$$ \sigma \mathbf{U}_{\text{in/out}}(z) = \mathbf{C}_{\text{in/out}} \times \mathbf{U}_{\text{in/out}}(z). \qquad (S3) $$

The time-reversal invariant $\mathbf{U}_{\text{in/out}}$ corresponds to the first eigenvector of the correlation matrix $\mathbf{C}_{\text{in/out}}$ and $\sigma$ is the associated eigenvalue. The time-reversal invariant yields an estimator of the input/output transmission matrix between the probe and each focal plane:

$$ \mathbf{T}_{\text{in/out}}(z) = \mathbf{U}_{\text{in/out}}(z) \circ \mathbf{T}_0 \qquad (S4) $$

The result of the ITR process is illustrated in Fig. ?? by applying it to the focused reflection matrix at depth $z = 56$ mm in sheep n°6, acquisition 4. Even if the phase of the time-reversal invariant (Fig. ??E) shares some similarities with the phase law resulting from the IPR process
(Fig. ??B), it also exhibits random fluctuations on some parts of the probe. Those areas coincide with the low amplitude displayed in the time-reversal invariant (Fig. ??D), especially on the left of the probe. This strong variation in amplitude is a manifestation of the angular filtering process induced by iterative time reversal on the received signals. When the focused reflection matrix is beamformed (Eqs. ?? and ??) using the time reversal invariant (Eq. ??), this angular filtering of ultrasound data results in enlarged focal spots, as illustrated by the resulting map of RPSFs shown in Fig. ??F. The comparison with the IPR process (Fig. ??C) therefore demonstrates the benefit of our approach compared to an ITR approach [54]. While the IPR process narrows the confocal component of each RPSF and lowers the multiple scattering background (Fig. ??C), ITR acts as an adaptive filter and enlarges each focal spot. This loss of resolution would be detrimental to ULM since it would affect the SNR and the bubble localization precision.
<figure><img src="image_15.png" /><figcaption>Figure S6: Comparison with other aberration correction methods. (A) RPSF measured at depth z = 56 mm in sheep n°6. (B) Corresponding aberration phase law extracted by an iterative phase reversal process (Eq. ??). (C) RPSFs after iterative phase reversal processing (Eqs. ??, ?? and ??). (D) Modulus of the time-reversal invariant (Eq. ??) extracted in the transducer basis. (E) Phase of the time-reversal invariant U (Eq. ??) extracted in the transducer basis. (F) RPSFs after iterative time reversal processing (Eq. ??) operated in the transducer basis. (G) Modulus of the first eigenvector V extracted by SVD of the distortion matrix expressed in the virtual source basis. (H) Phase of the first eigenvector extracted by SVD in the virtual source basis. (I) RPSFs after SVD beamforming (Eq. ??).</figcaption></figure>
The same observation holds for the SVD beamformer proposed by Bendjador et al. (?). This approach is mathematically equivalent to iterative time reversal but operated from the in-sonification basis, i.e here the virtual source basis (?). This basis and, equivalently, the plane wave basis are actually not adequate since most of aberrations are induced by the skull heterogeneities, i.e close to the probe. Again, we can show it by implementing the SVD beamformer in the UMI framework. To do so, a distortion matrix $D_{v\rho}$ can be built in the virtual source basis, following the same recipe as in the transducer basis (Eq. ??), and its first singular vector $V_{in}$ directly provides an estimator of the transmission matrix $T_{in}$, such that:

$$T_{in}(z) = V_{in}(z) \circ T_0 \qquad (S5)$$

Because aberrations are from being isoplanatic from the virtual source basis, the aberration phase law exhibited by $V_{in}$ is almost flat (Fig. ??H), resulting in an extremely limited aberration correction effect. Moreover, as for the ITR approach, the SVD provides an eigenvector $V_{in}$ that shows fluctuations in amplitude. It therefore acts as an angular filter on the corrected wave-field (Fig. ??G). As for ITR (Fig. ??F), the SVD beamformer results in enlarged focal spots (Figure ??I). The SVD beamformer would therefore be inefficient for improving the ULM process. Note that an even worse result would be obtained when the SVD beamformer is operated from the plane wave basis.
## S7 Quantitative comparison between MRA and ULM images

<figure><img src="image_16.png" /><figcaption>Figure S7: **Quantitative comparison between ULM and MRA cross-sections in the sagittal plane.** (A) MRA. (B) Raw ULM image. (C) ULM image after UMI correction. (D,E,F) Zoom on the junction between middle cerebral and posterior communicating arteries displayed as a white rectangle in panels A, B and C, respectively. For panels E and F, a structural similarity index (SI) is computed between the corresponding ULM image and the MRA image displayed in panel D.</figcaption></figure>

For a quantitative comparison between MRA and ULM images, we will consider an inset at the junction of the posterial cerebral arteries and the Willis polygon already visible on the sagittal slice displayed in Fig. ??B and D. This part of the vascular network is nicely resolved by MRA (Fig. ??D). It can therefore be used here as a ground truth for corresponding ULM images displayed in Fig. ??E and F. The effect of UMI can be quantitatively evaluated by considering the structural similarity index (SI) of each ULM image with MRA: SI=0.17 for raw ULM vs.
SI=0.33 with UMI. The SI index is therefore multiplied by a factor 2 in this specific part of the image.

ULM can do better than MRA by revealing the presence of smaller vessels (?). This higher resolving power of ULM is reinforced by UMI, both in: (i) the coronal plane with the emergence of the posterior cerebral artery in Fig. ??C that was not visible in the MRA image Fig. ??A; (ii) the sagittal plane with the resolution of smaller vessels on top of the Willis polygon in Fig. ??C that are not revealed by MRA (Fig. ??A). The latter cross-section highlights the drastic gain in contrast and resolution provided by the combination of ULM and UMI (Figs. ??C and F) compared to gold standard MRA (Figs. ??A and D) and raw ULM (Figs. ??B and E). This improvement cannot be here quantified with a structural similarity index since the MRA image is of too low quality to be used as a ground truth. However, the image quality can be quantified by its entropy *H* (?) that measures its degree of randomness. The value of *H* is reported on each image displayed in Fig. ???. While ULM only leads to a slight decrease of *H* by 0.2 compared to MRA (Figs. ??D,E), the entropy loss is increased by almost a factor 5 with UMI (Fig. ??F).
<figure><img src="image_17.png" /><figcaption>Figure S8: **Quantitative comparison between ULM and MRA cross-sections in the coronal plane.** (A) MRA. (B) ULM. (C) ULM combined with UMI. (D,E,F) Zoom on the right branch of the Willis circle depicted by the white rectangle in panels A, B and C, respectively. The entropy *H* of each image is indicated on each panel.</figcaption></figure>
S8 ULM images for each dataset

Figure ?? compares the ULM images obtained with and without UMI for different acquisitions on the three sheeps (Tab. ??). In all cases, a clear resolution and contrast improvement is observed for the different viewing angle (i.e. different probe position) in each acquisition. To better visualize such data in 3D, corresponding films in which the viewing angle changes are available in Movie S2.

The image improvement can be quantified by the loss of entropy $\Delta H$ and the contrast improvement after UMI correction. Both quantities are indicated in each image of Fig. ?? for depths ranging from 20 to 50 mm. The contrast is defined as the ratio between the variance and the mean intensity of each ULM image. Those quantities are averaged over the pixels of each ULM image for which bubbles have been detected. Zero-value pixels are discarded since they would bias the estimation of entropy and contrast.
<figure><img src="image_18.png" /><figcaption>Figure S9: **Impact of UMI on different ULM acquisitions.** ULM volumes with (left) and without (right) application of the UMI correction. Several data acquisitions on three different sheeps are presented here (see also Movie S2).</figcaption></figure>
<table><caption>Table S1: Experimental conditions for each acquisition.</caption><thead><tr><th rowspan="2">Sheep number</th><th rowspan="2">Acq. number</th><th rowspan="2">Date (day/month)</th><th colspan="2">UMI</th><th colspan="2">ULM</th><th rowspan="2">Figures</th></tr><tr><th>Time</th><th>Time</th><th>Injection volume (mL/min)</th><th>Nb of injections</th></tr></thead><tbody><tr><td>4</td><td>1</td><td>28/04</td><td>11:27</td><td>10:35</td><td>2</td><td>6</td><td>??,??,??,??</td></tr><tr><td rowspan="2">S5</td><td>1</td><td>31/05</td><td>11:21</td><td>11:00</td><td>2</td><td>3</td><td>??, ??</td></tr><tr><td>2</td><td>01/06</td><td>11:43</td><td>11:04</td><td>1.5</td><td>6</td><td>??, ??, ??,??</td></tr><tr><td rowspan="6">S6</td><td>1</td><td>17/06</td><td>10:33</td><td>10:40</td><td>1.5</td><td>5</td><td>??</td></tr><tr><td>2</td><td>17/06</td><td>10:33</td><td>10:58</td><td>2</td><td>5</td><td>??</td></tr><tr><td>3</td><td>17/06</td><td>13:33</td><td>13:11</td><td>2.5</td><td>5</td><td>??</td></tr><tr><td>4</td><td>07/07</td><td>11:49</td><td>13:00</td><td>1.75</td><td>7</td><td>??,??,??,??, ??,??, ??,??</td></tr><tr><td>5</td><td>08/07</td><td>12:36</td><td>13:33</td><td>1.5</td><td>10</td><td>??,??, ??</td></tr></tbody></table>
## Captions of Supplementary Movies

Movie S1: **Result of UMI at each depth.** Left: Selected depth. Middle Top: Local RPSFs before and after correction. Middle Bottom: Averaged RPSF before and after correction. Top right: Extracted aberration phase law. Bottom right: Resolution enhancement provided by UMI. The data shown here correspond to sheep n°6 (acquisition 4)

Movie S2: **ULM images for different acquisitions before and after UMI correction.** 3D images are shown as maximum intensity projections.