Flavien Bureau

Justine Robin

Arthur Le Ber

William Lambert

Mathias Fink

Alexandre Aubry

Three-Dimensional Ultrasound Matrix Imaging

Matrix imaging paves the way towards a next revolution in wave physics.

Based on the response matrix recorded between a set of sensors, it enables an optimized compensation of aberration phenomena and multiple scattering events that usually drastically hinder the focusing process in heterogeneous media. Although it gave rise to spectacular results in optical microscopy or seismic imaging, the success of matrix imaging has been so far relatively limited with ultrasonic waves because wave control is generally only performed with a linear array of transducers. In this paper, we extend ultrasound matrix imaging to a 3D geometry. Switching from a 1D to a 2D probe enables a much sharper estimation of the transmission matrix that links each transducer and each medium voxel. Here, we first present an experimental proof of concept on a tissue-mimicking phantom through ex-vivo tissues and then, show the potential of 3D matrix imaging for transcranial applications.

Introduction

The resolution of a wave imaging system can be defined as the ability to discern small details of an object. In conventional imaging, this resolution cannot overcome the diffraction limit of a half wavelength and may be further limited by the maximum collection angle of the imaging device. However, even with a perfect imaging system, the image quality is affected by the inhomogeneities of the propagation medium. Large-scale spatial variations of the wave velocity introduce aberrations as the wave passes through the medium of interest. Strong concentration of scatterers also induces multiple scattering events that randomize the directions of wave propagation, leading to a strong degradation of the image resolution and contrast. Such problems are encountered in all domains of wave physics, in particular for the inspection of biological tissues, whether it be by ultrasound imaging [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF] or optical microscopy 2 , or for the probing of natural resources or deep structure of the Earth's crust with seismic waves 3 .

To mitigate those problems, the concept of adaptive focusing has been adapted from astronomy where it was developed decades ago 4,5 . Ultrasound imaging employs array of transducers that allows to control and record the amplitude and phase of broadband wave-fields. Wave-front distortions can be compensated for by adjusting the time-delays added to each emitted and/or detected signal in order to focus ultrasonic waves at a certain position inside the medium [6][7][8][START_REF] Ali | Aberration correction in diagnostic ultrasound: A review of the prior field and current directions[END_REF] . The estimation of those time delays implies an iterative time-consuming focusing process that should be ideally repeated for each point in the field-of-view [START_REF] Måsøy | Iteration of transmit-beam aberration correction in medical ultrasound imaging[END_REF][START_REF] Montaldo | Time Reversal of Speckle Noise[END_REF] . Such a complex adaptive focusing scheme cannot be implemented in real time since it is extremely sensitive to motion [START_REF] Pernot | 3-D real-time motion correction in highintensity focused ultrasound therapy[END_REF] whether induced by the operator holding the probe or by the movement of tissues.

