# Three-Dimensional Ultrasound Matrix Imaging

Flavien Bureau, Justine Robin, Arthur Le Ber, William Lambert, Mathias Fink, Alexandre Aubry

► To cite this version:

Flavien Bureau, Justine Robin, Arthur Le Ber, William Lambert, Mathias Fink, et al. 3D-HAL: Three-Dimensional Ultrasound Matrix Imaging. Nature Communications, 2023, 14, pp.6793. 10.1038/s41467-023-42338-8. http://dx.doi.org/10.1038/s41467-023-42338-8

HAL Id: hal-04027631

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

Submitted on 17 Oct 2023

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

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

<figure><img src="image_1.png" /></figure>
# Three-Dimensional Ultrasound Matrix Imaging

Flavien Bureau,¹ Justine Robin,¹ Arthur Le Ber,¹ William

Lambert,¹,² Mathias Fink,¹ and Alexandre Aubry¹

¹Institut Langevin, ESPCI Paris,

PSL University, CNRS, 75005 Paris, France

²Hologic / SuperSonic Imagine,

135 Rue Emilien Gautier, 13290 Aix-en-Provence, France

(Dated: October 17, 2023)

## Abstract

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.
## 9 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¹ or optical microscopy², or for the probing of natural resources or deep structure of the Earth’s crust with seismic waves³.

To mitigate those problems, the concept of adaptive focusing has been adapted from astronomy where it was developed decades ago⁴,⁵. 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⁶⁻⁹. 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¹⁰,¹¹. Such a complex adaptive focusing scheme cannot be implemented in real time since it is extremely sensitive to motion¹² 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¹³,¹⁴ thanks to the tremendous progress made in terms of computational power and memory capacity during the last decade. To optimize the focusing process and
37 image formation, a matrix formalism can be fruitful$^{15-18}$. Indeed, once the reflection matrix **R** of the impulse responses between each transducer is known, any 38 physical experiment can be achieved numerically, either in a causal or anti-causal 39 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 40 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 41 medium$^{19}$. Once the **T**-matrix is known, a local compensation of aberrations can 42 be performed for each voxel, thereby providing a confocal image of the medium 43 with a close to ideal resolution and an optimized contrast everywhere. 44

47 Although it gave rise to striking results in optical microscopy$^{20-24}$ or seismic 48 imaging$^{25,26}$, the experimental demonstration of matrix imaging has been, so far, 49 less spectacular with ultrasonic waves$^{17,18,27,28}$. Indeed, the first proof-of-concept 50 experiments employed a linear array of transducers. Yet, aberrations in the human 51 body are 3D-distributed and a 1D control of the wave-field is not sufficient for a fine 52 compensation of wave-distortions as already shown by previous works$^{29-32}$. More- 53 over, 2D imaging limits the density of independent speckle grains which controls 54 the spatial resolution of the **T**-matrix estimator$^{28}$.

