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 1 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][9] . 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 10,11 . Such a complex adaptive focusing scheme cannot be implemented in real time since it is extremely sensitive to motion 12 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 13,14 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 [15][16][17][18] . 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 19 . 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 [20][21][22][23][24] or seismic imaging 25,26 , the experimental demonstration of matrix imaging has been, so far, less spectacular with ultrasonic waves 17,18,27,28 . 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 [29][30][31][32] . Moreover, 2D imaging limits the density of independent speckle grains which controls the spatial resolution of the T-matrix estimator 28 .
In this work, we extend the ultrasound matrix imaging (UMI) framework to 3D using a fully populated matrix array of transducers [33][34][35] . 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 1,19 . Inspired by the CLASS method developed in optical microscopy 20,22 , aberrations are here compensated by a novel iterative phase reversal algorithm more efficient for 3D UMI than a singular value decomposition [16][17][18] . 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: 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 1 : 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): 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 : 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): 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 17,28 , 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 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 16,19 . The bias made on our T-matrix estimator actually scales as (see Supplementary Section 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 28 . 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) Fig. 3h].
The question that now arises is how we can, in practice, know if the convergence ofT is fulfilled without any a priori knowledge on T. An answer can be found by comparing the estimated input and output aberration phase laws,T in (u, r p ) andT out (u, r p ), at a given point r p as shown in Figs The normalized scalar product P in/out is displayed as a function of w and shows the convergence of the IPR process [ ], 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,  where much better agreement betweenT in andT out is observed for a 3D box [third panel of Fig. 3d] than for a 2D patch [right panel of Fig. 3d] of same dimension w.
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 The second difference is that our spatial reciprocity criterion P in/out is very low [see the blue box plot in Fig. 4e]. This is the manifestation of a bad convergence of our T−matrix estimator. The incoherent background exhibited by the original  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.

PSFs [
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 19 .
This filter has a Gaussian shape, with a width l c (z) that scales as 3δρ 0 (z) 19 . 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.  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 17, 28 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, 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 20-22 , seismic imaging 25,26 and also radar 55 . 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
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  To improve the signal-to-noise ratio, the R-matrix is here acquired using a set of plane waves 56 . 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 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: In the plane wave basis, τ (θ, r) is given by 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: where the symbol ⟨· · · ⟩ denotes here a spatial average over the variable in the 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 23,58,59 .
In the transducer basis (c = u), the coefficients of G ρu (z) correspond to the z−derivative of the Green's function 19 : 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 17 : 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] 17,19 : where the symbol • stands for a Hadamard product. 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: 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, 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: whereT (n) out is the estimator of T out at the n th iteration of the phase reversal process. T 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, 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 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: 2D beamforming along (y,z )-plane out , y out , z) The images displayed in Fig. 6b The contrast, F, is computed locally by decomposing the normalised RPSF as the sum of three components 28 : α 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 , 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 : where the superscript ⊤ stands for matrix transpose. In a common midpoint representation, (Eq. 28) re-writes: A local degree of anti-symmetry β is then computed as follows: 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): 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.

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.   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. 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.

Supplementary Information
This

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: (S1) This operation is described schematically in Supplementary Figure S2  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 17,19 .
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: 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 16,17,19 : 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 17 : or in terms of matrix coefficients, .

(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 17,19 , the T-matrix was estimated by performing a singular value decomposition of D rc : or, equivalently, the eigenvalue decomposition of C out : Σ 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 To circumvent that problem, one can take advantage of the analogy with iterative time reversal (ITR). The eigenvector U (1) 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: out is the estimator of T out at the n th iteration of IPR.T out is an arbitrary wave-front that initiates IPR (typically a plane wave).T out = lim n→∞T When applied to the whole field-of-view, the IPR algorithm is mathematically equivalent to the CLASS algorithm developed in optical microscopy 20 . 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 betweenT out and T out , one can expressT out as follows: By injecting Eq. S2 into the last expression,T out can be expressed, at first order, as the sum of its expected value T out and a perturbation term δT out : (S12) The bias intensity can be expressed as follows: Using Eq. S3, the numerator of the last equation can be expressed as follows: with N u the number of transducers.
The denominator of Eq. S13 can be expressed as follows: The bias intensity is thus given by: 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,16 : 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 61 . The bias intensity can thus be rewritten as: 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.

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:T 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  factor C is thus much smaller in this area, which induces a strong bias on T when wave distortions are investigated over a reduced isoplanatic patch. On the contrary, a multi-scale analysis enables a gradual enhancement of this coherence factor in this area and finally leads to an unbiased estimation of T.
Supplementary Figure S5 shows the performance of UMI by comparing the RPSFs before and after aberration compensation. In the most aberrated area (top right of the field-of-view), the resolution is improved by almost a factor two, while the contrast is increased by 4.2 dB. 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. where D(∆ρ) is a de-scanned window function that eliminates the confocal peak and W is a spatial average window function around the targeted focal point r p .
The background can be decomposed as the sum of a fully symmetric matrix associated to multiple scattering (due to spatial reciprocity) and a fully random matrix associated to the electronic noise as follows: Projecting the B−matrix onto its anti-symmetric subspace directly holds the antisymmetric part of the electronic noise such that: Assuming equi-repartition of the electronic noise onto its symmetric and antisymmetric subspace leads to: The norm of the background can be expressed as follows: 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): with β the anti-symmetric rate of the B−matrix. Wave-front of the ITR process at iteration n r m Common midpoint r p Central point of a patch ∆ρ = ρ out − ρ in Distance input/output focusing points D(∆ρ)

S12. NOTATION AND SYMBOLS
De-scanned window function