Fortunately, this tedious process can now be performed in post-processing [START_REF] Jaeger | Full correction for spatially distributed speed-of-sound in echo ultrasound based on measuring aberration delays via transmit beam steering[END_REF][START_REF] Chau | A Locally Adaptive Phase Aberration Correction (LAPAC) Method for Synthetic Aperture Sequences[END_REF] thanks to the tremendous progress made in terms of computational power and memory capacity during the last decade. To optimize the focusing process and image formation, a matrix formalism can be fruitful [START_REF] Varslot | Eigenfunction analysis of stochastic backscatter for characterization of acoustic aberration in medical ultrasound imaging[END_REF][START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF][START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Bendjador | The SVD Beamformer: Physical Principles and Application to Ultrafast Adaptive Ultrasound[END_REF] . Indeed, once the reflection matrix R of the impulse responses between each transducer is known, any physical experiment can be achieved numerically, either in a causal or anti-causal way, for any incident beam and as many times as desired. More specifically, assuming that the medium remains fixed during the acquisition, a multi-scale analysis of the wave distortions can be performed to build an estimator of the transmission matrix T between each transducer of the probe and each voxel inside the medium [START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] . Once the T-matrix is known, a local compensation of aberrations can be performed for each voxel, thereby providing a confocal image of the medium with a close to ideal resolution and an optimized contrast everywhere.

Although it gave rise to striking results in optical microscopy [START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF][START_REF] Kwon | Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin[END_REF][START_REF] Najar | Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF] or seismic imaging [START_REF] Blondel | Matrix Approach of Seismic Imaging: Application to the Erebus Volcano, Antarctica[END_REF][START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-D imaging: Application to the San Jacinto fault zone, California[END_REF] , the experimental demonstration of matrix imaging has been, so far, less spectacular with ultrasonic waves [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Bendjador | The SVD Beamformer: Physical Principles and Application to Ultrafast Adaptive Ultrasound[END_REF][START_REF] Sommer | Pixel-reassignment in ultrasound imaging[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] . Indeed, the first proof-of-concept experiments employed a linear array of transducers. Yet, aberrations in the human body are 3D-distributed and a 1D control of the wave-field is not sufficient for a fine compensation of wave-distortions as already shown by previous works [START_REF] Ivancevich | Phase-aberration correction with a 3-D ultrasound scanner: Feasibility study[END_REF][START_REF] Lacefield | Time-shift estimation and focusing through distributed aberration using multirow arrays[END_REF][START_REF] Lindsey | Pitch-catch phase aberration correction of multiple isoplanatic patches for 3-D transcranial ultrasound imaging[END_REF][START_REF] Liu | Estimation and correction of ultrasonic wavefront distortion using pulse-echo data received in a two-dimensional aperture[END_REF] . Moreover, 2D imaging limits the density of independent speckle grains which controls the spatial resolution of the T-matrix estimator [START_REF] Lambert | Ultrasound matrix imaging -Part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] .

In this work, we extend the ultrasound matrix imaging (UMI) framework to 3D using a fully populated matrix array of transducers [START_REF] Ratsimandresy | A 3 MHz two dimensional array based on piezocomposite for medical imaging[END_REF][START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF][START_REF] Provost | 3-D ultrafast Doppler imaging applied to the noninvasive mapping of blood vessels in Vivo[END_REF] . The overall method is first validated by means of a well-controlled experiment combining ex-vivo pork tissues as aberrating layer on top of a tissue-mimicking phantom. 3D UMI is then applied to a head phantom whose skull induces a strong attenuation, aberration and multiple scattering of the ultrasonic wave-field, phenomena that UMI can quantify independently of each other [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] . Inspired by the CLASS method developed in optical microscopy [START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF] , aberrations are here compensated by a novel iterative phase reversal algorithm more efficient for 3D UMI than a singular value decomposition [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF][START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Bendjador | The SVD Beamformer: Physical Principles and Application to Ultrafast Adaptive Ultrasound[END_REF] . In contrast with previous works, the convergence of this algo-rithm is ensured by investigating the spatial reciprocity between the T-matrices in transmission and reception. Throughout the paper, we will compare the gain in terms of resolution and contrast provided by 3D UMI with respect to its 2D counterpart. In particular, we will demonstrate how 3D UMI can be a powerful tool for optimizing the focusing process inside the brain through the skull.

Results

Beamforming the reflection matrix in a focused basis.

3D UMI starts with the acquisition of the reflection matrix (see Methods) by means of a 2D array of transducers (32 × 32 elements, see Fig. 1a,b). It was performed first on a tissue-mimicking phantom with nylon rods through a layer of pork tissue of fat and muscle (obtained from a chop rib piece), acting as an aberrating layer [Fig. 2a], and then on a head phantom including brain and skull-mimicking tissue, to reproduce transcranial imaging (see below). In the first experiment, the reflection matrix R uu (t) is recorded in the transducer basis [Fig. 1a,c], i.e. by acquiring the impulse responses, R(u in , u out , t), between each transducer (u) of the probe. In the head phantom experiment, skull attenuation imposes a plane wave insonification sequence [Fig. 1b] to improve the signal-to-noise ratio. The reflection matrix R θu then contains the reflected wave-field R(θ in , u out , t) recorded by the transducers u out [Fig. 1c] for each incident plane wave of angle θ in .

Whatever the illumination sequence, the reflectivity of a medium at a given point r can be estimated in post-processing by a coherent compound of incident waves delayed to virtually focus on this point, and coherently summing the echoes recorded by the probe coming from that same point [Fig. 1d]. UMI basically consists in decoupling the input (r in ) and output (r out ) focusing points [Fig. 1e]. By applying appropriate time delays to the transmission (u in /θ in ) and reception (u out ) channels (see Methods), R uu (t) and R θu (t) can be projected at each depth z in a focused basis, thereby forming a broadband focused reflection matrix, R ρρ (z) ≡ [R(ρ in , ρ out , z)].

Since the focal plane is bi-dimensional, each matrix R ρρ (z) has a four-dimension structure: R(ρ in , ρ out , z) = R({x in , y in }, {x out , y out }, z). R ρρ (z) is thus concate-nated in 2D as a set of block matrices to be represented graphically [Fig. 1g]. In such a representation, every sub-matrix of R corresponds to the reflection matrix between lines of virtual transducers located at y in and y out , whereas every element in the given sub-matrix corresponds to a specific couple (x in , x out ) [Fig. 1e].

Each coefficient R(x in , y in , x out , y out , z) corresponds to the complex amplitude of the echoes coming from the point r out = (x out , y out , z) in the focal plane when focusing at point r in = (x in , y in , z) (or conversely, since R ρρ (z) is a symmetric matrix due to spatial reciprocity).

As already shown with 2D UMI, the diagonal of R ρρ (z) directly provides the transverse cross-section of the confocal ultrasound image:

I(ρ, z) = |R(ρ in = ρ out , z)| 2 (1)
where ρ = ρ in = ρ out is the transverse coordinate of the confocal point. The corresponding 3D image is displayed in Fig. 2e for the pork tissue experiment.

Longitudinal and transverse cross-sections illustrate the effect of the aberrations induced by the pork layer by highlighting the distortion exhibited by the image of the deepest nylon rod.

Probing the focusing quality.

We now show how to quantify aberrations in ultrasound speckle (without any guide star) by investigating the antidiagonals of R ρρ (z). In the single scattering regime, the focused R-matrix coefficients can be expressed as follows [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF] :

R(ρ out , ρ in , z) = dρH out (ρ -ρ out , ρ out , z)γ(ρ, z)H in (ρ -ρ in , ρ in , z) (2)
with H in/out , the input/output point spread function (PSF); and γ the medium reflectivity. This last equation shows that each pixel of the ultrasound image (diagonal elements of R ρρ (z)) results from a convolution between the sample reflectivity and an imaging PSF, which is itself a product of the input and output PSFs. The off-diagonal points in R ρρ (z) can be exploited for a quantification of the focusing quality at any pixel of the ultrasound image by extracting each antidiagonal. Such an operation is mathematically equivalent to a change of variable to express the focused R-matrix in a common midpoint basis 1 (see Supplementary Section 2):

R M (∆ρ, r m ) = R ρ m - ∆ρ 2 , ρ m + ∆ρ 2 , z , (3)
where the subscript M stands for the common midpoint basis. r m = {ρ m , z} = {(ρ in + ρ out )/2, z} is the common midpoint between the input and output focal spots, with the two separated by a distance ∆ρ = ρ out -ρ in .

In the speckle regime (random reflectivity), this quantity probes the local focusing quality as its ensemble average intensity, which we refer to as the reflection point spread function (RPSF), scales as an incoherent convolution between the input and output PSFs 1 :

RP SF (∆ρ, r m ) = |R M (∆ρ, r m )| 2 ∝ |H in | 2 ∆ρ ⊛ |H out | 2 (∆ρ, r m ), (4)
where ⟨• • • ⟩ denotes an ensemble average, which, in practice, is performed by a local spatial average (see Methods).

Figure 1h displays the mean RPSF associated with the focused R-matrix displayed in Fig. 1g (pork tissue experiment). It clearly shows a distorted RPSF which spreads well beyond the diffraction limit (black dashed line in Fig. 1h):

δρ 0 (z) ∼ λ c 2 sin {arctan [∆u/(2z)]} (5)
with ∆u the lateral extension of the probe. The RSPF also exhibits a strong anisotropy that could not have been grasped by 2D UMI. As we will see in the next section, this kind of aberrations can only be compensated through a 3D control of the wave-field.

Adaptive focusing by iterative phase reversal.

Aberration compensation in the UMI framework is performed using the distortion matrix concept. Introduced for 2D UMI [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] , the distortion matrix can be obtained by: (i) projecting the focused R-matrix either at input or output in a correction basis (here the transducer basis, see Fig. This matrix exhibits long-range correlations that can be understood in light of isoplanicity. If in a first approximation, the pork tissue layer can be considered as a phase screen aberrator, then the input and output PSFs can be considered as

spatially invariant: H in/out (ρ -ρ in/out , r in/out ) = H(ρ -ρ in/out
). UMI consists in exploiting those correlations to determine the transfer function T (u) of the phase screen. In practice, this is done by considering the correlation matrix

C = D × D † .
The correlation between distorted wave-fields enables a virtual reflector synthesized from the set of output focal spots [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF] [Fig. 1k]. While, in previous works [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] , an iterative time-reversal process (or equivalently a singular value decomposition of D) was performed to converge towards the incident wavefront that focuses perfectly through the medium heterogeneities onto this virtual scatterer, here an iterative phase reversal algorithm is employed to build an estimator T (u) of the transfer function (see Methods). Supplementary Figure 3 demonstrates the superiority of this algorithm compared to SVD for 3D UMI.

Iterative phase reversal provides an estimation of aberration transmittance [Fig. 1k] whose phase conjugate is used to compensate for wave distortions (see Methods). The resulting mean RPSF is displayed in Fig. 1m. Although it shows a clear improvement compared with the initial RPSF, high-order aberrations still subsist. Because of its 3D feature, the pork tissue layer cannot be fully reduced to an aberrating phase screen in the transducer basis.

Spatial reciprocity as a guide star.

The 3D distribution of the speed-of-sound breaks the spatial invariance of input and output PSFs. To that aim, a local correlation matrix C(r p ) should be considered around each point r p over a sliding box W(r -r p ) (see Methods), commonly called patches, whose choice of spatial extent w is subject to the following dilemma: On the one hand, the spatial window should be as small as possible to grasp the rapid variations of the PSFs across the field of view; on the other hand, these areas should be large enough to encompass a sufficient number of independent realizations of disorder [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] . The bias made on our T-matrix estimator actually scales as (see Supplementary Section 6):

|δT (u, r p )| 2 ∼ 1 C 2 N W . (6)
C is the so-called coherence factor that is a direct indicator of the focusing quality 8

but that also depends on the multiple scattering rate and noise background [START_REF] Lambert | Ultrasound matrix imaging -Part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] . N W is the number of diffraction-limited resolution cells in each spatial window.

The validity of the T-matrix estimator in a region W 1 (Fig. 3c) is investigated by examining the corrected RPSF in a neighbour region W 2 (yellow box).

): |δT | 2 ≃ 1 -P in/out . (7)
The Multi-scale compensation of wave distortions.

The scaling of the bias intensity |δT | 2 with the coherence factor C has not been discussed yet. This dependence is however crucial since it indicates that a gradual compensation of aberrations shall be favored rather than a direct partition of the field-of-view into small boxes 22 (see Supplementary Fig. 4). An optimal UMI process should proceed as follows: first, compensate for input and output wave distortions at a large scale to increase the coherence factor C; then, decrease the spatial window W and improve the resolution of the T-matrix estimator.

The whole process can be iterated, leading to a multi-scale compensation of wave distortions (see Methods). As explained above, the convergence of the process is monitored using spatial reciprocity (P in/out >0.9).

The result of 3D UMI is displayed in Fig.

Overcoming multiple scattering for trans-cranial imaging

The same UMI process is now applied to the ultrasound data collected on the head phantom [Fig. 4a]. The parameters of the multi-scale analysis are provided in the Methods section (see also Supplementary Fig. 6). The first difference with the pork tissue experiment lies in our choice of correction basis. Given the multi-layer configuration in this experiment, the D-matrix is investigated in the plane wave basis [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF] .

The second difference is that our spatial reciprocity criterion P in/out is very low PSFs [Fig. 5c] drastically affects the coherence factor C 28 , which, in return, gives rise to a strong bias on the T-matrix estimator (Eq. 6). The incoherent background is due to multiple scattering events in the skull and electronic noise, whose relative weight can be estimated by investigating the spatial reciprocity symmetry of the R-matrix (see Methods). Fig. 5b shows the depth evolution of the single and multiple scattering contributions, as well as electronic noise. While single scattering dominates at shallow depths (z < 20 mm), multiple scattering quickly reaches 35% and remains relatively constant until electronic noise increases, so that the three contributions are almost equal at depths of 75 mm.

Beyond the depth evolution, 3D imaging even allows the study of multiple scattering in the transverse plane, as shown in Figure 5a. Two areas are examined, marked with black boxes, corresponding to the RPSFs shown in [Fig. 5c] (z = 32 mm). In the center, the RPSFs exhibits a low background due to the presence of a spherical target, resulting in a single scattering rate of 90%. The second box on the right, however, is characterized by a much higher background, leading to a multiple-to-single scattering ratio slightly larger than one. This high level of multiple scattering highlights the difficult task of trans-cranial imaging with ultrasonic waves.

In order to overcome these detrimental effects, an adaptive confocal filter can be applied to the focused R-matrix [START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] .

R ′ (ρ in , ρ out , z) = R(ρ in , ρ out , z) exp - |ρ out -ρ in | 2 2l c (z) 2 (8)
This filter has a Gaussian shape, with a width l c (z) that scales as 3δρ 0 (z) [START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] . The application of a confocal filter drastically improves the correlation between input and output aberration phase laws (see Fig. 4e and Supplementary Fig. 7), proof that a satisfying convergence towards the T-matrix is obtained. of the isoplanatic patch size: 5 mm in the transverse direction (Fig. 4f) and 2 mm in depth (Fig. 4g). This rapid variation of the aberration phase law across the field of view confirms a posteriori the necessity of a local compensation of aberrations induced by the skull. It also confirms the importance of 3D UMI with a fully sampled 2D array, as previous work recommended that the array pitch should be no more than 50% of the aberrator correlation length to properly sample the corresponding adapted focusing law [START_REF] Lacefield | Examples of design curves for multirow arrays used with time-shift compensation[END_REF] .

The phase conjugate of the T-matrix at input and output enables a fine compensation of aberrations. A set of corrected RPSFs are shown in Fig. 5d.

The comparison with their initial values demonstrates the success of 3D UMI: a diffraction-limited resolution is obtained almost everywhere [Fig. 5e)], whether it be in ultrasound speckle or in the neighborhood of bright targets, at shallow or high depths, which proves the versatility of UMI.