55 In this work, we extend the ultrasound matrix imaging (UMI) framework to 56 3D using a fully populated matrix array of transducers$^{33-35}$. The overall method 57 is first validated by means of a well-controlled experiment combining ex-vivo pork 58 tissues as aberrating layer on top of a tissue-mimicking phantom. 3D UMI is 59 then applied to a head phantom whose skull induces a strong attenuation, aber- 60 ration and multiple scattering of the ultrasonic wave-field, phenomena that UMI 61 can quantify independently of each other$^{1,19}$. Inspired by the CLASS method de- 62 veloped in optical microscopy$^{20,22}$, aberrations are here compensated by a novel 63 iterative phase reversal algorithm more efficient for 3D UMI than a singular value 64 decomposition$^{16-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.

## 70 Results

### 72 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 per-
formed 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_{\theta u}$ then contains the reflected wave-field $R(\theta_{in}, u_{out}, t)$ recorded
by the transducers $u_{out}$ [Fig. 1c] for each incident plane wave of angle $\theta_{in}$.

Whatever the illumination sequence, the reflectivity of a medium at a given
point $\mathbf{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 ($\mathbf{r}_{\text{in}}$) and output ($\mathbf{r}_{\text{out}}$) focusing points [Fig. 1e]. By
applying appropriate time delays to the transmission ($\mathbf{u}_{\text{in}}/\theta_{\text{in}}$) and reception ($\mathbf{u}_{\text{out}}$)
channels (see Methods), $R_{uu}(t)$ and $R_{\theta u}(t)$ can be projected at each depth $z$ in a
focused basis, thereby forming a broadband focused reflection matrix, $R_{\rho\rho}(z) =$
<figure><img src="image_2.png" /><figcaption>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 ($\mathbf{r}_{\text{in}}$) and output ($\mathbf{r}_{\text{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 $\mathbf{R}_{\rho\rho}(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.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Ultrasound matrix imaging of a tissue-mimicking phantom through a pork tissue. (a) Schematic of the experiment. (b) Maps of original RPSFs (in amplitude) at depth $z = 29$ mm. (c) Aberration phase laws extracted at the different steps of the UMI process. (d) Corresponding RPSFs after aberration compensation at each step. (e,f) 3D confocal and UMI images with one longitudinal and transverse cross-section.</figcaption></figure>

93 $[R(\rho_{in}, \rho_{out}, z)]$.

94 Since the focal plane is bi-dimensional, each matrix $R_{\rho\rho}(z)$ has a four-dimension

95 structure: $R(\rho_{in}, \rho_{out}, z) = R\{\{x_{in}, y_{in}\}, \{x_{out}, y_{out}\}, z\}$. $R_{\rho\rho}(z)$ is thus concated.
nated in 2D as a set of block matrices to be represented graphically [Fig. 1g]. In such a representation, every sub-matrix of $\mathbf{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 $\mathbf{r}_{out} = (x_{out}, y_{out}, z)$ in the focal plane when focusing at point $\mathbf{r}_{in} = (x_{in}, y_{in}, z)$ (or conversely, since $\mathbf{R}_{\rho\rho}(z)$ is a symmetric matrix due to spatial reciprocity).

As already shown with 2D UMI, the diagonal of $\mathbf{R}_{\rho\rho}(z)$ directly provides the transverse cross-section of the confocal ultrasound image:

$$ \mathcal{I}(\rho, z) = |R(\rho_{in} = \rho_{out}, z)|^2 \quad (1) $$

where $\rho = \rho_{in} = \rho_{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 $\mathbf{R}_{\rho\rho}(z)$. In the single scattering regime, the focused $\mathbf{R}$-matrix coefficients can be expressed as follows¹:

$$ R(\rho_{out}, \rho_{in}, z) = \int d\rho H_{out}(\rho - \rho_{out}, \rho_{out}, z)\gamma(\rho, z)H_{in}(\rho - \rho_{in}, \rho_{in}, z) \quad (2) $$

with $H_{in/out}$, the input/output point spread function (PSF); and $\gamma$ the medium reflectivity. This last equation shows that each pixel of the ultrasound image (diagonal elements of $\mathbf{R}_{\rho\rho}(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
119 off-diagonal points in $R_{\rho\rho}(z)$ can be exploited for a quantification of the focusing
120 quality at any pixel of the ultrasound image by extracting each antidiagonal. Such
121 an operation is mathematically equivalent to a change of variable to express the
122 focused $R$-matrix in a common midpoint basis¹ (see Supplementary Section 2):

$$R_M(\Delta\rho, \mathbf{r}_m) = R\left(\rho_m - \frac{\Delta\rho}{2}, \rho_m + \frac{\Delta\rho}{2}, z\right), \quad (3)$$

123 where the subscript $M$ stands for the common midpoint basis. $\mathbf{r}_m = \{\rho_m, z\} = 124 \{(\rho_{in} + \rho_{out})/2, z\}$ is the common midpoint between the input and output focal
125 spots, with the two separated by a distance $\Delta\rho = \rho_{out} - \rho_{in}$.

126 In the speckle regime (random reflectivity), this quantity probes the local fo-
127 cusing quality as its ensemble average intensity, which we refer to as the reflection
128 point spread function (RPSF), scales as an incoherent convolution between the
129 input and output PSFs¹:

$$RPSF(\Delta\rho, \mathbf{r}_m) = \langle |R_M(\Delta\rho, \mathbf{r}_m)|^2 \rangle \propto |H_{in}|^2 \stackrel{\Delta\rho}{\otimes} |H_{out}|^2 (\Delta\rho, \mathbf{r}_m), \quad (4)$$

130 where $\langle \cdots \rangle$ denotes an ensemble average, which, in practice, is performed by a
131 local spatial average (see Methods).

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

$$\delta\rho_0(z) \sim \frac{\lambda_c}{2 \sin \{\arctan[\Delta u/(2z)]\}} \quad (5)$$

135 with $\Delta u$ the lateral extension of the probe. The RSPF also exhibits a strong
136 anisotropy that could not have been grasped by 2D UMI. As we will see in the
137 next section, this kind of aberrations can only be compensated through a 3D
138 control of the wave-field.
139 Adaptive focusing by iterative phase reversal.

140 Aberration compensation in the UMI framework is performed using the distor-
141 tion matrix concept. Introduced for 2D UMI<sup>17,28</sup>, the distortion matrix can be
142 obtained by: (i) projecting the focused **R**-matrix either at input or output in
143 a correction basis (here the transducer basis, see Fig. 1i); (ii) extracting wave
144 distortions exhibited by **R** when compared to a reference matrix that would have
145 been obtained in an ideal homogeneous medium of wave velocity c<sub>0</sub> [Fig. 1j]. The
146 resulting distortion matrix **D** = [D(**u**, **r**)] contains the aberrations induced when
147 focusing on any point **r**, expressed in the correction basis.

148 This matrix exhibits long-range correlations that can be understood in light of
149 isoplanicity. If in a first approximation, the pork tissue layer can be considered as
150 a phase screen aberrator, then the input and output PSFs can be considered as
151 spatially invariant: H<sub>in/out</sub>(ρ - ρ<sub>in/out</sub>, r<sub>in/out</sub>) = H(ρ - ρ<sub>in/out</sub>). UMI consists in
152 exploiting those correlations to determine the transfer function T(u) of the phase
153 screen. In practice, this is done by considering the correlation matrix C = D × D†.
154 The correlation between distorted wave-fields enables a virtual reflector synthesized
155 from the set of output focal spots<sup>17</sup> [Fig. 1k]. While, in previous works<sup>17,19</sup>, an
156 iterative time-reversal process (or equivalently a singular value decomposition of
157 D) was performed to converge towards the incident wavefront that focuses perfectly
158 through the medium heterogeneities onto this virtual scatterer, here an iterative
159 phase reversal algorithm is employed to build an estimator T̂(u) of the transfer
160 function (see Methods). Supplementary Figure 3 demonstrates the superiority of
161 this algorithm compared to SVD for 3D UMI.

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

168 Spatial reciprocity as a guide star.

169 The 3D distribution of the speed-of-sound breaks the spatial invariance of input 170 and output PSFs. Figure 2b illustrates this fact by showing a map of local RPSFs 171 (see Methods). The RPSF is more strongly distorted below the fat layer of the 172 pork tissue ($c_f \approx 1480 \pm 10 \text{ m/s}^{36}$) than below the muscle area ($c_m \approx 1560 \pm 50 \text{ m/s}$). A full-field compensation of aberrations similar to adaptive focusing does 173 not allow a fine compensation of aberrations [Fig. 2d1]. Access to the transmission 174 matrix $\mathbf{T} = [T(\mathbf{u}, \mathbf{r})]$ linking each transducer and each medium voxel is required 175 rather than just a simple aberration transmittance $T(\mathbf{u})$. 176

177 To that aim, a local correlation matrix $\mathbf{C}(\mathbf{r}_p)$ should be considered around each 178 point $\mathbf{r}_p$ over a sliding box $\mathcal{W}(\mathbf{r} - \mathbf{r}_p)$ (see Methods), commonly called patches, 179 whose choice of spatial extent $w$ is subject to the following dilemma: On the 180 one hand, the spatial window should be as small as possible to grasp the rapid 181 variations of the PSFs across the field of view; on the other hand, these areas should 182 be large enough to encompass a sufficient number of independent realizations of 183 disorder$^{16,19}$. The bias made on our $\mathbf{T}$-matrix estimator actually scales as (see 184 Supplementary Section 6):

$$ |\delta T(\mathbf{u}, \mathbf{r}_p)|^2 \sim \frac{1}{C^2 N_W} \qquad (6) $$

185 $C$ is the so-called coherence factor that is a direct indicator of the focusing quality$^8$ 186 but that also depends on the multiple scattering rate and noise background$^{28}$. $N_W$ 187 is the number of diffraction-limited resolution cells in each spatial window.

188 The validity of the $\mathbf{T}$-matrix estimator in a region $\mathcal{W}_1$ (Fig. 3c) is investigated 189 by examining the corrected RPSF in a neighbour region $\mathcal{W}_2$ (yellow box). $\mathcal{W}_1$ and 190 $\mathcal{W}_2$ are sufficiently close to assume, in a first approximation, that they belong to the 191 same isoplanatic patch. If the box is too small (left of Fig. 3d), our estimator has
192 not converged yet and the correction is not valid, as shown by the degraded quality
193 of the RPSF in $W_2$ [left panel of Fig. 3h] compared to its initial value [Fig. 3g]. With
194 sufficient spatial averaging [third panel of Fig. 3d], a valid aberration law can be
195 extracted, as shown by a corrected RPSF now close to be only diffraction-limited
196 [third panel of Fig. 3h].

197 The question that now arises is how we can, in practice, know if the convergence
198 of $\hat{T}$ is fulfilled without any a priori knowledge on $T$. An answer can be found
199 by comparing the estimated input and output aberration phase laws, $\hat{T}_{\text{in}}(\mathbf{u}, \mathbf{r}_p)$
200 and $\hat{T}_{\text{out}}(\mathbf{u}, \mathbf{r}_p)$, at a given point $\mathbf{r}_p$ as shown in Figs. 3e and f. Spatial reciprocity
201 implies that $\hat{T}_{\text{in}}$ and $\hat{T}_{\text{out}}$ shall be equal when the convergence of the estimator is
202 reached [third panel of Figs. 3e and f]. Their normalized scalar product, $P_{\text{in/out}} =$
203 $N_u^{-1}\hat{T}_{\text{in}}\hat{T}_{\text{out}}^{\dagger}$, can thus be used to probe the error made on the aberration phase law
204 $|\delta T|^2$. Both quantities are actually related as follows (see Supplementary Section
205 7):

$$|\delta T|^2 \approx 1 - P_{\text{in/out}}. \qquad (7)$$

206 The normalized scalar product $P_{\text{in/out}}$ is displayed as a function of $w$ and shows
207 the convergence of the IPR process [Fig. 3a]. For a sufficiently large box [third
208 panel of Fig. 3d], $\hat{T}$ is supposed to have converged towards $T$ when $\hat{T}_{\text{in}}$ and $\hat{T}_{\text{out}}$
209 are almost equal [third panel of Fig. 3e,f], while, for a small box [left panel of
210 Fig. 3d], a large discrepancy can be found between them. In the following, the
211 parameter $P_{\text{in/out}}$ will thus be used as a guide star for monitoring the convergence
212 of the UMI process.

213 The scaling law of Eq. 6 with respect to $N_W$ is checked in Fig. 3b. The in-
214 verse scaling of the bias with $N_W$ shows the advantage of 3D UMI over 2D UMI,
215 since $N_W \sim w^d$, with $d$ the imaging dimension. This superiority is evident in
216 Fig. 3a, which shows a faster convergence with 3D boxes (green curve) than with
217 2D patches (orange curve). For a given precision, 3D UMI thus provides a better
218 spatial resolution for our $T$-matrix estimator as shown by right panels of Figs. 3f,
<figure><img src="image_4.png" /><figcaption>FIG. 3. Convergence of the UMI process towards the T-matrix. (a) Normalized scalar product $P_{in/out}$ extracted at a point $\mathbf{r}_1$ (c) as a function of the size $w_1$ of the considered spatial window $W_1$ for 2D (orange) and 3D (green) imaging. (b) Corresponding bias intensity estimator, $|\delta T|^2 = 1 - P_{in/out}$, as a function of the number of resolution cells $N_W$ contained in the window $W_1$. The plot is in log-log scale and the theoretical power law (Eq. 6) is shown with a dashed black line for comparison. (c) Cross-section of the confocal volume showing the location of $W_1$ in green and $W_2$ in yellow. The green box $W_1$, centered around the point $\mathbf{r}_1 = (5, -5, 41)$ mm, denotes the region where the $\hat{T}$-matrix is extracted, while the yellow box $W_2$, of fixed size $w_2 = 2$ mm and centered around the point $\mathbf{r}_2 = (5, -5, 45)$ mm, is the area where the effect of aberration correction is investigated by means of the RPSF. (d) Spatial windows $W_1$ considered for the calculation of $C(\mathbf{r}_1)$. From left to right: Boxes of dimension $w = 0$ mm, $w = 0.75$ mm, $w = 1.25$ mm, rectangle of dimension $w = 1.25$ mm. (e,f) Corresponding input $\hat{T}_{in}$ and output $\hat{T}_{out}$ aberration laws, respectively. The scalar product $P_{in/out}$ is displayed in each sub-panel of (f). (g) Original RPSF associated with the yellow box $W_2$ before correction and (h) after correction using the corresponding $\hat{T}$-matrices displayed in panels (e) and (f).</figcaption></figure>

<footnote n="219">219 where much better agreement between $\hat{T}_{in}$ and $\hat{T}_{out}$ is observed for a 3D box [third</footnote>
220 panel of Fig. 3d] than for a 2D patch [right panel of Fig. 3d] of same dimension w.

221 Multi-scale compensation of wave distortions.

222 The scaling of the bias intensity |δT|^2 with the coherence factor C has not
223 been discussed yet. This dependence is however crucial since it indicates that a
224 gradual compensation of aberrations shall be favored rather than a direct partition
225 of the field-of-view into small boxes²² (see Supplementary Fig. 4). An optimal
226 UMI process should proceed as follows: first, compensate for input and output
227 wave distortions at a large scale to increase the coherence factor C; then, decrease
228 the spatial window W and improve the resolution of the T-matrix estimator.
229 The whole process can be iterated, leading to a multi-scale compensation of wave
230 distortions (see Methods). As explained above, the convergence of the process is
231 monitored using spatial reciprocity (P<sub>in/out</sub> >0.9).

232 The result of 3D UMI is displayed in Fig. 2. It shows the evolution of the
233 **T**-matrix at each step [Fig. 2c] and the corresponding local RPSFs [Fig. 2d].
234 In the most aberrated area (i.e. under the fat), the phase fluctuations of the
235 aberration law corresponds to a time delay spread of 56 ns (rms). This value is
236 comparable with past measurements through the human abdominal wall³⁷. The
237 pork tissue layer thus induces a level of aberrations typical of standard ultrasound
238 diagnosis. The comparison with the initial and full-field maps of RPSF highlights
239 the benefit of a local compensation via the **T**-matrix, with a diffraction-limited
240 resolution reached everywhere. The local aberration phase laws exhibited by **T̂**
241 perfectly match with the distribution of muscle and fat in the pork tissue layer.
242 The comparison of the final 3D image [Fig. 2f] and its cross-sections with their
243 initial counterparts [Fig. 2e] show the success of the UMI process, in particular
244 for the deepest nylon rod, which has retrieved its straight shape. The local RPSF
245 on the top right of Fig.2 shows a contrast improvement by 4.2 dB and resolution
246 enhancement by a factor 2 (see Methods and Supplementary Fig. 5).
<figure><img src="image_5.png" /><figcaption>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 $\hat{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.</figcaption></figure>

247 Overcoming multiple scattering for trans-cranial imaging

248 The same UMI process is now applied to the ultrasound data collected on the 249 head phantom [Fig. 4a]. The parameters of the multi-scale analysis are provided in 250 the Methods section (see also Supplementary Fig. 6). The first difference with the 251 pork tissue experiment lies in our choice of correction basis. Given the multi-layer 252 configuration in this experiment, the $D$-matrix is investigated in the plane wave 253 basis$^{17}$.

254 The second difference is that our spatial reciprocity criterion $P_{in/out}$ is very low 255 [see the blue box plot in Fig. 4e]. This is the manifestation of a bad convergence 256 of our $T$-matrix estimator. The incoherent background exhibited by the original
<figure><img src="image_6.png" /><figcaption>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 $\delta\rho_{(-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).</figcaption></figure>

257 PSFs [Fig. 5c] drastically affects the coherence factor $C^{28}$, which, in return, gives 258 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 259 relative weight can be estimated by investigating the spatial reciprocity symmetry 260 of the R-matrix (see Methods). Fig. 5b shows the depth evolution of the single 261 and multiple scattering contributions, as well as electronic noise. While single 262 scattering dominates at shallow depths ($z < 20$ mm), multiple scattering quickly 263 reaches 35% and remains relatively constant until electronic noise increases, so 264 that the three contributions are almost equal at depths of 75 mm. 265

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

275 In order to overcome these detrimental effects, an adaptive confocal filter can 276 be applied to the focused $\mathbf{R}$-matrix$^{19}$.

$$R'(\rho_{\text{in}}, \rho_{\text{out}}, z) = R(\rho_{\text{in}}, \rho_{\text{out}}, z) \exp \left( -\frac{|\rho_{\text{out}} - \rho_{\text{in}}|^2}{2l_c(z)^2} \right) \quad (8)$$

277 This filter has a Gaussian shape, with a width $l_c(z)$ that scales as $3\delta\rho_0(z)^{19}$. The 278 application of a confocal filter drastically improves the correlation between input 279 and output aberration phase laws (see Fig. 4e and Supplementary Fig. 7), proof 280 that a satisfying convergence towards the $\mathbf{T}$-matrix is obtained.

281 Figure 4d shows the $\mathbf{T}$-matrix obtained at different depths in the brain phan- 282 tom. Its spatial correlation function displayed in Figs. 4f,g provides an estimation 283 of the isoplanatic patch size: 5 mm in the transverse direction (Fig. 4f) and 2 mm 284 in depth (Fig. 4g). This rapid variation of the aberration phase law across the 285 field of view confirms a *posteriori* the necessity of a local compensation of aber- 286 rations induced by the skull. It also confirms the importance of 3D UMI with a 287 fully sampled 2D array, as previous work recommended that the array pitch should 288 be no more than 50% of the aberrator correlation length to properly sample the 289 corresponding adapted focusing law$^{38}$.

290 The phase conjugate of the $\mathbf{T}$-matrix at input and output enables a fine 291 compensation of aberrations. A set of corrected RPSFs are shown in Fig. 5d. 292 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 three-
dimensional image of the head phantom before and after UMI. [Figs. 4b and
298 c, respectively]. The different targets were initially strongly distorted by the skull,
299 and are now nicely resolved with UMI. In particular, the first target, located at
300 $z = 19$ mm and originally duplicated, has recovered its true shape. In addition,
301 two targets laterally spaced by 10 mm are observed at 42 mm depth, as expected
302 [Fig. 4a]. The image of the target observed at 54 mm depth is also drastically
303 improved in terms of contrast and resolution but is not found at the expected
304 transverse position. One potential explanation is the size of this target (2 mm
305 diameter) larger than the resolution cell. The guide star is thus far from being
306 point-like, which can induce an uncertainty on the absolute transverse position of
307 the target in the corrected image.

Finally, an isolated target can be leveraged to highlight the gain in contrast
provided by 3D UMI with respect to its 2D counterpart. To that aim, a linear
1D array is emulated from the same raw data by collimating the incident beam
in the y-direction [Fig. 6]. The ultrasound image is displayed before and after
UMI in Figs. 6b and c, respectively. The radial average of the corresponding focal
spots is displayed in Figs. 6d. Even though 2D UMI enables a diffraction-limited
resolution, the contrast gain *G* is quite moderate (*G*<sub>2D</sub> ~ 8dB) as it scales with the
number *N* of coherence grains exhibited by the 1D aberration phase law [Figs. 6a]:
*N*<sub>2D</sub> ~ 6.2. On the contrary, as expected, 3D UMI provides a strong enhancement
of the target echo (see the comparison between Figs. 6e,f and g): *G*<sub>3D</sub> ~ 18 dB.
The 2D aberration phase law actually provides a much larger number of spatial
degrees of freedom than its 1D counterpart: *N*<sub>3D</sub> ~ 63. The gain in contrast is
accompanied by a drastic increase of the transverse resolution (> 8× for *z* > 40
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

321 mm in Fig. 5e). Figure 6 demonstrates the necessity of a 2D ultrasonic probe for

322 trans-cranial imaging. Indeed, the complexity of wave propagation in the skull can

323 only be harnessed with a 3D control of the incident and reflected wave fields.
324 **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¹⁷,²⁸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³⁹ or ULM⁴⁰,⁴¹. These methods are strongly sensitive to aberrations⁴²,⁴³ 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⁴⁴ but also in ultrasound
352 speckle.

353 However, due to spatial aliasing, the number of illuminations required for UMI
354 scales with the number of resolution cells covered by the RPSF (see Supplemen-
355 tary Fig. 8). Because the aberration level through the skull is important, the
356 illumination basis should thus be fully sampled. It limits 3D transcranial UMI to
357 a compounded frame rate of only a few hertz, which is much too slow for ultrafast
358 imaging<sup>45</sup>. Moreover, a reduced number of illuminations breaks the symmetry of
359 the reflection matrix. It would therefore also affect the accuracy of our monitoring
360 parameter based on spatial reciprocity.

361 Soft tissues usually exhibit much slower movement, and provide signals several
362 dB higher than blood. Ultrasound imaging of tissues is generally discarded for the
363 brain because of the strong level of aberrations and reverberations. Interestingly,
364 UMI can open a new route towards quantitative brain imaging since a matrix
365 framework can also enable the mapping of physical parameters such as the speed-
366 of-sound<sup>1,46-48</sup>, attenuation and scattering coefficients<sup>49,50</sup>, or fiber anisotropy<sup>51,52</sup>.
367 Those various observables can be extremely enlightening for the characterization
368 of cerebral tissues.

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

374 Interestingly, if an ultrafast 3D UMI acquisition is possible (in cases with less
375 aberrations, or at shallow depths), the quickly decorrelating speckle observed in
376 blood flow can be an opportunity since it provides a large number of speckle real-
377 izations in a given voxel. A high resolution **T**-matrix could thus be, in principle,
378 extracted without spatial averaging and relying on any isoplanatic assumption<sup>53,54</sup>.

379 So far, one limit of UMI concerns the strong aberration regime in which extreme
380 time delay fluctuations can occur. Indeed, our approach relies on a broadband 381 focused reflection matrix that consists in a coherent time gating of singly-scattered 382 echoes. If time delay fluctuations are larger than the time resolution $\delta t$ of our 383 measurement, the angular components of each echo will not necessarily emerge in 384 the same time gate and aberration compensation will be imperfect.

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

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

398 To conclude, 3D UMI is general and can be applied to any insonification se- 399 quence (plane wave or virtual source illumination) or array configuration (ran- 400 dom or periodic, sparse or dense). Matrix imaging can be also extended to any 401 field of wave physics for which a multi-element technology is available: optical 402 imaging$^{20-22}$, seismic imaging$^{25,26}$ and also radar$^{55}$. All the conclusions raised in 403 that paper can be extended to each of these fields. The matrix formalism is thus 404 a powerful tool for the big data revolution coming in wave imaging.
## 405 Methods

406 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 \sim 1560$ m/s) and a part of fat tissue ($c_f \sim 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.

422

<table><caption>TABLE I. Matrix array datasheet.</caption><tr><td>Number of transducers</td><td>32 × 32 = 1024 (with 6 dead elements)</td></tr><tr><td>Geometry (y-axis)</td><td>3 inactive rows between each block of 256 elements</td></tr><tr><td>Pitch</td><td>δu = 0.5 mm (≈ λ at c = 1540 m/s)</td></tr><tr><td>Aperture</td><td>Δu = (Δu<sub>x</sub> &nbsp; Δu<sub>y</sub>) = (16 mm &nbsp; 17.5 mm)</td></tr><tr><td>Central frequency</td><td>f<sub>c</sub> = 3 MHz</td></tr><tr><td>Bandwidth (at -6dB)</td><td>80% → Δf = [1.8 - 4.2] MHz</td></tr><tr><td>Transducer directivity</td><td>θ<sub>max</sub> = 28° at c = 1400 m/s</td></tr></table>

423

424

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

425

426

427
428 burst of three half periods at the central frequency $f_c$. For each excitation $\mathbf{u}_{\text{in}}$, the
429 back-scattered wave-field is recorded by all probe elements $\mathbf{u}_{\text{out}}$ over a time length
430 $\Delta t = 139$ µs. This set of impulse responses is stored in the canonical reflection
431 matrix $\mathbf{R}_{\mathbf{u}\mathbf{u}}(t) = [R(\mathbf{u}_{\text{in}}, \mathbf{u}_{\text{out}}, t)]$.

432 **Description of the head phantom experiment.**

433 In this second experiment, the same probe [Tab. I] is placed slightly above
434 the temporal window of a mimicking head phantom, whose characteristics are
435 described in Tab. II. To investigate the performance of UMI in terms of resolution
436 and contrast, the manufacturer (True Phantom Solutions) was asked to place small
437 spherical targets made of bone-mimicking material inside the brain. They are
438 arranged crosswise, evenly spaced in the 3 directions with a distance of 1 cm
439 between two consecutive targets, and their diameter increases with depth: 0.2,
440 0.5, 1, 2, 3 mm [Fig. 4a]. Skull thickness is of ~ 6 mm on average at the position
441 where the probe is placed and the first spherical target is located at $z \approx 20$ mm
442 depth, while the center of the cross is at $z \approx 40$ mm depth. The transverse size of
443 the head is ~ 14 cm.

<table><caption>TABLE II. **Head phantom characteristics.**</caption><thead><tr><th></th><th>Speed-of-sound<br>[m/s]</th><th>Density<br>[g/cm<sup>3</sup>]</th><th>Attenuation<br>@2.25 MHz [dB/cm]</th></tr></thead><tbody><tr><td>Cortical bone</td><td>3000 ± 30</td><td>2.31</td><td>6.4 ± 0.3</td></tr><tr><td>Trabecular bone</td><td>2800 ± 50</td><td>2.03</td><td>21 ± 2</td></tr><tr><td>Brain tissue</td><td>1400 ± 10</td><td>0.99</td><td>1.0 ± 0.2</td></tr><tr><td>Skin tissue</td><td>1400 ± 10</td><td>1.01</td><td>1.7 ± 0.2</td></tr></tbody></table>

444 To improve the signal-to-noise ratio, the **R**-matrix is here acquired using a set
445 of plane waves⁵⁶. For each plane wave of angles of incidence $\theta_{\text{in}} = (\theta_x, \theta_y)$, the
446 time-dependent reflected wave field $R(\theta_{\text{in}}, \mathbf{u}_{\text{out}}, t)$ is recorded by each transducer
447 $\mathbf{u}_{\text{out}}$. This set of wave-fields forms a reflection matrix acquired in the plane wave
448 basis, $\mathbf{R}_{\theta\mathbf{u}}(t) = [R(\theta_{\text{in}}, \mathbf{u}_{\text{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. 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:

$$ \Delta\tau(\theta_{in}, \mathbf{u}_{in}) = [u_x \sin\theta_x + u_y \sin\theta_y]/c_0. \quad (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 delay-and-sum beamforming scheme that consists in applying appropriate time-delays in order to focus at different points at input $\mathbf{r}_{in} = (ρ_{in}, z) = (\{x_{in}, y_{in}\}, z)$ and output $\mathbf{r}_{out} = (ρ_{out}, z) = (\{x_{out}, y_{out}\}, z)$:

$$ R(ρ_{in}, ρ_{out}, z) = \sum_{\mathbf{i}_{in}} \sum_{\mathbf{u}_{out}} A(\{\mathbf{i}_{in}, \mathbf{r}_{in}\}, \{\mathbf{u}_{out}, \mathbf{r}_{out}\}) R(\mathbf{i}_{in}, \mathbf{u}_{out}, τ(\mathbf{i}_{in}, \mathbf{r}_{in}) + τ(\mathbf{u}_{out}, \mathbf{r}_{out})) \quad (10) $$

where $\mathbf{i} = \mathbf{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} \sim 28^{\circ}57$.

In the transducer basis, the time-of-flights, $τ(\mathbf{u}, \mathbf{r})$, writes:

$$ τ(\mathbf{u}, \mathbf{r}) = \frac{|\mathbf{u} - \mathbf{r}|}{c_0} = \frac{\sqrt{(x - u_x)^2 + (y - u_y)^2 + z^2}}{c_0}. \quad (11) $$
In the plane wave basis, $\tau(\theta, \mathbf{r})$ is given by

$$ \tau(\theta, \mathbf{r}) = \left[ x \sin \theta_x + y \sin \theta_y + z \sqrt{1 - \sin^2 \theta_x - \sin^2 \theta_y} \right] / c_0. \quad (12) $$

**Local average of the reflection point spread function.** To probe the local RPSF, the field-of-view is divided into spatial regions $W(\mathbf{r}_m - \mathbf{r}_p)$, defined by their center $\mathbf{r}_p$ and their extent $\mathbf{w} = (w_\rho, w_z)$, where $w_\rho$ 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:

$$ RPSF(\Delta\rho, \mathbf{r}_p) = \langle |R_M(\Delta\rho, \mathbf{r}_m)|^2 W(\mathbf{r}_m - \mathbf{r}_p) \rangle_{\mathbf{r}_m} \quad (13) $$

where the symbol $\langle \cdots \rangle$ denotes here a spatial average over the variable in the subscript. $W(\mathbf{r}_m - \mathbf{r}_p) = 1$ for $|\rho_m - \rho_p| < w_\rho/2$ and $|z_m - z_p| < w_z/2$, and zero otherwise. The dimensions of $W$ used for [Fig. 2b,d] are: $\mathbf{w} = (w_\rho, w_z) = (3.2, 3)$ mm. The dimensions of $W$ to obtain [Figs. 5c,d] are: $\mathbf{w} = (w_\rho, w_z) = (4, 5.5)$ mm.

**Distortion Matrix in 3D UMI.** The first step consists in projecting the focused $\mathbf{R}$-matrix $\mathbf{R}_{\rho\rho}(z)$ [Fig. 1e] onto a dual basis $\mathbf{c}$ at output [Fig. 1i]:

$$ \mathbf{R}_{\rho c}(z) = \mathbf{R}_{\rho \rho}(z) \times \mathbf{G}_{\rho c}(z) \quad (14) $$

where the symbol $\times$ stands for the matrix product. $\mathbf{G}_{\rho c}(z)$ is the propagation matrix predicted by the homogeneous propagation model between the focused basis ($\rho$) and the correction basis ($\mathbf{c}$) at each depth $z$. $\mathbf{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 ($\mathbf{c} = \mathbf{u}$), the coefficients of $\mathbf{G}_{\rho u}(z)$ correspond to the
478 z-derivative of the Green's function¹⁹:

$$G(\rho, \mathbf{u}, z) = \frac{ze^{ik_c \sqrt{|\mathbf{u}-\rho|^2 + z^2}}}{4\pi(|\mathbf{u}-\rho|^2 + z^2)} \quad (15)$$

479 where $k_c$ is the wavenumber at the central frequency. In the Fourier basis ($\mathbf{c} = \mathbf{k}$), 480 $\mathbf{G}_{\rho\mathbf{k}}$ simply corresponds to the Fourier transform operator¹⁷:

$$G(\rho, \mathbf{k}) = \exp(j\mathbf{k}\cdot\rho) = \exp(j(k_x x + k_y y)). \quad (16)$$

481 At each depth $z$, the reflected wave-fronts contained in $\mathbf{R}_{\rho\mathbf{c}}$ are then decomposed 482 into the sum of a geometric component $\mathbf{G}_{\rho\mathbf{c}}$, that would be ideally obtained in 483 absence of aberrations, and a distorted component that corresponds to the gap 484 between the measured wave-fronts and their ideal counterparts [Fig. 1j]¹⁷,¹⁹:

$$\mathbf{D}_{\rho\mathbf{c}}(z) = \mathbf{G}_{\rho\mathbf{c}}^*(z) \circ \mathbf{R}_{\rho\mathbf{c}}(z) \quad (17)$$

485 where the symbol $\circ$ stands for a Hadamard product. $\mathbf{D}_{\mathbf{r}\mathbf{c}} = \mathbf{D}_{\rho\mathbf{c}}(z) = [D(\{\rho_{\text{in}}, z\}, \mathbf{c}_{\text{out}})]$ 486 is the so-called distortion matrix, here expressed at the output. Note that the 487 same operations can be performed by exchanging input and output to obtain the 488 input distortion matrix $\mathbf{D}_{\mathbf{c}\mathbf{r}} = [D(\mathbf{c}_{\text{in}}, \mathbf{r}_{\text{out}})] = [D(\mathbf{c}_{\text{in}}, \{\rho_{\text{out}}, z\})]$. 489

489 Local correlation analysis of the D-matrix. The next step is to exploit 490 local correlations in $\mathbf{D}_{\mathbf{r}\mathbf{c}}$ to extract the T-matrix. To that aim, a set of output 491 correlation matrices $\mathbf{C}_{\text{out}}(\mathbf{r}_p)$ shall be considered between distorted wave-fronts in 492 the vicinity of each point $\mathbf{r}_p$ in the field-of-view:

$$C(\mathbf{c}_{\text{out}}, \mathbf{c}'_{\text{out}}, \mathbf{r}_p) = \langle D(\mathbf{r}_{\text{in}}, \mathbf{c}_{\text{out}}) D^*(\mathbf{r}_{\text{in}}, \mathbf{c}'_{\text{out}}) W(\mathbf{r}_{\text{in}} - \mathbf{r}_p) \rangle_{\mathbf{r}_{\text{in}}} \quad (18)$$

493 An equivalent operation can be performed in input in order to extract a local 494 correlation matrix $\mathbf{C}_{\text{in}}(\mathbf{r}_p)$ from the input distortion matrix $\mathbf{D}_{\mathbf{c}\mathbf{r}}$.
495 Iterative phase reversal algorithm. The iterative phase reversal algorithm is
496 a computational process that provides an estimator of the transmission matrix,

$$
\mathbf{T}_{\text{out}}(z) = \mathbf{G}_{\rho c}^{\top}(z) \times \mathbf{H}_{\text{out}}(z), \qquad (19)
$$

497 where the superscript $\top$ stands for matrix transpose. $\mathbf{T}_{\text{out}}=[\text{T}(\mathbf{c}_{\text{out}}, \mathbf{r}_{\text{p}})]$ links each
498 point $\mathbf{c}_{\text{out}}$ in the dual basis and each voxel $\mathbf{r}_{\text{p}}$ of the medium to be imaged [Fig. 1k].

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

$$
\hat{\mathbf{T}}_{\text{out}}^{(n)}(\mathbf{r}_p) = \exp \left[ i \arg \left\{ \mathbf{C}_{\text{out}}(\mathbf{r}_p) \times \hat{\mathbf{T}}_{\text{out}}^{(n-1)}(\mathbf{r}_p) \right\} \right] \quad (20)
$$

500 where $\hat{\mathbf{T}}_{\text{out}}^{(n)}$ is the estimator of $\mathbf{T}_{\text{out}}$ at the $n^{\text{th}}$ iteration of the phase reversal process.

501 $\hat{\mathbf{T}}_{\text{out}}^{(0)}$ is an arbitrary wave-front that initiates the iterative phase reversal process
502 (typically a flat phase law) and $\hat{\mathbf{T}}_{\text{out}} = \lim_{n \to \infty} \hat{\mathbf{T}}_{\text{out}}^{(n)}$ is the result of this iterative
503 phase reversal process.

504 This iterative phase reversal algorithm, repeated for each point $\mathbf{r}_p$, yields an
505 estimator $\hat{\mathbf{T}}_{\text{out}}$ of the $\mathbf{T}$-matrix. Its digital phase conjugation enables a local
506 compensation of aberrations [Fig. 1l]. The focused $\mathbf{R}$-matrix can be updated as
507 follows:

$$
\mathbf{R}_{\rho\rho}^{(\text{corr})}(z) = \left[ \mathbf{D}_{\rho c}(z) \circ \hat{\mathbf{T}}_{\text{out}}^{\dagger}(z) \right] \times \mathbf{G}_{\rho c}^{\dagger}(z) \quad (21)
$$

508 where the symbol $\dagger$ stands for transpose conjugate and $\circ$ for the Hadamard product. The same process is then applied to the input correlation matrix $\mathbf{C}_{\text{in}}$ for the
509 estimation of the input transmission matrix, $\mathbf{T}_{\text{in}}(z) = \mathbf{G}_{\rho c}^{\top}(z) \times \mathbf{H}_{\text{in}}(z)$.
510

511 **Multi-scale analysis of wave distortions.** To ensure the convergence of the
512 IPR algorithm, several iterations of the aberration correction process are performed
513 while reducing the size of the patches $\mathcal{W}$ with an overlap of 50% between them.

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

<table><caption>TABLE III. Parameters of UMI in both experiments.</caption><thead><tr><th rowspan="2">Correction step</th><th colspan="3">Pork tissue</th><th colspan="6">Head phantom</th></tr><tr><th>1°</th><th>2°</th><th>3°</th><th>1°</th><th>2°</th><th>3°</th><th>4°</th><th>5°</th><th>6°</th></tr></thead><tbody><tr><td>Number of transverse patches</td><td>1 × 1</td><td>2 × 2</td><td>4 × 4</td><td>1 × 1</td><td>2 × 2</td><td>3 × 3</td><td>4 × 4</td><td>5 × 5</td><td>6 × 6</td></tr><tr><td>w<sub>ρ</sub> = (w<sub>x</sub>, w<sub>y</sub>) [mm]</td><td>16</td><td>12</td><td>8</td><td>20</td><td>15</td><td>13.3</td><td>10</td><td>8</td><td>6.6</td></tr><tr><td>w<sub>z</sub> [mm]</td><td>3</td><td>3</td><td>3</td><td>5.5</td><td>5.5</td><td>5.5</td><td>5.5</td><td>5.5</td><td>5.5</td></tr></tbody></table>

519

520

521 **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:

$$ \tau'(\theta^{(s)}, s, z) = \frac{s \sin \theta^{(s)} + z \cos \theta^{(s)}}{c_0} \qquad (22) $$

525

$$ \tau'(u^{(s)}, s, z) = \frac{\sqrt{(s - u^{(s)})^2 + z^2}}{c_0}. \qquad (23) $$

526 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_{\text{in}}, y_{\text{out}}, z) = \sum_{\theta_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}} R \left( \underbrace{\theta_{\text{in}}, \mathbf{u}_{\text{out}}, \tau'(\theta_{\text{in}}^{(y)}, y_{\text{in}}, z) + \tau'(u_{\text{out}}^{(y)}, y_{\text{out}}, z)}_{\text{2D beamforming along }(y,z)\text{-plane}} + \underbrace{\tau'(\theta_{\text{in}}^{(x)}, x_{\text{f}}, z_{\text{f}}) + \tau'(u_{\text{out}}^{(x)}, x_{\text{f}}, z_{\text{f}}) - 2z_{\text{f}}/c_0}_{\text{Cylindrical law to focus at }(x_{\text{f}}, z_{\text{f}})} \right). \quad (24) $$

527 The images displayed in Fig. 6b,c are obtained by synthesizing input and output
528 beams collimated in the $(y, z)$-plane by focusing on a line located at $(x_f = 0 \text{ mm},$ 
529 $z_f = 37.25 \text{ mm})$, thereby mimicking the beamforming process by a conventional 
530 linear array of transducers.

531 **Estimation of contrast and resolution.** Contrast and resolution are evaluated
532 by means of the RPSF. Equivalent to the full width at half maximum commonly
533 used in 2D UMI, the transverse resolution $\delta\rho$ is assessed in 3D based on the area
534 $A_{(-3dB)}$ at half maximum of the RPSF amplitude:

$$
\delta\rho_{(-3dB)} = \sqrt{\mathcal{A}_{(-3dB)}/\pi} \qquad (25)
$$

535 The contrast, $\mathcal{F}$, is computed locally by decomposing the normalised RPSF as the
536 sum of three components$^{28}$:

$$
\overline{\text{RPSF}}(\mathbf{r}_p, \Delta\boldsymbol{\rho}) = \frac{\text{RPSF}(\mathbf{r}_p, \Delta\boldsymbol{\rho})}{\text{RPSF}(\mathbf{r}_p, \Delta\boldsymbol{\rho} = \mathbf{0})} = \alpha_S(\mathbf{r}_p) + \alpha_M(\mathbf{r}_p) + \alpha_N(\mathbf{r}_p). \quad (26)
$$

537 $\alpha_S$ is the single scattering rate that corresponds to the confocal peak. $\alpha_M$ is a
538 multiple scattering rate that gives rise to a diffuse halo; $\alpha_N$ corresponds to the
539 electronic noise rate which results in a flat plateau. A local contrast can then be
540 deduced from the ratio between $\alpha_S$ and the incoherent background $\alpha_B = \alpha_M + \alpha_N$,

$$
\mathcal{F}(\mathbf{r}_p) = \frac{\alpha_S(\mathbf{r}_p)}{\alpha_B(\mathbf{r}_p)} = \frac{1 - \alpha_B(\mathbf{r}_p)}{\alpha_B(\mathbf{r}_p)} \quad (27)
$$

541 **Single and multiple scattering rates.** The single scattering, multiple scatter-
542 ing and noise rates can be directly computed from the decomposition of the RPSF
543 (Eq. 26). However, at large depths, multiple scattering and noise are difficult to
544 discriminate since they both give rise to a flat plateau in the RPSF. In that case,
545 the spatial reciprocity symmetry can be invoked to differentiate their contribution.
546 The multiple scattering component actually gives rise to a symmetric **R**-matrix
547 while electronic noise is associated with a fully random matrix. The relative part 548 of the two components can thus be estimated by computing the degree of anti- 549 symmetry $\beta$ in the $\mathbf{R}$-matrix. To that aim, the $\mathbf{R}$-matrix is first projected onto 550 its anti-symmetric subspace at each depth :

$$ \mathbf{R}_{\rho\rho}^{(A)}(z) = \frac{\mathbf{R}_{\rho\rho}(z) - \mathbf{R}_{\rho\rho}^{\top}(z)}{2} \qquad (28) $$

551 where the superscript $\top$ stands for matrix transpose. In a common midpoint 552 representation, (Eq. 28) re-writes:

$$ R_{\mathcal{M}}^{(A)}(\mathbf{r}_m, \Delta\rho) = \frac{R_{\mathcal{M}}(\mathbf{r}_m, \Delta\rho) - R_{\mathcal{M}}(\mathbf{r}_m, -\Delta\rho)}{2} \qquad (29) $$

553 A local degree of anti-symmetry $\beta$ is then computed as follows:

$$ \beta(\mathbf{r}_p) = \frac{\left\langle \left| R_{\mathcal{M}}^{(A)}(\mathbf{r}_m, \Delta\rho) \right|^2 W(\mathbf{r}_m - \mathbf{r}_p) \mathcal{D}(\Delta\rho) \right\rangle_{[\mathbf{r}_m, \Delta\rho]}}{\left\langle |R_{\mathcal{M}}(\mathbf{r}_m, \Delta\rho)|^2 W(\mathbf{r}_m - \mathbf{r}_p) \mathcal{D}(\Delta\rho) \right\rangle_{[\mathbf{r}_m, \Delta\rho]}} \qquad (30) $$

where $\mathcal{D}(\Delta\rho)$ is a de-scanned window function that eliminates the confocal peak such that the computation of $\beta$ is only made by considering the incoherent background. Typically, we chose $\mathcal{D}(\Delta\rho) = 1$ for $\Delta\rho > 6\delta\rho_0(z)$, and zero otherwise. Assuming equi-partition of the electronic noise between its symmetric and anti-symmetric subspace, the multiple scattering rate $\alpha_M$ and noise ratio $\alpha_N$ can then be deduced (see Supplementary Section 11):

$$ \alpha_M(\mathbf{r}_p) = (1 - 2\beta(\mathbf{r}_p)) \alpha_B(\mathbf{r}_p) \qquad (31) $$

$$ \alpha_N(\mathbf{r}_p) = 2\beta(\mathbf{r}_p)\alpha_B(\mathbf{r}_p) \qquad (32) $$

554 In the head phantom experiment [Fig. 5b], these rates are estimated at each depth 555 by averaging over a window of size $\mathbf{w} = (w_\rho, w_z) = (20, 5.5)$ mm.
556 **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 involves 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.

565

<table><caption>TABLE IV. **Computational insights.** 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.</caption><thead><tr><th></th><th></th><th>2D imaging</th><th>3D imaging</th></tr></thead><tbody><tr><td>Number of channels [Input × Output]</td><td></td><td>32 × 32 ≈ 10<sup>3</sup></td><td>1024 × 1024 ≈ 10<sup>6</sup></td></tr><tr><td>Field-of-view (Δx, Δy, Δz)</td><td></td><td>(20, 0, 80) mm</td><td>(20, 20, 80) mm</td></tr><tr><td></td><td>Data</td><td>6 Mo</td><td>8 ms</td></tr><tr><td>Reflection matrix acquisition: <b>R</b><sub>uu</sub>(t)</td><td></td><td>6 Go</td><td>260 ms</td></tr><tr><td>Confocal image <i>I</i>(<b>r</b>)</td><td></td><td>53 ko</td><td>5.1 ms</td></tr><tr><td>Matrix Imaging</td><td>Focused <b>R</b>-matrix: <b>R</b><sub>ρρ</sub>(z)</td><td>2.2 Mo</td><td>2.2 Mo</td></tr><tr><td></td><td>Estimation of <b>T</b> &amp; correction</td><td>15 ms</td><td>3.6 Go</td></tr><tr><td></td><td></td><td>0.15 s</td><td>2.3 h</td></tr><tr><td></td><td></td><td></td><td>4.5 min</td></tr></tbody></table>
566 **Data availability.** The ultrasound data generated in this study is available at
567 Zenodo<sup>60</sup> (https://zenodo.org/record/8159177).

568

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

571

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

577

578 **Author Contributions.** A.A. and M.F. initiated the project. A.A. supervised
579 the project. F.B. and A.L.B. coded the ultrasound acquisition sequences. F.B.
580 and J.R. performed the experiments. F.B., A.L.B., and W.L. developed the post-
581 processing tools. F.B., J.R. and A.A. analyzed the experimental results. A.A.
582 performed the theoretical study. F.B. prepared the figures. F.B., J.R. and A.A.
583 prepared the manuscript. F.B., J.R., A.L.B., W.L., M.F., and A.A. discussed the
584 results and contributed to finalizing the manuscript.

585

586 **Competing interests.** A.A., M.F., and W.L. are inventors of a patent related
587 to this work held by CNRS (no. US11346819B2, published May 2022). W.L. had
588 his PhD funded by the SuperSonic Imagine company and is now an employee of
589 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.

## 597 S1. WORKFLOW

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

## 601 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_{\rho\rho}(z) = [R(\rho_{in}, \rho_{out}, z)]$, consists in the following change of variable to project the data into a common midpoint basis:

$$ \underbrace{\begin{bmatrix} \rho_{in} \\ \rho_{out} \\ z \end{bmatrix}}_{\text{Focused}} \rightarrow \underbrace{\begin{bmatrix} \Delta\rho \\ \rho_{m} \\ z \end{bmatrix}}_{\text{Common midpoint}} = \begin{bmatrix} \rho_{out} - \rho_{in} \\ \frac{\rho_{in} + \rho_{out}}{2} \\ z \end{bmatrix}. \qquad (S1) $$

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°. In this representation,
<figure><img src="image_8.png" /><figcaption>FIG. S1. Flowchart of the UMI process.</figcaption></figure>

610 $x_m = (x_{in} + x_{out})/2$ is the common midpoint between the input and output focal 611 spot, with the two separated by a distance $\Delta x = x_{out} - x_{in}$. These considerations 612 can be extended to 3D imaging, so that the transverse coordinate, previously $x$, 613 now becomes $\rho = (x, y)$.

614 S3. CORRELATION MATRIX OF WAVE DISTORTIONS

615 In the accompanying paper, an iterative phase reversal (IPR) process and a 616 multi-scale analysis of **D** have been implemented to retrieve the **T**-matrix. In
<figure><img src="image_9.png" /><figcaption>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 $\Delta x$ and their common midpoint $x_m$. (f) Extracted RPSF by spatial averaging over all midpoints $x_m$ at depth $z = 30$ mm.</figcaption></figure>

617 the following, we provide a theoretical framework to justify this process, outline 618 its limits and conditions of success. For sake of lighter notation, the dependence 619 over $\mathbf{r}_p$ will be omitted in the following.

620 At each step of the aberration correction process, a local correlation matrix 621 of $\mathbf{D}$ is computed. The UMI process assumes the convergence of the correlation 622 matrix $\mathbf{C}$ towards its ensemble average $\langle \mathbf{C} \rangle$, the so-called covariance matrix$^{17,19}$. 623 In fact, this convergence is never fully realized and $\mathbf{C}$ should be decomposed as
624 the sum of this covariance matrix $\langle C \rangle$ and a perturbation term $\delta C$:

$$ C = \langle C \rangle + \delta C. \qquad (S2) $$

625 The intensity of the perturbation term scales as the inverse of the number $N_W$ = 626 $(w_\rho^2 w_z)/(\delta\rho_0^2 \delta z_0)$ of resolution cells in each sub-region$^{16,17,19}$:

$$ \langle |\delta C(\mathbf{c}, \mathbf{c}', \mathbf{r}_p)|^2 \rangle = \frac{\langle |C(\mathbf{c}, \mathbf{c}', \mathbf{r}_p)|^2 \rangle}{N_W} \qquad (S3) $$

627 This perturbation term can thus be reduced by increasing the size of the spatial 628 window $W$, but at the cost of a resolution loss. In the following, we express 629 theoretically the bias induced by this perturbation term on the estimation of $T$- 630 matrices. In particular, we will show how it scales with $N_W$ in each spatial window 631 $W$ and the focusing quality. To that aim, we will consider the output correlation 632 matrix $C_{\text{out}}$ but a similar demonstration can be performed at input.

## 633 S4. COVARIANCE MATRIX: SYNTHESIS OF A VIRTUAL GUIDE 634 STAR

635 Under assumptions of local isoplanicity in each spatial window and random 636 reflectivity, the covariance matrix can be expressed as follows$^{17}$:

$$ \langle C_{\text{out}} \rangle = T_{\text{out}} \times C_H \times T_{\text{out}}^{\dagger}, \qquad (S4) $$

637 or in terms of matrix coefficients,

$$ \langle C(\mathbf{c}, \mathbf{c}') \rangle = T_{\text{out}}(\mathbf{c}) T_{\text{out}}^*(\mathbf{c}') \underbrace{\int d\boldsymbol{\rho} |H_{\text{in}}(\boldsymbol{\rho})|^2 e^{-i2\pi \frac{(\mathbf{c}-\mathbf{c}') \cdot \boldsymbol{\rho}}{\lambda z_p}}}_{=C_H(\mathbf{c},\mathbf{c}')} \qquad (S5) $$
638 $C_H$ is a reference correlation matrix associated with a virtual reflector whose scattering distribution corresponds to the input focal spot intensity $|H_{in}(\rho)|^2$. This scatterer plays the role of virtual guide star in the UMI process (Fig. 1k of the accompanying paper).

642 S5. COMPARISON BETWEEN ITERATIVE TIME REVERSAL AND PHASE REVERSAL

644 In previous works on 2D UMI$^{17,19}$, the T-matrix was estimated by performing a singular value decomposition of $D_{rc}$:

$$D_{rc} = V_{in}^{\dagger} \times \Sigma \times U_{out}, \qquad (S6)$$

646 or, equivalently, the eigenvalue decomposition of $C_{out}$:

$$C_{out} = U_{out}^{\dagger} \times \Sigma^2 \times U_{out}. \qquad (S7)$$

647 $\Sigma$ is a diagonal matrix containing the singular values $\sigma_i$ in descending order: $\sigma_1 > \sigma_2 > .. > \sigma_N$. $U_{out}$ and $V_{in}$ are unitary matrices that contain the orthonormal set of output and input eigenvectors, $U_{out}^{(i)} = [U_{out}^{(i)}(c)]$ and $V_{in}^{(i)} = [V_{in}^{(i)}(r)]$.

650 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 651 case, Eq. S7 becomes $C_{out}(c, c') = T_{out}(c)T_{out}^*(c')$. $D_{rc}$ is then of rank 1 – the first 652 output singular vector $U_{out}^{(1)}$ yields the aberration transmittance $T_{out}$. 653

654 However, in reality, the input PSF $H_{in}$ is of course far from being point-like. 655 The spectrum of $D_{rc}$ displays a continuum of singular values [Supplementary Fig- 656 ure S3d]. The effective rank of $C_{out}$ is shown to scale as the number of resolution
<figure><img src="image_10.png" /><figcaption>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 $\mathbf{w} = (w_x, w_y, w_z) = (2, -5, 2)$mm centered around point $\mathbf{r}_p = (3, -5.5, 23)$ mm. (e,f) Modulus of the two first eigenvectors $\mathbf{U}_{\text{out}}^{(i)}$. (g) Modulus of the vector $\mathbf{C}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}}$. (h) Delimitation of muscle and fat over the probe surface. (i,j,k) Phase of $\mathbf{U}_{\text{out}}^{(1)}, \mathbf{U}_{\text{out}}^{(2)}$ and $\hat{\mathbf{T}}_{\text{out}}$.</figcaption></figure>
657 cells covered by the input PSF $H_{in}^{19}$:

$$M_{\delta} \sim (\delta\rho_{in}/\delta\rho_{0})^{2}. \qquad (S8)$$

658 with $\delta\rho_{in}$ the spatial extension of $H_{in}$. The amplitude of the corresponding eigen- 659 vectors $U_{out}^{(i)}$ depends on the exact shape of the virtual guide star, that is to say, 660 on aberrations induced by the incident wave-front.

661 Supplementary Figures S3e and f show the modulus of two first eigenvectors, 662 $U_{out}^{(1)}$ and $U_{out}^{(2)}$. They clearly show a complementary feature. While $U_{out}^{(1)}$ is associ- 663 ated with the fat layer, $U_{out}^{(2)}$ maps onto the muscle part of the pork chop [Supple- 664 mentary Figure S3h]. This result can be understood by the discontinuity of the 665 speed-of-sound between the muscle and fat parts of the pork chop that breaks the 666 spatial invariance and isoplanicity. As a consequence, the SVD process tends to 667 converge onto eigenstates associated with the most isoplanatic components of $D_{rc}$.

668 This property is not satisfactory in the present case since each eigenvector only 669 covers a part of the probe aperture. In other words, the phases of $U^{(1)}$ [Supplemen- 670 tary Figure S3i] and $U^{(2)}$ [Supplementary Figure S3j] are only satisfying estimators 671 of $T$ over some parts of the probe. Therefore, they cannot independently lead to 672 an aberration compensation over the full numerical aperture.

673 To circumvent that problem, one can take advantage of the analogy with itera- 674 tive time reversal (ITR). The eigenvector $U_{out}^{(1)}$ can actually be seen as the result of 675 the following fictitious experiment that consists in illuminating the virtual scatterer 676 by an arbitrary wave-front and recording the reflected wave-field [Supplementary 677 Figure S3a]. This wave-field is time-reversed and back-emitted towards the virtual 678 scatterer [Supplementary Figure S3b]. This process can then be iterated many 679 times and each step can be mathematically written as:

$$\sigma \mathbf{W}^{(n+1)} = \mathbf{C}_{\text{out}} \times \mathbf{W}^{(n)} \qquad (S9)$$
680 with $\mathbf{W}^{(n)}$, the wave-front at iteration $n$ of the ITR process and $\sigma$, the scatterer
681 reflectivity. ITR is shown to converge towards a time-reversal invariant that is
682 nothing other than the first eigenvector, $\mathbf{U}_{\text{out}}^{(1)} = \lim_{n \to +\infty} \mathbf{W}^{(n)}$.

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

$$ \hat{\mathbf{T}}_{\text{out}}^{(n+1)} = \exp \left[ i \arg \left\{ \mathbf{C}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}}^{(n)} \right\} \right] \quad (\text{S10}) $$

687 where $\hat{\mathbf{T}}_{\text{out}}^{(n)}$ is the estimator of $\mathbf{T}_{\text{out}}$ at the $n^{\text{th}}$ iteration of IPR. $\hat{\mathbf{T}}_{\text{out}}^{(0)}$ is an arbitrary wave-front that initiates IPR (typically a plane wave). $\hat{\mathbf{T}}_{\text{out}} = \lim_{n \to +\infty} \hat{\mathbf{T}}_{\text{out}}^{(n)}$ is the result of this IPR process. Unlike ITR, IPR equally addresses each angular component of the imaging process to reach a diffraction-limited resolution. Supplementary Figure S3g illustrates this fact by showing the modulus of $\mathbf{C}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}}$.
688
689 Compared with $\mathbf{U}_{\text{out}}^{(1)}$ [Supplementary Figure S3e] and $\mathbf{U}_{\text{out}}^{(2)}$ [Supplementary Figure S3f], it clearly shows that the phase-reversed invariant $\hat{\mathbf{T}}_{\text{out}}$ simultaneously addresses each angular component of the aberrated wave-field. $\hat{\mathbf{T}}_{\text{out}}$ is thus a much better estimator of the $\mathbf{T}$-matrix [Supplementary Figure S3k] than the aberration phase laws extracted by the SVD process [Supplementary Figures S3i and j].
690
691
692
693
694
695
696

697 When applied to the whole field-of-view, the IPR algorithm is mathematically
698 equivalent to the CLASS algorithm developed in optical microscopy²⁰. However,
699 the IPR algorithm is much more efficient for a local compensation of aberrations.
700 For IPR, the angular resolution $\delta\theta$ of the aberration phase law is only limited by the
701 angular pitch of the plane wave illumination basis or the pitch $p$ of the transducer
702 array in the canonical basis: $\delta\theta_I \sim \lambda/p$. With CLASS, the resolution $\delta\theta_C$ of the
703 aberration law is governed by the size of the spatial window $\mathcal{W}$ on which the focused
704 reflection matrix is truncated: $\delta\theta_C \sim z/w_\rho$. It can be particularly detrimental
705 when high-order aberrations and small isoplanatic patches are targeted.

706 S6. BIAS ON THE T-MATRIX ESTIMATION

707 In practice, however, the T-matrix estimator is still impacted by the blurring
708 of the synthesized guide star and the presence of diffusive background and/or
709 noise. Therefore, the whole process shall be iterated at input and output in order
710 to gradually refine the guide star and reduce the bias on our T-matrix estimator.
711 Moreover, the spatial window W over which the C-matrix is computed shall be
712 gradually decreased in order to address the high-order aberration components, the
713 latter one being associated with smaller isoplanatic patches.

714 To understand the parameters controlling the bias $\delta \mathbf{T}_{\text{out}}$ between $\hat{\mathbf{T}}_{\text{out}}$ and $\mathbf{T}_{\text{out}}$,
715 one can express $\hat{\mathbf{T}}_{\text{out}}$ as follows:

$$
\hat{\mathbf{T}}_{\text{out}} = \exp \left( j \arg \left\{ \mathbf{C}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}} \right\} \right) = \frac{\mathbf{C}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}}}{\| \mathbf{C}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}} \|} \quad (\text{S11})
$$

716 By injecting Eq. S2 into the last expression, $\hat{\mathbf{T}}_{\text{out}}$ can be expressed, at first order,
717 as the sum of its expected value $\mathbf{T}_{\text{out}}$ and a perturbation term $\delta\hat{\mathbf{T}}_{\text{out}}$:

$$
\hat{\mathbf{T}}_{\text{out}} = \underbrace{\frac{\langle \mathbf{C}_{\text{out}} \rangle \times \mathbf{T}_{\text{out}}}{\| \langle \mathbf{C}_{\text{out}} \rangle \times \mathbf{T}_{\text{out}} \|}}_{=\mathbf{T}_{\text{out}}} + \underbrace{\frac{\delta \mathbf{C}_{\text{out}} \times \mathbf{T}_{\text{out}}}{\| \langle \mathbf{C}_{\text{out}} \rangle \times \mathbf{T}_{\text{out}} \|}}_{\simeq \delta \hat{\mathbf{T}}_{\text{out}}}. \quad (\text{S12})
$$

718 The bias intensity can be expressed as follows:

$$
|\delta \mathbf{T}_{\text{out}}|^2 = \frac{\mathbf{T}_{\text{out}}^{\dagger} \times \delta \mathbf{C}_{\text{out}}^{\dagger} \times \delta \mathbf{C}_{\text{out}} \times \mathbf{T}_{\text{out}}}{\mathbf{T}_{\text{out}}^{\dagger} \times \langle \mathbf{C}_{\text{out}} \rangle^{\dagger} \times \langle \mathbf{C}_{\text{out}} \rangle \times \mathbf{T}_{\text{out}}} \quad (\text{S13})
$$

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

$$
\mathbf{T}_{\text{out}}^{\dagger} \times \delta \mathbf{C}_{\text{out}}^{\dagger} \times \delta \mathbf{C}_{\text{out}} \times \mathbf{T}_{\text{out}} = N_{u}^{2} \langle |\delta \mathbf{C}(\mathbf{c}, \mathbf{c}')|^{2} \rangle = N_{u}^{2} |\mathbf{C}(\mathbf{c}, \mathbf{c}')|^{2} / N_{W}. \quad (\text{S14})
$$
720 with $N_u$ the number of transducers.

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

$$ \mathbf{T}_{\text{out}}^{\dagger} \times \langle \mathbf{C}_{\text{out}} \rangle^{\dagger} \times \langle \mathbf{C}_{\text{out}} \rangle \times \mathbf{T}_{\text{out}} = M^2 \left| \sum_{\mathbf{c}} T_{\text{in}} \stackrel{\mathbf{c}}{\otimes} T_{\text{in}}(\mathbf{c}) \right|^2 \quad (\text{S15}) $$

722 The bias intensity is thus given by:

$$ |\delta T_{\text{out}}(\mathbf{c})|^2 = \frac{\left| T_{\text{in}} \stackrel{\mathbf{c}}{\otimes} T_{\text{in}}(\mathbf{0}) \right|^2}{N_W \left| \sum_{\mathbf{c}} T_{\text{in}} \stackrel{\mathbf{c}}{\otimes} T_{\text{in}}(\mathbf{c}) \right|^2} \quad (\text{S16}) $$

723 In the last expression, we recognize the ratio between the coherent intensity (energy 724 deposited exactly at focus) and the mean incoherent input intensity. This quantity 725 is known as the coherence factor in ultrasound imaging$^{8,16}$:

$$ C_{\text{in}} = \frac{\sum_{\mathbf{c}} T_{\text{in}} \stackrel{\mathbf{c}}{\otimes} T_{\text{in}}(\mathbf{c})}{T_{\text{in}} \stackrel{\mathbf{c}}{\otimes} T_{\text{in}}(\mathbf{0})} = \frac{|H_{\text{in}}(\boldsymbol{\rho} = \mathbf{0})|^2}{\Delta\rho_{\text{max}}^{-2} \int d\boldsymbol{\rho} |H_{\text{in}}(\boldsymbol{\rho})|^2} \quad (\text{S17}) $$

726 In the speckle regime and for a 2D probe, the coherence factor $C$ ranges from 0, 727 for strong aberrations and/or multiple scattering background, to 4/9 in the ideal 728 case$^{61}$. The bias intensity can thus be rewritten as:

$$ |\delta T_{\text{out}}(\mathbf{c})|^2 = \frac{1}{C_{\text{in}}^2 N_W} \quad (\text{S18}) $$

729 This last expression justifies the multi-scale analysis proposed in the accompanying 730 paper. A gradual increase of the focusing quality, quantified by $C$, is required to 731 address smaller spatial windows that scale as $N_W$. Following this scheme, the bias 732 made of our $\mathbf{T}$-matrix estimator can be minimized.
733 S7. PROBING THE BIAS INTENSITY WITH SPATIAL RECIPROCITY

734 In the accompanying paper, we use the scalar product $P_{in/out}$ between input and 735 output aberration phase laws to monitor the bias $|\delta T|^2$ of our $\mathbf{T}$-matrix estimator.

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

$$ \hat{T}(\mathbf{c}, \mathbf{r}_p) = \exp[j\{\phi(\mathbf{c}, \mathbf{r}_p) + \delta\phi(\mathbf{c}, \mathbf{r}_p)\}] \quad (S19) $$

738 with $T(\mathbf{c}, \mathbf{r}_p) = \exp[j\phi(\mathbf{c}, \mathbf{r}_p)]$ and $\delta\phi(\mathbf{c}, \mathbf{r}_p)$ the phase error of the estimator.

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

$$ |\delta T(\mathbf{c}, \mathbf{r}_p)|^2 = |1 - \exp[j\delta\phi(\mathbf{c}, \mathbf{r}_p)]|^2 = 4 \sin^2 \left[ \frac{\delta\phi(\mathbf{c}, \mathbf{r}_p)}{2} \right] \stackrel{\delta\phi \ll 1}{\sim} [\delta\phi(\mathbf{c}, \mathbf{r}_p)]^2 \quad (S20) $$

740 On the other hand, the scalar product $P_{in/out}$ is given by

$$ P_{in/out} = N_c^{-1} \sum_{\mathbf{c}} \exp[j\{\delta\phi_{in}(\mathbf{c}, \mathbf{r}_p) - \delta\phi_{out}(\mathbf{c}, \mathbf{r}_p)\}] \quad (S21) $$

741 In the previous equation, the sum over $\mathbf{c}$ can be replaced by an ensemble average 742 since $N_c = N_u \gg 1$:

$$ P_{in/out} = \langle \exp[j\{\delta\phi_{in}(\mathbf{c}, \mathbf{r}_p) - \delta\phi_{out}(\mathbf{c}, \mathbf{r}_p)\}] \rangle. \quad (S22) $$

743 Assuming a small phase error ($\delta\phi_{in/out} \ll 1$), the last equation can be rewritten 744 as follows

$$ P_{in/out} \approx 1 + j \langle \delta\phi_{in}(\mathbf{c}, \mathbf{r}_p) - \delta\phi_{out}(\mathbf{c}, \mathbf{r}_p) \rangle - \frac{\langle [\delta\phi_{in}(\mathbf{c}, \mathbf{r}_p) - \delta\phi_{out}(\mathbf{c}, \mathbf{r}_p)]^2 \rangle}{2}. \quad (S23) $$
745 Since $\langle \delta\phi_{\text{in/out}} \rangle = 0$ and $\langle \delta\phi_{\text{in}}\delta\phi_{\text{out}} \rangle = 0$, the last expression simplifies into

$$P_{\text{in/out}} \approx 1 - \frac{\langle |\delta\phi_{\text{in}}(\mathbf{c}, \mathbf{r}_{\text{p}})|^2 \rangle}{2} - \frac{\langle |\delta\phi_{\text{out}}(\mathbf{c}, \mathbf{r}_{\text{p}})|^2 \rangle}{2}. \qquad (\text{S24})$$

746 Assuming an equivalent phase error at input and output ($\langle |\delta\phi_{\text{in}}(\mathbf{c}, \mathbf{r}_{\text{p}})|^2 \rangle = \langle |\delta\phi_{\text{out}}(\mathbf{c}, \mathbf{r}_{\text{p}})|^2 \rangle$) finally leads to:

$$P_{\text{in/out}} \approx 1 - \langle |\delta\phi(\mathbf{c}, \mathbf{r}_{\text{p}})|^2 \rangle. \qquad (\text{S25})$$

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

$$P_{\text{in/out}} \approx 1 - \langle |\delta T(\mathbf{c}, \mathbf{r}_{\text{p}})|^2 \rangle. \qquad (\text{S26})$$

749 $P_{\text{in/out}}$ is thus a relevant quantity to estimate the bias intensity (see Fig. 3b of the 750 accompanying paper).

751 S8. MULTI-SCALE ANALYSIS OF WAVE DISTORTIONS

752 Supplementary Figure S4 demonstrates the benefit of a multi-scale analysis of 753 wave distortions with a gradual decrease of spatial windows $\mathcal{W}$ at each step of 754 the UMI process [Supplementary Figure S4a]. To that aim, this aberration correction scheme is compared with a direct estimation of the $\mathbf{T}$-matrix over the 755 smallest patches $\mathcal{W}$ [Supplementary Figure S4d]. The estimated transmission matrices $\hat{\mathbf{T}}$ differ in both cases (see comparison between Supplementary Figures S4b 756 and e) especially in the fat layer. The RPSFs obtained after phase conjugation of 757 $\hat{\mathbf{T}}$ demonstrate the benefit of the multi-scale analysis [Supplementary Figure S4c] 758 compared with a direct local investigation of wave distortions [Supplementary Figure S4f]. The fat area is actually the most aberrated in the field-of-view (see initial 759 RPSFs displayed by Fig. 2b of the accompanying paper). The initial coherence
LittlebyLittlevsDirectALex.png

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 ($\hat{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 ($\hat{T}$) extracted by a local analysis of D. (f) RPSFs after local aberration compensation.

<footnote n="763">763 factor *C* is thus much smaller in this area, which induces a strong bias on *T* when</footnote>

<footnote n="764">764 wave distortions are investigated over a reduced isoplanatic patch. On the con-</footnote>
trary, 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.

<figure><img src="image_11.png" /><figcaption>FIG. S5. Contrast & resolution enhancement in the pork chop experiment.
(a) Maps of local RPSF ($z = 29$ mm). (b) Local RPSF on the top right of the field-of-view. (c) Map of RPSF after the UMI process. (d) Corrected RPSF on the top right of the field-of-view. The resolution is evaluated at -3dB (see Methods in the accompanying paper). The contrast $F$ is the ratio between the confocal peak and the multiple scattering/noise background (see also Methods).</figcaption></figure>

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.
<figure><img src="image_12.png" /><figcaption>FIG. S6. Multi-scale compensation of wave distortions in the head phantom. (a) Successive patches used to perform a multi-scale analysis of wave distortions. (b) Radial profile of the RPSF amplitude at each step for three different depths (From left to right: z = 20, z = 32 and z = 60 mm). (c) Resolution as a function of depth at each step of correction (from red to green). At large depth (red dashed line), initial resolution can not be extracted as the incoherent background is larger than 1/2 as shown in panel (b) for z = 60 mm.</figcaption></figure>

S9. CONFOCAL FILTER

Supplementary Figure S7 shows the effect of the confocal filter on the $T$-matrix 779 estimation. The output aberration phase laws contained in $\hat{T}_{\text{out}}$ look much more 780 noisy in absence of an adaptive confocal filter (see the comparison between Supplementary Figures S7a and b). As shown by the scalar product between input and 781 output aberration phase laws [Supplementary Figure S7c], this “noise” comes from 782 the imperfect convergence of $\hat{T}$ towards $T$. Without any confocal filter, multiple 783 scattering drastically reduces the coherence factor and induces a strong bias on 784 estimation of $T$ (see Supplementary Section S5). On the contrary, the adaptive 785 confocal filter enables an enhancement of this coherence factor $C$ to ensure a sat- 786 isfactory estimation of $T$. The high degree of correlation between $\hat{T}_{\text{in}}$ and $\hat{T}_{\text{out}}$ 787
<figure><img src="image_13.png" /><figcaption>FIG. S7. Confocal filter in transcranial imaging. (a,b) Output aberration phase laws ($\hat{T}_{\text{out}}$) extracted without and with a confocal filter. (c,d) Normalised scalar products $P_{\text{in/out}}$ without and with a confocal filter, respectively. (e,f) RPSFs obtained with UMI without and with a confocal filter. (g) Box plot corresponding to the panels (c,d). Experimental data shown in this figure correspond to the head phantom experiment described in the accompanying paper ($z = 50mm$).</figcaption></figure>

789 proves this last assertion [Supplementary Figure S7d]. The effect of the confocal filter is also particularly obvious when looking at the RPSF obtained at the end of the UMI process. While a strong incoherent background subsists on the lateral parts of the field-of-view when no confocal filter is applied [Supplementary Figure S7e], a homogeneous focusing quality is obtained with the confocal filter [Supplementary Figure S7f].
<figure><img src="image_14.png" /><figcaption>FIG. S8. Illumination sequence. (a-e) Representation of different plane wave illumination sequence in the k-space. (f-j) Aliasing effect exhibited by the RPSFs due to incompleteness of illumination sequence displayed in panels a-e, respectively. These RPSFs have been measured in a speckle area of a tissue-mimicking phantom.</figcaption></figure>

S10. ILLUMINATION BASIS

Supplementary Figure S8 shows the impact of the illumination sequence on UMI. If the input illumination basis is complete [Supplementary Figure S8a], the RSPF exhibits the expected diffraction-limited resolution [Supplementary Figure S8f]. The side lobes along the y-axis are due to the probe geometry made of four blocks of transducers separated by a distance of 0.5 mm (three inactive rows of transducers along the y-axis).

When the number of illuminating plane waves is reduced [Supplementary Figures S8b-e], spatial aliasing occurs on corresponding RPSFs [Supplementary Figures S8g-j]. The maximal extension $\Delta\rho_{max}$ of the RPSF has to be fixed to avoid the spatial aliasing induced by the incompleteness of the plane wave illumination basis; $\Delta\rho_{max}$ is inversely proportional to the angular step $\delta\theta$ of the plane wave
807 illumination basis:

$$ \Delta\rho_{\max} \sim \lambda_c/(2\delta\theta) \qquad (S27) $$

808 with $\lambda_c$ the central wavelength and $\delta\theta$ the angular pitch used for the illumination 809 sequence. Thus, to avoid spatial aliasing, the coefficients $R(\rho_{in}, \rho_{out}, z)$ associated 810 with a transverse distance $|\rho_{out} - \rho_{in}|$ larger than the superior bound $\Delta\rho_{max}$ should 811 be filtered via a confocal filter.

812 Equation S27 implies the necessity of recording a high-dimension $\mathbf{R}$-matrix for 813 transcranial imaging, as aberrations are particularly important in that configura- 814 tion (see Fig. 5 of the accompanying paper). The number of independent incident 815 waves should scale as the number of resolution cells over which the RPSF spreads.

816 S11. DISCRIMINATE MULTIPLE SCATTERING FROM ELECTRONIC 817 NOISE

818 We consider here the background of the focused reflection matrix for a given 819 point $\mathbf{r}_p$:

$$ B(\Delta\rho, \mathbf{r}_p) = \langle R_M(\Delta\rho, \mathbf{r}_m) D(\Delta\rho) W(\mathbf{r}_m - \mathbf{r}_p) \rangle_{\mathbf{r}_m} \qquad (S28) $$

820 where $D(\Delta\rho)$ is a de-scanned window function that eliminates the confocal peak 821 and $W$ is a spatial average window function around the targeted focal point $\mathbf{r}_p$.

822 The background can be decomposed as the sum of a fully symmetric matrix 823 associated to multiple scattering (due to spatial reciprocity) and a fully random 824 matrix associated to the electronic noise as follows:

$$ \underbrace{\mathbf{B}}_{\text{Background}} = \underbrace{\mathbf{M}}_{\text{Multiple scattering}} + \underbrace{\mathbf{N}}_{\text{Noise}} \qquad (S29) $$

825 Projecting the $\mathbf{B}$-matrix onto its anti-symmetric subspace directly holds the anti-
826 symmetric part of the electronic noise such that:

$$ \mathbf{B}^{(A)} = \frac{\mathbf{B} - \mathbf{B}^{\top}}{2} = \mathbf{N}^{(A)} \qquad (S30) $$

Assuming equi-repartition of the electronic noise onto its symmetric and anti-symmetric subspace leads to:

$$ \|\mathbf{B}^{(A)}\|^2 = \|\mathbf{N}^{(A)}\|^2 = \frac{1}{2}\|\mathbf{N}\|^2 \qquad (S31) $$

827 The norm of the background can be expressed as follows:

$$ \|\mathbf{B}\|^2 = \|\mathbf{M}\|^2 + \|\mathbf{N}\|^2 + 2 \underbrace{\langle \mathbf{M} | \mathbf{N} \rangle}_{\sim 0} \qquad (S32) $$

Assuming that the scalar product between the electronic noise and the multiple scattering is zero on average, the multiple scattering rate $\alpha_M$ can be derived by combining equations (S31) & (S32):

$$ \alpha_M = \frac{\|\mathbf{M}\|^2}{\|\mathbf{B}\|^2} = 1 - 2 \underbrace{\frac{\|\mathbf{B}^{(A)}\|^2}{\|\mathbf{B}\|^2}}_{\beta} \qquad (S33) $$

828 with $\beta$ the anti-symmetric rate of the B-matrix.

## 829 S12. NOTATION AND SYMBOLS

## 830 References.

831 [1] W. Lambert, L. A. Cobus, M. Couade, M. Fink, and A. Aubry, Reflection Matrix

832 Approach for Quantitative Imaging of Scattering Media, Phys. Rev. X 10, 021048

833 (2020).
<table><caption>TABLE S5. **List of symbols for matrix imaging.**</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td><b>R</b></td><td>Reflection matrix</td></tr><tr><td><b>H</b></td><td>Point spread function matrix</td></tr><tr><td><i>RPSF</i></td><td>Reflection point spread function</td></tr><tr><td><b>G</b></td><td>Propagation matrix</td></tr><tr><td><b>D</b></td><td>Distortion matrix</td></tr><tr><td><b>C</b></td><td>Correlation matrix</td></tr><tr><td><i>δC</i></td><td>Perturbation term of <b>C</b></td></tr><tr><td><b>T</b> and <b>T̂</b></td><td>Transmission matrix and its estimator</td></tr><tr><td>|<i>δT</i>|<sup>2</sup></td><td>Bias intensity of <b>T</b>-matrix estimator</td></tr><tr><td><i>P</i><sub>in/out</sub></td><td>Scalar product between <b>T̂</b><sub>in</sub> and <b>T̂</b><sub>out</sub></td></tr><tr><td><b>i</b></td><td>Illumination basis</td></tr><tr><td><b>c</b></td><td>Correction basis</td></tr><tr><td><b>u</b></td><td>Transducer basis</td></tr><tr><td><b>k</b></td><td>Fourier basis</td></tr><tr><td><b>θ</b></td><td>Plane wave basis</td></tr><tr><td><i>l</i><sub>c</sub></td><td>Confocal filter size</td></tr><tr><td>ITR</td><td>Iterative Time Reversal</td></tr><tr><td>IPR</td><td>Iterative Phase Reversal</td></tr><tr><td><b>W</b><sup>(n)</sup></td><td>Wave-front of the ITR process at iteration <i>n</i></td></tr><tr><td><b>r</b><sub>m</sub></td><td>Common midpoint</td></tr><tr><td><b>r</b><sub>p</sub></td><td>Central point of a patch</td></tr><tr><td><b>Δρ</b> = <i>ρ</i><sub>out</sub> - <i>ρ</i><sub>in</sub></td><td>Distance input/output focusing points</td></tr><tr><td><i>D</i>(<b>Δρ</b>)</td><td>De-scanned window function</td></tr></tbody></table>

834 [2] V. Ntziachristos, Going deeper than microscopy: The optical imaging frontier in biology, Nat. Methods 7, 603 (2010).

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

837 [4] H. W. Babcock, The possibility of compensating astronomical seeing, Publ. As- tron. Soc. Pac. 65, 229 (1953).

839 [5] F. Roddier, ed., *Adaptive Optics in Astronomy* (Cambridge University Press, Cam- bridge, 1999).

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

842
<table><caption>TABLE S6. Matrix notations.</caption><thead><tr><th>Basis</th><th>Symbol</th><th>Adapted for</th></tr></thead><tbody><tr><td>Acquisition basis</td><td>R<sub>iu</sub>(t) = [R(i<sub>in</sub>, u<sub>out</sub>, t)]</td><td>Data recording</td></tr><tr><td>Focused basis</td><td>R<sub>ρρ</sub>(z) = [R(ρ<sub>in</sub>, ρ<sub>out</sub>, z)]</td><td rowspan="2">Focusing quality and multiple scattering quantification<sup>19</sup></td></tr><tr><td>Common midpoint</td><td>R<sub>M</sub>(z) = [R(Δρ, ρ<sub>m</sub>, z)]</td></tr><tr><td>Dual basis (input)</td><td>R<sub>cr</sub> = [R(c<sub>in</sub>, r<sub>out</sub>)]<br>D<sub>cr</sub> = [D(c<sub>in</sub>, r<sub>out</sub>)]<br>C<sub>in</sub> = [C(c<sub>in</sub>, c'<sub>in</sub>)]<br>T̂<sub>in</sub> = [T̂(c<sub>in</sub>, r<sub>p</sub>)]</td><td rowspan="2">Local aberration compensation<sup>28</sup></td></tr><tr><td>Dual basis (output)</td><td>R<sub>rc</sub> = [R(r<sub>in</sub>, c<sub>out</sub>)]<br>D<sub>rc</sub> = [D(r<sub>in</sub>, c<sub>out</sub>)]<br>C<sub>out</sub> = [C(c<sub>out</sub>, c'<sub>out</sub>)]<br>T̂<sub>out</sub> = [T̂(r<sub>p</sub>, c<sub>out</sub>)]</td></tr></tbody></table>

<table><caption>TABLE S7. Mathematical symbols.</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td>×</td><td>Matrix product</td></tr><tr><td>○</td><td>Hadamard product</td></tr><tr><td>⊗</td><td>Convolution product</td></tr><tr><td>†</td><td>Transpose conjugate of a matrix</td></tr><tr><td>⊤</td><td>Matrix transpose</td></tr><tr><td>^</td><td>Estimator of a physical quantity</td></tr><tr><td>SVD</td><td>Singular Value Decomposition</td></tr><tr><td>U<sup>(i)</sup></td><td>i<sup>th</sup> right singular vector of a matrix</td></tr><tr><td>V<sup>(i)</sup></td><td>i<sup>th</sup> left singular vector of a matrix</td></tr><tr><td>σ<sub>i</sub></td><td>i<sup>th</sup> singular value of a matrix</td></tr><tr><td>⟨· · ·⟩</td><td>Ensemble average</td></tr></tbody></table>

<footnote n="843">843 Freq. Control **35**, 768 (1988).</footnote>

<footnote n="844">844 [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).</footnote>

<footnote n="847">847 [8] R. Mallart and M. Fink, Adaptive focusing in scattering media through sound-speed inhomogeneities: The van Cittert Zernike approach and focusing criterion, J.</footnote>
<table><caption>TABLE S8. **List of general symbols.**</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td><i>I</i></td><td>Image ⇔ Estimation of the reflectivity</td></tr><tr><td><b>r</b> = (x, y, z)</td><td>Focal point</td></tr><tr><td><b>ρ</b> = (x, y)</td><td>Transverse coordinate</td></tr><tr><td>λ<sub>c</sub></td><td>Wavelength at the central frequency</td></tr><tr><td>f<sub>s</sub></td><td>Sampling frequency</td></tr><tr><td>f<sub>c</sub></td><td>Central frequency</td></tr><tr><td>c<sub>0</sub></td><td>Speed-of-sound hypothesis</td></tr><tr><td><b>u</b> = (u<sub>x</sub>, u<sub>y</sub>, 0)</td><td>Transducer position</td></tr><tr><td>δρ<sub>0</sub></td><td>Transverse ideal resolution</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>τ</td><td>Time-of-flight</td></tr><tr><td>Δτ</td><td>Time-delay</td></tr><tr><td>γ</td><td>Medium reflectivity</td></tr><tr><td><b>θ</b> = [θ<sub>x</sub>, θ<sub>y</sub>]</td><td>Plane wave</td></tr><tr><td><b>k</b> = [k<sub>x</sub>, k<sub>y</sub>]</td><td>Fourier basis</td></tr><tr><td>β</td><td>Anti-symmetric rate of a matrix</td></tr><tr><td>θ<sub>max</sub></td><td>Directivity of transducers</td></tr><tr><td>δθ</td><td>Plane wave sampling</td></tr><tr><td><b>Δu</b> = (Δu<sub>x</sub>, Δu<sub>y</sub>)</td><td>Probe dimension</td></tr><tr><td><i>C</i></td><td>Coherence factor</td></tr><tr><td><i>A</i><sub>(-3dB)</sub></td><td>Area above 1/2 on RPSF amplitude</td></tr><tr><td>δρ<sub>(-3dB)</sub></td><td>Experimental RPSF resolution</td></tr><tr><td>δρ<sub>0</sub></td><td>Diffraction-limited resolution</td></tr><tr><td><i>F</i></td><td>RPSF contrast</td></tr><tr><td>α<sub>S</sub></td><td>RPSF single scattering rate</td></tr><tr><td>α<sub>M</sub></td><td>RPSF multiple scattering rate</td></tr><tr><td>α<sub>N</sub></td><td>RPSF electronic noise rate</td></tr><tr><td>α<sub>B</sub></td><td>RPSF background rate</td></tr><tr><td><i>W</i></td><td>Spatial average window function</td></tr><tr><td><i>N</i><sub>W</sub></td><td>Number of resolution cells in <i>W</i></td></tr><tr><td><b>w</b> = (w<sub>ρ</sub>, w<sub>z</sub>) = ({w<sub>x</sub>, w<sub>y</sub>}, w<sub>z</sub>)</td><td>Dimension of <i>W</i></td></tr><tr><td><i>A</i></td><td>Apodization term of synthetic aperture</td></tr></tbody></table>
[9] R. Ali, T. Brevett, L. Zhuang, H. Bendjador, A. S. Podkowa, S. S. Hsieh, W. Simson, S. J. Sanabria, C. D. Herickhoff, and J. J. Dahl, Aberration correction in diagnostic ultrasound: A review of the prior field and current directions, Z. Med. Phys. (in press) 10.1016/j.zemedi.2023.01.003 (2023).

[10] S.-E. Måsøy, T. Varslot, and B. Angelsen, Iteration of transmit-beam aberration correction in medical ultrasound imaging, J. Acoust. Soc. Am. **117**, 450 (2005).

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

[12] M. Pernot, M. Tanter, and M. Fink, 3-D real-time motion correction in high-intensity focused ultrasound therapy, Ultrasound Med. Biol. **30**, 1239 (2004).

[13] M. Jaeger, E. Robinson, H. G. Akarøy, and M. Frenz, Full correction for spatially distributed speed-of-sound in echo ultrasound based on measuring aberration delays via transmit beam steering, Phys. Med. Biol. **60**, 4497 (2015).

[14] G. Chau, M. Jakovljevic, R. Lavarello, and J. Dahl, A Locally Adaptive Phase Aberration Correction (LAPAC) Method for Synthetic Aperture Sequences, Ultrason. Imaging **41**, 3 (2019).

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

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

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

[18] H. Bendjador, T. Deffieux, and M. Tanter, The SVD Beamformer: Physical Principles and Application to Ultrafast Adaptive Ultrasound, IEEE Trans. Med. Imag.
877 **39**, 3100 (2020).

878 [19] W. Lambert, J. Robin, L. A. Cobus, M. Fink, and A. Aubry, Ultrasound matrix imaging – Part I: The focused reflection matrix, the F-factor and the role of multiple scattering, IEEE Trans. Med. Imag. **41**, 3907 (2022).

881 [20] S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K. Song, J. H. Park, J. H. Lee, M. J. Kim, K. H. Kim, and W. Choi, High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering, Nat. Commun. **8**, 2157 (2017).

885 [21] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, Distortion matrix concept for deep optical imaging in scattering media, Sci. Adv. **6**, eaay7170 (2020).

888 [22] S. Yoon, H. Lee, J. H. Hong, Y.-S. Lim, and W. Choi, Laser scanning reflection-matrix microscopy for aberration-free imaging through intact mouse skull, Nat. Commun. **11**, 5721 (2020).

891 [23] Y. Kwon, J. H. Hong, S. Kang, H. Lee, Y. Jo, K. H. Kim, S. Yoon, and W. Choi, Computational conjugate adaptive optics microscopy for longitudinal through-skull imaging of cortical myelin, Nat. Commun. **14**, 105 (2023).

894 [24] U. Najar, V. Barolle, P. Balondrade, M. Fink, A. C. Boccara, M. Fink, and A. Aubry, Non-invasive retrieval of the transmission matrix for optical imaging deep inside a multiple scattering medium, arXiv: 2303.06119 (2023).

897 [25] T. Blondel, J. Chaput, A. Derode, M. Campillo, and A. Aubry, Matrix Approach of Seismic Imaging: Application to the Erebus Volcano, Antarctica, J. Geophys. Res. : Solid Earth **123**, 10936 (2018).

900 [26] R. Touma, T. Blondel, A. Derode, M. Campillo, and A. Aubry, A distortion matrix framework for high-resolution passive seismic 3-D imaging: Application to the San Jacinto fault zone, California, Geophy. J. Int. **226**, 780 (2021).
[27] T. I. Sommer and O. Katz, Pixel-reassignment in ultrasound imaging, Appl. Phys. Lett. **119**, 123701 (2021).

[28] W. Lambert, L. A. Cobus, J. Robin, M. Fink, and A. Aubry, Ultrasound matrix imaging – Part II: The distortion matrix for aberration correction over multiple isoplanatic patches, IEEE Trans. Med. Imag. **41**, 3921 (2022).

[29] N. M. Ivancevich, J. J. Dahl, G. E. Trahey, and S. W. Smith, Phase-aberration correction with a 3-D ultrasound scanner: Feasibility study, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **53**, 1432 (2006).

[30] J. Lacefield and R. Waag, Time-shift estimation and focusing through distributed aberration using multirow arrays, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **48**, 1606 (2001).

[31] B. D. Lindsey and S. W. Smith, Pitch-catch phase aberration correction of multiple isoplanatic patches for 3-D transcranial ultrasound imaging, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **60**, 463 (2013).

[32] D.-L. Liu and R. Waag, Estimation and correction of ultrasonic wavefront distortion using pulse-echo data received in a two-dimensional aperture, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **45**, 473 (1998).

[33] L. Ratsimandresy, P. Mauchamp, D. Dinet, N. Felix, and R. Dufait, A 3 MHz two dimensional array based on piezocomposite for medical imaging, in *2002 IEEE Ultrasonics Symposium, 2002. Proceedings.*, Vol. 2 (IEEE, Munich, Germany, 2002) pp. 1265–1268.

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

[35] J. Provost, C. Papadacci, C. Demene, J.-L. Gennisson, M. Tanter, and M. Pernot, 3-D ultrafast Doppler imaging applied to the noninvasive mapping of blood vessels in Vivo, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **62**, 1467 (2015).
[36] S. A. Goss, R. L. Johnston, and F. Dunn, Compilation of empirical ultrasonic properties of mammalian tissues. II, J. Acoust. Soc. Am. **68**, 93 (1980).

[37] L. M. Hinkelman, D. Liu, L. A. Metlay, and R. C. Waag, Measurements of ultrasonic pulse arrival time and energy level variations produced by propagation through abdominal wall, J. Acoust. Soc. Am. **95**, 530 (1994).

[38] J. Lacefield and R. Waag, Examples of design curves for multirow arrays used with time-shift compensation, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **49**, 1340 (2002).

[39] N. M. Ivancevich, G. F. Pinton, H. A. Nicoletto, E. Bennett, D. T. Laskowitz, and S. W. Smith, Real-time 3-D contrast-enhanced transcranial ultrasound and aberration correction, Ultrasound Med. Biol. **34**, 1387 (2008).

[40] A. Bertolo, M. Nouhoum, S. Cazzanelli, J. Ferrier, J.-C. Mariani, A. Kliewer, B. Belliard, B.-F. Osmanski, T. Deffieux, S. Pezet, Z. Lenkei, and M. Tanter, Whole-Brain 3D Activation and Functional Connectivity Mapping in Mice using Transcranial Functional Ultrasound Imaging, J. Vis. Exp. **168**, e62267 (2021).

[41] A. Chavignon, B. Heiles, V. Hingot, C. Orset, D. Vivien, and O. Couture, 3D Transcranial Ultrasound Localization Microscopy in the Rat Brain With a Multiplexed Matrix Probe, IEEE Trans. Biomed. Eng. **69**, 2132 (2022).

[42] C. Demené, J. Robin, A. Dizeux, B. Heiles, M. Pernot, M. Tanter, and F. Perren, Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients, Nat. Biomed. Imag. **5**, 219 (2021).

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

[44] J. Robin, C. Demené, B. Heiles, V. Blanvillain, L. Puke, F. Perren-Landis, and M. Tanter, In vivo Adaptive Focusing for Clinical Contrast-Enhanced Transcranial Ultrasound Imaging in Human, Phys. Med. Biol **68**, 025019 (2023).
957 [45] M. Tanter and M. Fink, Ultrafast imaging in biomedical ultrasound, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **61**, 102 (2014).

958 [46] M. Jaeger, G. Held, S. Peeters, S. Preisser, M. Gr´unig, and M. Frenz, Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: proof of principle, Ultrasound Med. Biol **41**, 235 (2015).

962 [47] M. Imbault, A. Faccinetto, B.-F. Osmanski, A. Tissier, T. Deffieux, J.-L. Gennisson, V. Vilgrain, and M. Tanter, Robust sound speed estimation for ultrasound-based hepatic steatosis assessment, Phys. Med. Biol **62**, 3582 (2017).

965 [48] M. Jakovljevic, S. Hsieh, R. Ali, G. Chau Loo Kung, D. Hyun, and J. J. Dahl, Local speed of sound estimation in tissue using pulse-echo ultrasound: Model-based approach, J. Acoust. Soc. Am. **144**, 254 (2018).

968 [49] A. Aubry and A. Derode, Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues, J. Acoust. Soc. Am. **129**, 225 (2011).

971 [50] C. Brütt, A. Aubry, B. Gérardin, A. Derode, and C. Prada, Weight of single and recurrent scattering in the reflection matrix of complex media, Phys. Rev. E **106**, 025001 (2022).

974 [51] C. Papadacci, M. Tanter, M. Pernot, and M. Fink, Ultrasound backscatter tensor imaging (BTI): analysis of the spatial coherence of ultrasonic speckle in anisotropic soft tissues, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **61**, 986 (2014).

977 [52] A. Rodriguez-Molares, A. Fatemi, L. Lovstakken, and H. Torp, Specular Beam-forming, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **64**, 1285 (2017).

979 [53] D. Zhao, L. N. Bohs, and G. E. Trahey, Phase aberration correction using echo signals from moving targets i: Description and theory, Ultrason. Imaging **14**, 97 (1992).

982 [54] B.-F. Osmanski, G. Montaldo, M. Tanter, and M. Fink, Aberration correction by time reversal of moving speckle noise, IEEE Trans. Ultrason. Ferroelectr. Freq.
984 Control **59**, 1575 (2012).

985 [55] F. Berland, T. Fromenteze, D. Boudescoque, P. Di Bin, H. H. Elwan, C. Aupetit-Berthelemot, and C. Decroze, Microwave photonic mimo radar for short-range 3d imaging, IEEE Access **8**, 107326 (2020).

988 [56] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink, Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **56**, 489 (2009).

991 [57] V. Perrot, M. Polichetti, F. Varray, and D. Garcia, So you think you can DAS? A viewpoint on delay-and-sum beamforming, Ultrasonics **111**, 106309 (2021).

993 [58] M. Fink and C. Dorme, Aberration correction in ultrasonic medical imaging with time-reversal techniques, Int. J. Imaging Syst. Technol. **8**, 110 (1997).

995 [59] J. Mertz, H. Paudel, and T. G. Bifano, Field of view advantage of conjugate adaptive optics in microscopy applications, Appl. Opt. **54**, 3498 (2015).

997 [60] F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, and A. Aubry, Ultrasound matrix imaging [data]. Zenodo (2023).

999 [61] S. Silverstein, Ultrasound scattering model: 2-d cross-correlation and focusing criteria-theory, simulations, and experiments, IEEE Trans. Ultrason. Ferroelectr. Freq. Control **48**, 1023 (2001).