The performance of 3D UMI is also striking when comparing the threedimensional image of the head phantom before and after UMI. mm in Fig. 5e). Figure 6 demonstrates the necessity of a 2D ultrasonic probe for trans-cranial imaging. Indeed, the complexity of wave propagation in the skull can only be harnessed with a 3D control of the incident and reflected wave fields.

Discussion

In this experimental proof-of-concept, we demonstrated the capacity of 3D UMI to correct strong aberrations such as those encountered in trans-cranial imaging. This work is not only a 3D extension of previous studies [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] since several crucial elements have been introduced to make UMI more robust.

First, the proposed iterative phase reversal algorithm outperforms the SVD for local compensation of aberrations because it can evaluate the aberration law on a larger angular support (see Supplementary Fig. 3), resulting in a sharper compensation of aberrations. Second, the bias of our T-matrix estimator has been expressed analytically (Eq. 6) as a function of the coherence factor that grasps the detrimental effects of the virtual guide star blurring induced by aberrations, multiple scattering and noise. This led us to define a general strategy for UMI with: (i ) a multi-scale compensation of wave distortions to gradually reduce the blurring of the virtual guide star and tackle high-order aberrations associated with small isoplanatic lengths; (ii ) the application of an adaptive confocal filter to cope with multiple scattering and noise; (iii ) a fine monitoring of the convergence of our estimator by means of spatial reciprocity. The latter is a real asset, as it provides an objective criterion to check the physical significance of the extracted aberration laws and optimize the resolution of our T-matrix estimator.

Although the results presented in this paper are striking, they were obtained in vitro, and some challenges remain for in vivo brain imaging. Until now, UMI has only been applied to a static medium, while biological tissues are usually moving, especially in the case of vascular imaging, where blood flow makes the reflectivity vary quickly over time. A lot of 3D imaging modes are indeed designed to image blood flow, such as transcranial Doppler imaging 39 or ULM [START_REF] Bertolo | Whole-Brain 3D Activation and Functional Connectivity Mapping in Mice using Transcranial Functional Ultrasound Imaging[END_REF][START_REF] Chavignon | 3D Transcranial Ultrasound Localization Microscopy in the Rat Brain With a Multiplexed Matrix Probe[END_REF] . These methods are strongly sensitive to aberrations [START_REF] Demené | Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients[END_REF][START_REF] Soulioti | Super-Resolution Imaging Through the Human Skull[END_REF] and their coupling with matrix imaging would be rewarding to increase the signal-to-noise ratio and improve the image resolution, not only in the vicinity of bright reflectors 44 but also in ultrasound speckle.

However, due to spatial aliasing, the number of illuminations required for UMI scales with the number of resolution cells covered by the RPSF (see Supplementary Fig. 8). Because the aberration level through the skull is important, the illumination basis should thus be fully sampled. It limits 3D transcranial UMI to a compounded frame rate of only a few hertz, which is much too slow for ultrafast imaging [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF] . Moreover, a reduced number of illuminations breaks the symmetry of the reflection matrix. It would therefore also affect the accuracy of our monitoring parameter based on spatial reciprocity.

Soft tissues usually exhibit much slower movement, and provide signals several dB higher than blood. Ultrasound imaging of tissues is generally discarded for the brain because of the strong level of aberrations and reverberations. Interestingly, UMI can open a new route towards quantitative brain imaging since a matrix framework can also enable the mapping of physical parameters such as the speedof-sound 1,46-48 , attenuation and scattering coefficients [START_REF] Aubry | Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues[END_REF][START_REF] Brütt | Weight of single and recurrent scattering in the reflection matrix of complex media[END_REF] , or fiber anisotropy [START_REF] Papadacci | Ultrasound backscatter tensor imaging (BTI): analysis of the spatial coherence of ultrasonic speckle in anisotropic soft tissues[END_REF][START_REF] Rodriguez-Molares | Specular Beamforming[END_REF] .

Those various observables can be extremely enlightening for the characterization of cerebral tissues.

Alternatively, a solution to directly implement 3D UMI in vivo for ultrafast imaging, would be to design an imaging sequence in which the fully sampled R-matrix is acquired prior to the ultrafast acquisition itself, where the illumination basis can be drastically downsampled. The T-matrix obtained from R could then be used to correct the ultrafast images in post-processing.

Interestingly, if an ultrafast 3D UMI acquisition is possible (in cases with less aberrations, or at shallow depths), the quickly decorrelating speckle observed in blood flow can be an opportunity since it provides a large number of speckle realizations in a given voxel. A high resolution T-matrix could thus be, in principle, extracted without spatial averaging and relying on any isoplanatic assumption [START_REF] Zhao | Phase aberration correction using echo signals from moving targets i: Description and theory[END_REF][START_REF] Osmanski | Aberration correction by time reversal of moving speckle noise[END_REF] . So far, one limit of UMI concerns the strong aberration regime in which extreme time delay fluctuations can occur. Indeed, our approach relies on a broadband focused reflection matrix that consists in a coherent time gating of singly-scattered echoes. If time delay fluctuations are larger than the time resolution δt of our measurement, the angular components of each echo will not necessarily emerge in the same time gate and aberration compensation will be imperfect.

Beyond strong aberrations, another issue for transcranial imaging arises from multiple reflections caused by the skull. While such reverberations are not observed in the pork tissue experiment, their detrimental effects are much greater in a transcranial experiment because of the large impedance mismatch between the skull and brain tissues. In this work, such artefacts are not corrected and they drastically pollute the image at shallow depths (z < 20 mm).

To cope with those issues, a polychromatic approach to matrix imaging is required. Indeed, the aberration compensation scheme proposed in this paper is equivalent to a simple application of time delays on each transmit and receive channel. On the contrary, a full compensation of reverberation requires the tailoring of a complex spatio-temporal adaptive (or even inverse) filter. To that aim, 3D UMI provides an adequate framework to exploit, at best, all the spatio-temporal degrees of freedom provided by a high-dimension array of broadband transducers.

To conclude, 3D UMI is general and can be applied to any insonification sequence (plane wave or virtual source illumination) or array configuration (random or periodic, sparse or dense). Matrix imaging can be also extended to any field of wave physics for which a multi-element technology is available: optical imaging [START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF][START_REF] Badon | Distortion matrix concept for deep optical imaging in scattering media[END_REF][START_REF] Yoon | Laser scanning reflectionmatrix microscopy for aberration-free imaging through intact mouse skull[END_REF] , seismic imaging [START_REF] Blondel | Matrix Approach of Seismic Imaging: Application to the Erebus Volcano, Antarctica[END_REF][START_REF] Touma | A distortion matrix framework for high-resolution passive seismic 3-D imaging: Application to the San Jacinto fault zone, California[END_REF] and also radar [START_REF] Berland | Microwave photonic mimo radar for short-range 3d imaging[END_REF] . All the conclusions raised in that paper can be extended to each of these fields. The matrix formalism is thus a powerful tool for the big data revolution coming in wave imaging. The reflection matrix is acquired by recording the impulse response between each transducer of the probe using IQ modulation with a sampling frequency f s = 6 MHz. To that aim, each transducer u in emits successively a sinusoidal burst of three half periods at the central frequency f c . For each excitation u in , the back-scattered wave-field is recorded by all probe elements u out over a time length ∆t = 139 µs. This set of impulse responses is stored in the canonical reflection

Methods

matrix R uu (t) = [R(u in , u out , t)].
Description of the head phantom experiment.

In this second experiment, the same probe [Tab. I] is placed slightly above the temporal window of a mimicking head phantom, whose characteristics are described in Tab. II. To investigate the performance of UMI in terms of resolution and contrast, the manufacturer (True Phantom Solutions) was asked to place small spherical targets made of bone-mimicking material inside the brain. They are arranged crosswise, evenly spaced in the 3 directions with a distance of 1 cm between two consecutive targets, and their diameter increases with depth: 0.2, 0.5, 1, 2, 3 mm [Fig. 4a]. Skull thickness is of ∼ 6 mm on average at the position where the probe is placed and the first spherical target is located at z ≈ 20 mm depth, while the center of the cross is at z ≈ 40 mm depth. The transverse size of the head is ∼ 14 cm. To improve the signal-to-noise ratio, the R-matrix is here acquired using a set of plane waves [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF] . For each plane wave of angles of incidence θ in = (θ x , θ y ), the time-dependent reflected wave field R(θ in , u out , t) is recorded by each transducer u out . This set of wave-fields forms a reflection matrix acquired in the plane wave A set of 1225 plane waves are thus generated by applying appropriate time delays ∆τ (θ in , u in ) to each transducer u in = (u x , u y ) of the probe:

∆τ (θ in , u in ) = [u x sin θ x + u y sin θ y ]/c 0 . ( 9
)
Focused beamforming of the reflection matrix. The focused R-matrix,

R ρρ (z) = [R(ρ in , ρ out , z)],
is built in the time domain via a conventional delayand-sum beamforming scheme that consists in applying appropriate time-delays in order to focus at different points at input r in = (ρ in , z) = ({x in , y in }, z) and output

r out = (ρ out , z) = ({x out , y out }, z): R(ρ in , ρ out , z) = i in uout A({i in , r in }, {u out , r out })R (i in , u out , τ (i in , r in ) + τ (u out , r out )) (10)
where i = u or θ accounts for the illumination basis. A is an apodization factor that limit the extent of the synthetic aperture at emission and reception. This synthetic aperture is dictated by the transducers' directivity θ max ∼ 28 •57 .

In the transducer basis, the time-of-flights, τ (u, r), writes:

τ (u, r) = |u -r| c 0 = (x -u x ) 2 + (y -u y ) 2 + z 2 c 0 . (11)
In the plane wave basis, τ (θ, r) is given by

τ (θ, r) = x sin θ x + y sin θ y + z 1 -sin 2 θ x -sin 2 θ y /c 0 . ( 12
)
Local average of the reflection point spread function. To probe the local RPSF, the field-of-view is divided into spatial regions W(r m -r p ), defined by their center r p and their extent w = (w ρ , w z ), where w ρ and w z denote the lateral and axial extent, respectively. A local average of the back-scattered intensity can then be performed in each region:

RP SF (∆ρ, r p ) = |R M (∆ρ, r m )| 2 W(r m -r p ) rm (13)
where the symbol ⟨• • • ⟩ denotes here a spatial average over the variable in the

R ρc (z) = R ρρ (z) × G ρc (z) (14)
where the symbol × stands for the matrix product. G ρc (z) is the propagation matrix predicted by the homogeneous propagation model between the focused basis (ρ) and the correction basis (c) at each depth z. c can be either the plane wave, the transducer, or any other correction basis suitable for a particular experiment [START_REF] Kwon | Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin[END_REF][START_REF] Fink | Aberration correction in ultrasonic medical imaging with time-reversal techniques[END_REF][START_REF] Mertz | Field of view advantage of conjugate adaptive optics in microscopy applications[END_REF] .

In the transducer basis (c = u), the coefficients of G ρu (z) correspond to the z-derivative of the Green's function [START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] :

G(ρ, u, z) = ze ikc √ |u-ρ| 2 +z 2 4π(|u -ρ| 2 + z 2 ) ( 15
)
where k c is the wavenumber at the central frequency. In the Fourier basis (c = k), G ρk simply corresponds to the Fourier transform operator [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF] :

G(ρ, k) = exp (jk.ρ) = exp(j (k x x + k y y)) . (16)
At each depth z, the reflected wave-fronts contained in R ρc are then decomposed into the sum of a geometric component G ρc , that would be ideally obtained in absence of aberrations, and a distorted component that corresponds to the gap between the measured wave-fronts and their ideal counterparts [Fig. 1j] [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] :

D ρc (z) = G * ρc (z) • R ρc (z) (17)
where the symbol • stands for a Hadamard product.

D rc = D ρc (z) = [D({ρ in , z}, c out )]
is the so-called distortion matrix, here expressed at the output. Note that the same operations can be performed by exchanging input and output to obtain the input distortion matrix

D cr = [D(c in , r out )] = [D(c in , {ρ out , z})].
Local correlation analysis of the D-matrix. The next step is to exploit local correlations in D rc to extract the T-matrix. To that aim, a set of output correlation matrices C out (r p ) shall be considered between distorted wave-fronts in the vicinity of each point r p in the field-of-view:

C(c out , c ′ out , r p ) = ⟨D(r in , c out )D * (r in , c ′ out )W(r in -r p )⟩ r in (18)
An equivalent operation can be performed in input in order to extract a local correlation matrix C in (r p ) from the input distortion matrix D cr .

Iterative phase reversal algorithm. The iterative phase reversal algorithm is a computational process that provides an estimator of the transmission matrix,

T out (z) = G ⊤ ρc (z) × H out (z), (19)
where the superscript ⊤ stands for matrix transpose. T out =[T(c out , r p )] links each point c out in the dual basis and each voxel r p of the medium to be imaged [Fig. 1k].

Mathematically, the algorithm is based on the following recursive relation:

T(n) out (r p ) = exp i arg C out (r p ) × T(n-1) out (r p ) (20)
where T(n) out is the estimator of T out at the n th iteration of the phase reversal process.

T(0)

out is an arbitrary wave-front that initiates the iterative phase reversal process (typically a flat phase law) and Tout = lim n→∞ T(n) out is the result of this iterative phase reversal process.

This iterative phase reversal algorithm, repeated for each point r p , yields an estimator Tout of the T-matrix. Its digital phase conjugation enables a local compensation of aberrations [Fig. 1l]. The focused R-matrix can be updated as follows:

R (corr) ρρ (z) = D ρc (z) • T † out (z) × G † ρc (z) (21)
where the symbol † stands for transpose conjugate and • for the Hadamard product. The same process is then applied to the input correlation matrix C in for the estimation of the input transmission matrix,

T in (z) = G ⊤ ρc (z) × H in (z).
Multi-scale analysis of wave distortions. To ensure the convergence of the IPR algorithm, several iterations of the aberration correction process are performed while reducing the size of the patches W with an overlap of 50% between them.

Three correction steps are performed in the pork tissue experiment, whereas six are performed in the head phantom experiment [as described in Table III]. At each step, the correction is performed both at input and output and reciprocity between input and output aberration laws is checked. The correction process is stopped if the normalized scalar product P in/out does not reach 0.9.

Pork tissue

Head phantom Correction step 1 3 3 3 5.5 5.5 5.5 5.5 5.5 5.5 Synthesise a 1D linear array. To estimate the benefits of 3D imaging compared to 2D UMI, a simulation of a 1D array is performed on experimental ultrasound data acquired with our 2D matrix array. To that aim, cylindrical time delays are applied at input and output:

• 2 • 3 • 1 • 2 • 3 • 4 • 5 • 6 • Number of transverse patches 1 × 1 2 × 2 4 × 4 1 × 1 2 × 2 3 × 3 4 × 4 5 × 5 6 × 6 w ρ = (w x ,
τ ′ (θ (s) , s, z) = s sin θ (s) + z cos θ (s) c 0 (22)
τ ′ (u (s) , s, z) = (s -u (s) ) 2 + z 2 c 0 . ( 23
)
with s = x or y, depending on our focus plane choice.

The focused R-matrix is still built in the time domain but using this time the following delay-and-sum beamforming:

R (2D) (y in , y out , z) = θ in uout R    θ in , u out , 2D beamforming along (y,z )-plane τ ′ (θ (y) in , y in , z) + τ ′ (u (y)
out , y out , z)

+ τ ′ (θ (x) in , x f , z f ) + τ ′ (u (x) out , x f , z f ) -2z f /c 0 Cylindrical law to focus at (x f , z f )    . (24)
The images displayed in Fig. 6b,c are obtained by synthesizing input and output beams collimated in the (y, z)-plane by focusing on a line located at (x f = 0 mm, z f = 37.25 mm), thereby mimicking the beamforming process by a conventional linear array of transducers.

Estimation of contrast and resolution. Contrast and resolution are evaluated by means of the RPSF. Equivalent to the full width at half maximum commonly used in 2D UMI, the transverse resolution δρ is assessed in 3D based on the area A (-3dB) at half maximum of the RPSF amplitude:

δρ (-3dB) = A (-3dB) /π (25)
The contrast, F, is computed locally by decomposing the normalised RPSF as the sum of three components 28 :

RP SF (r p , ∆ρ) = RP SF (r p , ∆ρ) RP SF (r p , ∆ρ = 0) = α S (r p ) + α M (r p ) + α N (r p ). ( 26
)
α S is the single scattering rate that corresponds to the confocal peak. α M is a multiple scattering rate that gives rise to a diffuse halo; α N corresponds to the electronic noise rate which results in a flat plateau. A local contrast can then be deduced from the ratio between α S and the incoherent background α B = α M +α N ,

F(r p ) = α S (r p ) α B (r p ) = 1 -α B (r p ) α B (r p ) (27)
Single and multiple scattering rates. The single scattering, multiple scattering and noise rates can be directly computed from the decomposition of the RPSF (Eq. 26). However, at large depths, multiple scattering and noise are difficult to discriminate since they both give rise to a flat plateau in the RPSF. In that case, the spatial reciprocity symmetry can be invoked to differentiate their contribution.

The multiple scattering component actually gives rise to a symmetric R-matrix while electronic noise is associated with a fully random matrix. The relative part of the two components can thus be estimated by computing the degree of antisymmetry β in the R-matrix. To that aim, the R-matrix is first projected onto its anti-symmetric subspace at each depth :

R (A) ρρ (z) = R ρρ (z) -R ⊤ ρρ (z) 2 ( 28
)
where the superscript ⊤ stands for matrix transpose. In a common midpoint representation, (Eq. 28) re-writes:

R (A) M (r m , ∆ρ) = R M (r m , ∆ρ) -R M (r m , -∆ρ) 2 . ( 29
)
A local degree of anti-symmetry β is then computed as follows:

β(r p ) = R (A) M (r m , ∆ρ) 2 W(r m -r p )D(∆ρ) [rm,∆ρ] |R M (r m , ∆ρ)| 2 W(r m -r p )D(∆ρ) [rm,∆ρ] (30)
where D(∆ρ) is a de-scanned window function that eliminates the confocal peak such that the computation of β is only made by considering the incoherent background. Typically, we chose D(∆ρ) = 1 for ∆ρ > 6δρ 0 (z), and zero otherwise.

Assuming equi-partition of the electronic noise between its symmetric and antisymmetric subspace, the multiple scattering rate α M and noise ratio α N can then be deduced (see Supplementary Section 11):

α M (r p ) = (1 -2β(r p )) α B (r p ) (31)
α N (r p ) = 2β(r p )α B (r p ) (32)
In the head phantom experiment [Fig. 5b], these rates are estimated at each depth by averaging over a window of size w = (w ρ , w z ) = (20, 5.5) mm. Here, we compare the typical amount of data and computational time at each post-processing step of UMI. The comparison between 2D and 3D imaging is made using a single line of transducers versus all the transducers of our matrix array. In both cases, the pixel/voxel resolution is fixed at 0.5 mm, which corresponds approximately to one wavelength. The maximum distance between the input and output focusing points is set to 10 mm. The estimation of T is here investigated without a multi-scale analysis on a single iteration at input and output. x m = (x in + x out )/2 is the common midpoint between the input and output focal spot, with the two separated by a distance ∆x = x out -x in . These considerations can be extended to 3D imaging, so that the transverse coordinate, previously x, now becomes ρ = (x, y).

S3. CORRELATION MATRIX OF WAVE DISTORTIONS

In the accompanying paper, an iterative phase reversal (IPR) process and a multi-scale analysis of D have been implemented to retrieve the T-matrix. In the following, we provide a theoretical framework to justify this process, outline its limits and conditions of success. For sake of lighter notation, the dependence over r p will be omitted in the following.

At each step of the aberration correction process, a local correlation matrix of D is computed. The UMI process assumes the convergence of the correlation matrix C towards its ensemble average ⟨C⟩, the so-called covariance matrix [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] .

In fact, this convergence is never fully realized and C should be decomposed as the sum of this covariance matrix ⟨C⟩ and a perturbation term δC:

C = ⟨C⟩ + δC. (S2)
The intensity of the perturbation term scales as the inverse of the number N W = (w 2 ρ w z )/(δρ 2 0 δz 0 ) of resolution cells in each sub-region [START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF][START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] :

|δC(c, c ′ , r p )| 2 = |C(c, c ′ , r p )| 2 N W (S3)
This perturbation term can thus be reduced by increasing the size of the spatial window W, but at the cost of a resolution loss. In the following, we express theoretically the bias induced by this perturbation term on the estimation of Tmatrices. In particular, we will show how it scales with N W in each spatial window W and the focusing quality. To that aim, we will consider the output correlation matrix C out but a similar demonstration can be performed at input.

S4. COVARIANCE MATRIX: SYNTHESIS OF A VIRTUAL GUIDE

STAR

Under assumptions of local isoplanicity in each spatial window and random reflectivity, the covariance matrix can be expressed as follows [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF] :

⟨C out ⟩ = T out × C H × T † out , (S4)
or in terms of matrix coefficients,

⟨C(c, c ′ )⟩ = T out (c)T * out (c ′ ) dρ|H in (ρ)| 2 e -i2π (c-c ′ ).ρ λzp =C H (c,c ′ ) . ( S5
)
C H is a reference correlation matrix associated with a virtual reflector whose scattering distribution corresponds to the input focal spot intensity |H in (ρ)| 2 . This scatterer plays the role of virtual guide star in the UMI process (Fig. 1k of the accompanying paper).

S5. COMPARISON BETWEEN ITERATIVE TIME REVERSAL AND PHASE REVERSAL

In previous works on 2D UMI [START_REF] Lambert | Distortion matrix approach for ultrasound imaging of random scattering media[END_REF][START_REF] Lambert | Ultrasound matrix imaging -Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] , the T-matrix was estimated by performing a singular value decomposition of D rc :

D rc = V † in × Σ × U out , (S6)
or, equivalently, the eigenvalue decomposition of C out :

C out = U † out × Σ 2 × U out . (S7)
Σ is a diagonal matrix containing the singular values σ i in descending order: σ 1 > σ 2 > .. > σ N . U out and V in are unitary matrices that contain the orthonormal set of output and input eigenvectors, U

i) out = [U (i) out (c)] and V (i) in = [V (i) in (r)]. (
The reason of this eigenvalue decomposition can be intuitively understood by considering the asymptotic case of a point-like input focusing beam. In this ideal case, Eq. S7 becomes C out (c, c ′ ) = T out (c)T * out (c ′ ). D rc is then of rank 1 -the first output singular vector U with N u the number of transducers.

The denominator of Eq. S13 can be expressed as follows:

T † out × ⟨C out ⟩ † × ⟨C out ⟩ × T out = M 2 c T in c ⊛ T in (c) 2 (S15)
The bias intensity is thus given by:

|δT out (c)| 2 = T in c ⊛ T in (0) 2 N W c T in c ⊛ T in (c) 2 (S16)
In the last expression, we recognize the ratio between the coherent intensity (energy deposited exactly at focus) and the mean incoherent input intensity. This quantity is known as the coherence factor in ultrasound imaging 8,[START_REF] Robert | Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging[END_REF] :

C in = c T in c ⊛ T in (c) T in c ⊛ T in (0) = |H in (ρ = 0)| 2 ∆ρ -2 max dρ|H in (ρ)| 2 (S17)
In the speckle regime and for a 2D probe, the coherence factor C ranges from 0, for strong aberrations and/or multiple scattering background, to 4/9 in the ideal case [START_REF] Silverstein | Ultrasound scattering model: 2-d cross-correlation and focusing criteria-theory, simulations, and experiments[END_REF] . The bias intensity can thus be rewritten as:

|δT out (c)| 2 = 1 C 2 in N W (S18)
This last expression justifies the multi-scale analysis proposed in the accompanying paper. A gradual increase of the focusing quality, quantified by C, is required to address smaller spatial windows that scale as N W . Following this scheme, the bias made of our T-matrix estimator can be minimized. Projecting the B-matrix onto its anti-symmetric subspace directly holds the anti-symmetric part of the electronic noise such that:

B (A) = B -B ⊤ 2 = N (A) (S30)
Assuming equi-repartition of the electronic noise onto its symmetric and antisymmetric subspace leads to:

∥B (A) ∥ 2 = ∥N (A) ∥ 2 = 1 2 ∥N∥ 2 (S31)
The norm of the background can be expressed as follows:

∥B∥ 2 = ∥M∥ 2 + ∥N∥ 2 + 2 ⟨M|N⟩ ∼0 (S32)
Assuming that the scalar product between the electronic noise and the multiple scattering is zero on average, the multiple scattering rate α M can be derived by combining equations (S31) & (S32):

α M = ∥M∥ 2 ∥B∥ 2 =1 -2 ∥B (A) ∥ 2 ∥B∥ 2 β (S33)
with β the anti-symmetric rate of the B-matrix.

S12. NOTATION AND SYMBOLS

FIG. 1 .FIG. 2 .

 12 FIG. 1. 3D Ultrasound Matrix Imaging (UMI). (a,b) The R-matrix can be acquired in the transducer (a) or plane-wave (b) basis in transmit and (c) recording the back-scattered wave-field on each transducer in receive. (d) Confocal imaging consists in a simultaneous focusing of waves at input and output. (e) In UMI, the input (r in ) and output (r out ) focusing points are decoupled. (f ) x-cross-section of the (g) focused R-matrix. (h) UMI enables a quantification of aberrations by extracting a local RPSF (displayed here in amplitude) from each antidiagonal of R ρρ (z). (i) UMI then consists in a projection of the focused R-matrix in a correction (here transducer) basis at output. The resulting dual R-matrix connects each focusing point to its reflected wave-front. (j) UMI then consists in realigning those wave-fronts to isolate their distorted component from their geometrical counterpart, thereby forming the D-matrix. (k) An iterative phase reversal algorithm provides an estimator of the T-matrix between the correction basis and the mid-point of input focusing points considered in panel g. (l) The phase conjugate of the T-matrix provides a focusing law that improves the focusing process at output. (m) RPSF amplitude after the output UMI process. The ultrasound data shown in this figure corresponds to the pork tissue experiment at depth z = 40 mm.

  1i); (ii) extracting wave distortions exhibited by R when compared to a reference matrix that would have been obtained in an ideal homogeneous medium of wave velocity c 0 [Fig. 1j]. The resulting distortion matrix D = [D(u, r)] contains the aberrations induced when focusing on any point r, expressed in the correction basis.

  Figure 2b illustrates this fact by showing a map of local RPSFs (see Methods). The RPSF is more strongly distorted below the fat layer of the pork tissue (c f ≈ 1480 ± 10 m/s 36 ) than below the muscle area (c m ≈ 1560 ± 50 m/s). A full-field compensation of aberrations similar to adaptive focusing does not allow a fine compensation of aberrations [Fig. 2d1]. Access to the transmission matrix T = [T (u, r)] linking each transducer and each medium voxel is required rather than just a simple aberration transmittance T (u).

W 1 and W 2

 2 are sufficiently close to assume, in a first approximation, that they belong to the same isoplanatic patch. If the box is too small (left of Fig. 3d), our estimator has not converged yet and the correction is not valid, as shown by the degraded quality of the RPSF in W 2 [left panel of Fig. 3h] compared to its initial value[Fig. 3g]. With sufficient spatial averaging [third panel of Fig. 3d], a valid aberration law can be extracted, as shown by a corrected RPSF now close to be only diffraction-limited [third panel of Fig. 3h].The question that now arises is how we can, in practice, know if the convergence of T is fulfilled without any a priori knowledge on T. An answer can be found by comparing the estimated input and output aberration phase laws, Tin (u, r p ) and Tout (u, r p ), at a given point r p as shown in Figs.3e and f. Spatial reciprocity implies that Tin and Tout shall be equal when the convergence of the estimator is reached [third panel of Figs.3e and f]. Their normalized scalar product, P in/out = N -1 u Tin T † out , can thus be used to probe the error made on the aberration phase law |δT | 2 . Both quantities are actually related as follows (see Supplementary Section

  7

FIG. 3 .

 3 Fig.3d], a large discrepancy can be found between them. In the following, the parameter P in/out will thus be used as a guide star for monitoring the convergence of the UMI process.The scaling law of Eq. 6 with respect to N W is checked in Fig.3b. The inverse scaling of the bias with N W shows the advantage of 3D UMI over 2D UMI, since N W ∼ w d , with d the imaging dimension. This superiority is evident in Fig.3a, which shows a faster convergence with 3D boxes (green curve) than with 2D patches (orange curve). For a given precision, 3D UMI thus provides a better spatial resolution for our T-matrix estimator as shown by right panels of Figs.3f,

  where much better agreement between Tin and Tout is observed for a 3D box [third 219 panel of Fig.3d] than for a 2D patch [right panel of Fig.3d] of same dimension w.

  2. It shows the evolution of the T-matrix at each step [Fig. 2c] and the corresponding local RPSFs [Fig. 2d].In the most aberrated area (i.e. under the fat), the phase fluctuations of the aberration law corresponds to a time delay spread of 56 ns (rms). This value is comparable with past measurements through the human abdominal wall[START_REF] Hinkelman | Measurements of ultrasonic pulse arrival time and energy level variations produced by propagation through abdominal wall[END_REF] . The pork tissue layer thus induces a level of aberrations typical of standard ultrasound diagnosis. The comparison with the initial and full-field maps of RPSF highlights the benefit of a local compensation via the T-matrix, with a diffraction-limited resolution reached everywhere. The local aberration phase laws exhibited by T perfectly match with the distribution of muscle and fat in the pork tissue layer. The comparison of the final 3D image [Fig. 2f] and its cross-sections with their initial counterparts [Fig.2e] show the success of the UMI process, in particular for the deepest nylon rod, which has retrieved its straight shape. The local RPSF on the top right of Fig.2shows a contrast improvement by 4.2 dB and resolution enhancement by a factor 2 (see Methods and Supplementary Fig.5).

FIG. 4 .

 4 FIG. 4. Ultrasound Matrix Imaging (UMI) of the head phantom. (a) Top and oblique views of the experimental configuration. Image credits: Harryarts and kjpargeter on Freepik. (b,c) Original and UMI images, respectively. (d) Aberration laws at 3 different depths. From top to bottom: z = 20 mm, z = 32 mm, z = 60 mm. (e) Reciprocity criterion P in/out with or without the use of a confocal filter: Each box chart displays the median, lower and upper quartiles, and the minimum and maximum values. (f,g). Correlation function of the T-matrix in the (x, z)-plane (f ) and (x, y)plane (g), respectively. We attribute the sidelobes along the y-axis (g) to the inactive rows separating each block of 256 elements of the matrix array.

[FIG. 5 .

 5 FIG. 5. Aberrations and multiple scattering quantification. (a) Single scattering (green), multiple scattering (blue) and noise (red) rate at z = 32 mm. (b) Single scattering, multiple scattering, and noise rates as a function of depth. (c,d) Maps of local RPSFs (in amplitude) before and after correction, respectively, at three different depths (From left to right: z = 20 mm, 32 mm and 60 mm. Black boxes in panel (a) and (c) corresponds to the same area. (e) Resolution δρ (-3dB) as a function of depth. Initial resolution (red line) and its value after UMI (green line) are compared with the ideal (diffraction-limited) resolution (Eq. 5).

Figure 4d shows

  Figure 4d shows the T-matrix obtained at different depths in the brain phantom. Its spatial correlation function displayed in Figs. 4f,g provides an estimation

FIG. 6 .

 6 FIG. 6. 2D versus 3D matrix imaging in a head phantom. (a) Aberration law extracted with 2D UMI. (b,c) Original and corrected images of the same target with 2D UMI, respectively. (d) Aberration law extracted with 3D UMI for a target located at z = 38 mm. (e,f ) Original and corrected images of the same target with 3D UMI, respectively. (g) Imaging PSF before (red) and after (green) 2D (dotted line) and 3D (solid line) UMI. The depth range considered in each panel corresponds to the echo of the target located at z = 38 mm.

  Description of the pork tissue experiment. The first sample under investigation is a tissue-mimicking phantom (speed of sound: c 0 = 1540 m/s) composed of random distribution of unresolved scatterers which generate ultrasonic speckle characteristic of human tissue [Fig. 2a]. The system also contains nylon filaments placed at regular intervals, with a point-like cross-section, and, at a depth of 40 mm, a 10 mm-diameter hyperechoic cylinder, containing a higher density of unresolved scatterers. A 12-mm thick pork tissue layer is placed on top of the phantom. It is immersed in water to ensure its acoustical contact with the probe and the phantom. Since the pork layer contains a part of muscle tissue (c m ∼ 1560 m/s) and a part of fat tissue (c f ∼ 1480 m/s), it acts as an aberrating layer. This experiment mimics the situation of abdominal in vivo imaging, in which layers of fat and muscle tissues generate strong aberration and scattering at shallow depths. The acquisition of the reflection matrix is performed using a 2D matrix array of transducers (Vermon) whose characteristics are provided in Tab. I. The electronic hardware used to drive the probe was developed by Supersonic Imagine (member of Hologic group) in the context of collaboration agreement with Langevin Institute. Number of transducers 32 × 32 = 1024 (with 6 dead elements) Geometry (y-axis) 3 inactive rows between each block of 256 elements Pitch δu = 0.5 mm (≈ λ at c = 1540 m/s) Central frequency f c = 3 MHz Bandwidth (at -6dB) 80%→ ∆f = [1.8 -4.2] MHz Transducer directivity θ max = 28 • at c = 1400 m/s

  basis, R θu (t) = [R(θ in , u out , t)]. Since the transducer and plane wave bases are related by a simple Fourier transform at the central frequency, the array pitch δu and probe size ∆u dictate the angular pitch δθ and maximum angle θ max necessary to acquire a full reflection matrix in the plane wave basis such that: θ max = arcsin[λ c /(2δu)] ≈ 28 • ; δθ = arcsin [λ c /(2∆u y )] ≈ 0.8 • , with λ c = c 0 /f c the central wavelength and c 0 = 1400 m/s the speed-of-sound in the brain phantom.

  subscript. W(r m -r p ) = 1 for |ρ m -ρ p | < w ρ /2 and |z m -z p | < w z /2, and zero otherwise. The dimensions of W used for [Fig. 2b,d] are: w = (w ρ , w z ) = (3.2, 3) mm. The dimensions of W to obtain [Figs. 5c,d] are: w = (w ρ , w z ) = (4, 5.5) mm. Distortion Matrix in 3D UMI. The first step consists in projecting the focused R-matrix R ρρ (z) [Fig. 1e] onto a dual basis c at output [Fig. 1i]:

FIG. S1 .

 S1 FIG. S1. Flowchart of the UMI process.

FIG. S2 .

 S2 FIG. S2. Common midpoint representation. In 2D ultrasound imaging with linear or convex probes, the change from a (a,b) conventional to a (c,d) common midpoint representation corresponds to a 45 • rotation of the focused reflection matrix. Panels a and c show a schematic representation of such a transformation, while panels b and d show experimental ultrasound data in speckle of an ultrasound phantom with a linear probe. Note that the change from the focused to the common midpoint representation implies two new sampling grids, represented by solid and dashed lines. (e) Schematic representation of the position of the input (x in ) and output (x out ) focal spots, spaced by ∆x and their common midpoint x m . (f ) Extracted RPSF by spatial averaging over all midpoints x m at depth z = 30 mm.

( 1 )

 1 FIG. S3. Iterative Time Reversal vs. Iterative Phase Reversal. (a) The first step of ITR and IPR corresponds to the following fictitious experiment: Insonifying the medium by an arbitrary wave-front (here a plane wave) using an array of transducers and recording the reflected wave-front with the same probe. (b) The ITR process consists in time-reversing this wave-front in post-processing and sending it back into the medium, recording again the reflected wave-front, and so on. (c) The IPR process is similar but normalizes the amplitude of the time-reversed wavefront at every iteration. (d) Singular value distribution of D rc for a box W of dimension w = (w x , w y , w z ) = (2, -5, 2)mm centered around point r p = (3, -5.5, 23) mm. (e,f ) Modulus of the two first eigenvectors U (i) out . (g) Modulus of the vector C out × Tout . (h) Delimitation of muscle and fat over the probe surface. (i,j,k) Phase of U (1) out , U (2) out and Tout .

FIG. S4 .

 S4 FIG. S4. Multi-scale versus direct local analysis of wave distortions (pork chop experiment, z = 29 mm). (a) Representation of the spatial windows used at each step of UMI (see Tab. III of the accompanying paper). (b) Aberration phase laws ( T) extracted by a multi-scale analysis. (c) RPSFs after multi-scale aberration compensation. (d) Representation of the spatial windows used for a direct local compensation of wave distortions. (e) Aberration phase laws ( T) extracted by a local analysis of D. (f ) RPSFs after local aberration compensation.

TABLE I .

 I Matrix array datasheet.

TABLE II .

 II Head phantom characteristics.

		Speed-of-sound Density	Attenuation
		[m/s]	[g/cm 3 ] @2.25 MHz [dB/cm]
	Cortical bone	3000 ± 30	2.31	6.4 ± 0.3
	Trabecular bone	2800 ± 50	2.03	21 ± 2
	Brain tissue	1400 ± 10	0.99	1.0 ± 0.2
	Skin tissue	1400 ± 10	1.01	1.7 ± 0.2

TABLE III .

 III Parameters of UMI in both experiments.

  Computational insights. While the UMI process is close to real-time for 2D imaging (i.e. for linear, curve or phased array probes), 3D UMI (using a fully

	populated matrix array of transducers) is still far from it (see Tab. IV) as it in-
	volves the processing of much more ultrasound data. Even if computing a confocal
	3D image only requires a few minutes, building the focused R-matrix from the
	raw data takes a few hours (on GPU with CUDA language) while one step of
	aberration correction only lasts for a few minutes. All the post-processing was
	realized with Matlab (R2021a) on a working station with 2 processors @2.20GHz,
	128Go of RAM, and a GPU with 48 Go of dedicated memory.
			2D imaging	3D imaging
	Number of channels [Input × Output]	32 × 32 ≈ 10 3 1024 × 1024 ≈ 10 6
	Field-of-view (∆x, ∆y, ∆z)	(20, 0, 80) mm (20, 20, 80) mm
			Data Time Data	Time
	Reflection matrix acquisition: R uu (t)	6 Mo 8 ms 6 Go	260 ms
	Confocal image I(r)	53 ko 5.1 ms 2.2 Mo 1.3 min
	Matrix Imaging	Focused R-matrix: R ρρ (z) 2.2 Mo 15 ms 3.6 Go Estimation of T & correction 0.15 s	2.3 h 4.5 min

TABLE IV .

 IV Computational insights.

Acknowledgments. The authors wish to thank L. Marsac for providing initial ultrasound acquisition sequences. The authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 819261, REMI-NISCENCE project, AA).

Data availability. The ultrasound data generated in this study is available at Zenodo 60 (https://zenodo.org/record/8159177).

Code availability. Codes used to post-process the ultrasound data within this paper are available from the corresponding author upon request.

Competing interests. A.A., M.F., and W.L. are inventors of a patent related to this work held by CNRS (no. US11346819B2, published May 2022). W.L. had his PhD funded by the SuperSonic Imagine company and is now an employee of this company. All authors declare that they have no other competing interests.

Supplementary Information

This document provides further information on: (i ) the UMI workflow; (ii ) the RPSF and the common midpoint basis; (iii ) the comparison between iterative time reversal and phase reversal; (iv ) the bias of the T-matrix estimator; (v ) the comparison between a multi-scale and local analysis of wave distortions; (vi ) the impact of the confocal filter; (vii ) the effect of an incompleteness of the illumination basis.

S1. WORKFLOW

Supplementary Figure S1 shows a workflow that sums up the different steps of the UMI procedure performed in the accompanying paper.

S2. RPSF AND COMMON MIDPOINT

To probe the local focusing quality, the reflection point spread function (RPSF) can be investigated. Its extraction from the focused reflection matrix, R ρρ (z) = [R(ρ in , ρ out , z)], consists in the following change of variable to project the data into a common midpoint basis:

This operation is described schematically in Supplementary Figure S2 for the simple case of 2D imaging with a linear array of transducers. It consists in extracting each antidiagonal of the focused reflection matrix R xx (z) (red boxes in Supplementary Figure S2a), corresponding to a matrix rotation by 45

with δρ in the spatial extension of H in . The amplitude of the corresponding eigenvectors

out depends on the exact shape of the virtual guide star, that is to say, on aberrations induced by the incident wave-front.

Supplementary Figures S3e andf show the modulus of two first eigenvectors, U

out and U

(2) out . They clearly show a complementary feature. While U

out is associated with the fat layer, U

out maps onto the muscle part of the pork chop [Supplementary Figure S3h]. This result can be understood by the discontinuity of the speed-of-sound between the muscle and fat parts of the pork chop that breaks the spatial invariance and isoplanicity. As a consequence, the SVD process tends to converge onto eigenstates associated with the most isoplanatic components of D rc .

This property is not satisfactory in the present case since each eigenvector only covers a part of the probe aperture. In other words, the phases of U [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF] [Supplementary Figure S3i] and U (2) [Supplementary Figure S3j] are only satisfying estimators of T over some parts of the probe. Therefore, they cannot independently lead to an aberration compensation over the full numerical aperture.

To circumvent that problem, one can take advantage of the analogy with iterative time reversal (ITR). The eigenvector U

out can actually be seen as the result of the following fictitious experiment that consists in illuminating the virtual scatterer by an arbitrary wave-front and recording the reflected wave-field [Supplementary Figure S3a]. This wave-field is time-reversed and back-emitted towards the virtual scatterer [Supplementary Figure S3b]. This process can then be iterated many times and each step can be mathematically written as:

with W (n) , the wave-front at iteration n of the ITR process and σ, the scatterer reflectivity. ITR is shown to converge towards a time-reversal invariant that is nothing other than the first eigenvector, U

To optimize the estimation of aberrations over the full probe aperture, our idea is to modify the ITR process by still re-emitting a phase-reversed wave-field but with a constant amplitude on each probe element [Supplementary Figure S3c]. In practice, this operation is performed using the following IPR algorithm: andj].

When applied to the whole field-of-view, the IPR algorithm is mathematically equivalent to the CLASS algorithm developed in optical microscopy [START_REF] Kang | High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering[END_REF] . However, the IPR algorithm is much more efficient for a local compensation of aberrations.

For IPR, the angular resolution δθ of the aberration phase law is only limited by the angular pitch of the plane wave illumination basis or the pitch p of the transducer array in the canonical basis: δθ I ∼ λ/p. With CLASS, the resolution δθ C of the aberration law is governed by the size of the spatial window W on which the focused reflection matrix is truncated: δθ C ∼ z/w ρ . It can be particularly detrimental when high-order aberrations and small isopalanatic patches are targeted.

S6. BIAS ON THE T-MATRIX ESTIMATION

In practice, however, the T-matrix estimator is still impacted by the blurring of the synthesized guide star and the presence of diffusive background and/or noise. Therefore, the whole process shall be iterated at input and output in order to gradually refine the guide star and reduce the bias on our T-matrix estimator.

Moreover, the spatial window W over which the C-matrix is computed shall be gradually decreased in order to address the high-order aberration components, the latter one being associated with smaller isoplanatic patches.

To understand the parameters controlling the bias δT out between Tout and T out , one can express Tout as follows:

By injecting Eq. S2 into the last expression, Tout can be expressed, at first order, as the sum of its expected value T out and a perturbation term δ Tout :

The bias intensity can be expressed as follows:

Using Eq. S3, the numerator of the last equation can be expressed as follows:

S7. PROBING THE BIAS INTENSITY WITH SPATIAL RECIPROCITY

In the accompanying paper, we use the scalar product P in/out between input and output aberration phase laws to monitor the bias |δT | 2 of our T-matrix estimator.

Here we demonstrate the link between both quantities. To do so, the estimator can be written as:

with T (c, r p ) = exp [jϕ(c, r p )] and δϕ(c, r p ) the phase error of the estimator.

On the one hand, the bias intensity can be rewritten using Eq. 6 as follows:

On the other hand, the scalar product P in/out is given by

In the previous equation, the sum over c can be replaced by an ensemble average since N c = N u >> 1:

Assuming a small phase error (δϕ in/out << 1), the last equation can be rewritten as follows

Since ⟨δϕ in/out ⟩ = 0 and ⟨δϕ in δϕ out ⟩ = 0, the last expression simplifies into

Assuming an equivalent phase error at input and output (⟨|δϕ in (c, r p )| 2 ⟩ = ⟨|δϕ out (c, r p )| 2 ⟩) finally leads to:

Combining the latter expression with Eq. S20 leads to the final result:

P in/out is thus a relevant quantity to estimate the bias intensity (see Fig. 3b of the accompanying paper).

S8. MULTI-SCALE ANALYSIS OF WAVE DISTORTIONS

Supplementary Figure S4 demonstrates the benefit of a multi-scale analysis of wave distortions with a gradual decrease of spatial windows W at each step of the UMI process [Supplementary Figure S4a]. To that aim, this aberration correction scheme is compared with a direct estimation of the T-matrix over the smallest patches W [Supplementary Figure S4d]. The estimated transmission matrices T differ in both cases (see comparison between Supplementary Figures S4b ande Supplementary Figure S6 shows the evolution of the RPSF during the UMI process applied to the head phantom experiment. A gradual enhancement of the focusing process is observed at each step of UMI, which enables an estimation of the T-matrix at a higher resolution. [2] V. Ntziachristos, Going deeper than microscopy: The optical imaging frontier in biology, Nat. Methods 7, 603 (2010).

[3] 'Oz Yilmaz, Seismic Data Analysis (Society of Exploration Geophysicists, 2001).

[4] H. W. Babcock, The possibuility of compensating astronomical seeing, Publ. Astron. Soc. Pac. 65, 229 (1953).

[5] F. Roddier, ed., Adaptive Optics in Astronomy (Cambridge University Press, Cambridge, 1999).

[6] M. O'Donnell and S. Flax, Phase-aberration correction using signals from point reflectors and diffuse scatterers: measurements, IEEE Trans. Ultrason. Ferroelectr.

Basis

Symbol Adapted for

Local aberration compensation [START_REF] Lambert | Ultrasound matrix imaging -Part II: The distortion matrix for aberration correction over multiple isoplanatic patches[END_REF] Dual basis Freq. Control 35, 768 (1988).

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

[8] R. Mallart and M. Fink, Adaptive focusing in scattering media through sound?speed inhomogeneities: The van Cittert Zernike approach and focusing criterion, J.

Symbol