William Lambert

Laura A Cobus

Justine Robin

Mathias Fink

Alexandre Aubry

Ultrasound Matrix Imaging-Part II. The distortion matrix for aberration correction over multiple isoplanatic patches

I. INTRODUCTION

In most ultrasound imaging, the human body is insonified by a series of incident waves. The medium reflectivity is then estimated by detecting acoustic backscatter from short-scale variations of the acoustic impedance. An image (spatial map) of reflectivity is commonly constructed using delay-and-sum beamforming (DAS). In this process, echoes coming from a particular point, or image pixel, are selected by summing the signals generated by this echo at the aperture, thereby accounting -for each element -for the respective time-of-flight associated with the forward and return travel paths of the ultrasonic wave between the probe and that point. The resulting signal is allocated at the corresponding pixel of the image, and the procedure repeated for each pixel. The time-of-flight between any incident wave and focal point is calculated with the assumption that the medium is homogeneous with a constant speed-of-sound. However, in human tissue, long-scale fluctuations of the acoustic impedance can invalidate this assumption [1]. The resulting incorrectly calculated times-of-flight (also called focusing laws) can lead to aberration of the associated image, meaning that resolution and contrast are strongly degraded. While adaptive focusing methods have been developed to deal with this issue, they rely on the hypothesis that the * alexandre.aubry@espci.fr speed-of-sound variations occur only in a thin screen at the probe aperture. However, this assumption is simply incorrect in soft tissues [2] such as fat, skin and muscle, in which the order of magnitude of acoustic impedance fluctuations is around 5% [START_REF] Duck | Physical properties of tissue: A comprehensive reference book[END_REF]. This causes higher-order aberrations which are only invariant over small regions, often referred to as isoplanatic patches. To tackle this issue, recent studies [START_REF] Bendjador | [END_REF]5] extract an aberration law for each image voxel by probing the correlation of the time-delayed echoes coming from adjacent focal spots. The aberration laws are estimated either in the time domain [START_REF] Bendjador | [END_REF][6][7][8] or in the Fourier domain [5,9], for different insonification sequences (focused beams [6,8], single-transducer insonification [5] or plane wave illumination [START_REF] Bendjador | [END_REF]). In all of these techniques, a focusing law is estimated in either the receive [6,9] or transmit [5,10] mode, but this law is then used to compensate for aberrations in both reflection and transmission. However, spatial reciprocity between input and output is only valid if the emission and detection of waves are performed in the same basis; in other words, the distortion undergone by a wavefront travelling to and from a particular point is only the same if the wave has interacted with the same heterogeneities in both directions. If this condition is not fulfilled, applying the same aberration phase law in transmit and receive modes may improve the image quality to some degree, but will not be optimal.

To obtain optimized focusing laws both in transmit and receive modes, these two steps must be considered separately. Recent studies [11][12][13] have shown how to de-couple the location of the transmit and receive focal spots to build a focused reflection (FR) matrix R. Containing the medium responses between virtual sources and virtual sensors located within the medium, this matrix is the foundation of ultrasound matrix imaging (UMI). Firstly, a focusing criterion and a background rate can be built from the FR matrix, which allows the mapping of the focusing quality and contrast over all pixels of the ultrasound image [14] in both speckle and specular regimes [11]. Secondly, a distortion matrix D can be built from the FR matrix for a local correction of high-order aberrations; this concept was first presented in optical imaging [15], then in ultrasound [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF], and most recently in seismology [START_REF] Touma | [END_REF]. Whereas R holds the wavefronts which are reflected from the medium, D contains the deviations from an ideal reflected wavefront which would be obtained in the absence of heterogeneities. It has been shown that, for specular reflectors [15], in sparse media [START_REF] Touma | [END_REF] and in the speckle regime for a multi-layered speed-of-sound distribution [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF], a singular value decomposition (SVD) of D yields a one-to-one association between each isoplanatic patch in the focal plane and the corresponding wavefront distortion in the far-field. This information then enables a correction of aberrations over multiple isoplanatic patches.

For laterally-varying aberrations, the D-matrix should be analyzed locally and investigated over limited spatial windows [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. In this paper, we generalize this approach to cope with in-vivo applications in which multiple scattering and high-order aberrations can greatly reduce the size of isoplanatic patches. In that respect, the FR matrix will play a pivotal role. First, it enables the application of an adaptive confocal filter to: (i ) reduce the detrimental influence of multiple scattering and/or noise; (ii ) ensure a local isoplanaticity in order to converge towards a satisfying aberration law. Second, working with the FR matrix allows easy projection of the input or output wave-fields into different bases (far-field, transducer plane, some intermediate surface, etc.) for optimal aberration correction. In contrast with previous studies that compensate aberrations from a single (transducer or plane wave) basis [START_REF] Bendjador | [END_REF]5,9, 10], we will show that alternating correction bases is particularly relevant for spatially-distributed aberrations. In each basis, local aberration phase laws can be estimated by adjusting the field-of-view covered by each local D-matrix with each isoplanatic patch. By theoretically modelling the SVD process, we show how to optimize the size of the targeted isoplanatic patch and the numerical confocal pinhole. The whole process can then be iterated by gradually reducing the size of the isoplanatic patches, thereby compensating for more and more complex aberrations throughout the iteration process.

The D-matrix formalism was first introduced in the context of ultrasound imaging with simple proof-ofconcept experiments [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF] in which the points described above (confocal filter, projection in complementary correction bases, theoretical modelling of the SVD process)

were not tackled. Here, we show that these important steps make UMI much more robust for challenging invivo applications. The overall method is validated by means of numerical simulations that mimic aberrations through the abdominal wall. As an experimental proofof-concept, we apply UMI to the complex case of in-vivo imaging of a gallbladder. A set of optimized focusing laws is obtained for each point of the medium, enabling (i ) a mapping of the isoplanatic modes in the field-of-view, revealing the arrangement of the different tissues in the medium, and (ii ) the calculation of an ultrasound image with optimal contrast and enhanced resolution over a major part of the field-of-view. The drastic improvement compared to the conventional ultrasound image is quantified by means of the focusing F -factor and the incoherent background rate introduced in the first paper of the series [14].

II. THE FOCUSED REFLECTION MATRIX A. Experimental procedure

The experiment was performed on a healthy volunteer (this study is in conformation with the declaration of Helsinki). The probe was placed in a subcostal position in order to image the gallbladder along its short axis. In this configuration, the gallbladder appears as a circular ring. The experimental procedure consists in recording the reflection matrix R using a standard plane wave sequence [18]. The acquisition was performed using a medical ultrafast ultrasound scanner (Aixplorer Mach-30, Supersonic Imagine, Aix-en-Provence, France) driving a 2 -10 MHz linear transducer array containing N = 192 transducers with a pitch p = 0.2 mm (SL10-2, Supersonic Imagine). The ultrasound sequence consisted in transmitting 101 steering angles spanning from -25 o to 25 o with steps of δθ = 0.5 o , calculated assuming a constant speed of sound c 0 = 1540 m/s. This choice was made by optimizing the spatially-averaged focusing factor F obtained for a set of speed-of-sound hypotheses ranging from c 0 = 1450 -1600 m/s [11,14]. The emitted signal was a sinusoidal burst of three half periods of the central frequency f c = 5 MHz, with pulse repetition frequency 1000 Hz. For each excitation, the back-scattered signal was recorded by all probe elements over time ∆t = 80 µs with a sampling frequency f s = 40 MHz. Mathematically, we write this set of acoustic responses as R uθ (t) ≡ [R(u out , θ in , t)], where u out defines the coordinate of the receiving transducer, θ in the angle of incidence and t the time-of-flight. Subscripts 'in' and 'out' denote propagation in the forward and backward directions, respectively. Note that the coefficients of R uθ should be complex, as they contain the amplitude and phase of the medium response. If the responses R uθ (t) are not complex modulated RF signals, then the corresponding analytic signals should be considered.

B. Computing the focused reflection matrix

In conventional ultrasound imaging, the reflectivity of a medium at a given point is estimated by (i ) focusing a wave on this point, thereby creating a virtual source, and (ii ) coherently summing the echoes coming from that same point, thus synthesizing a virtual receiver at that location. In UMI, this focusing operation is performed in post-processing, and the input/output focusing points, r in /r out , are decoupled [Fig. 1(a)]. This is the principle of the broadband focused reflection matrix R rr = R(r in , r out ) containing the responses between virtual sources and sensors located throughout the medium. Note that the concept of virtual transducers is here mainly didactic and that, of course, a virtual transducer does not act as a real source or sink of energy. Moreover, they are highly directive; the directivity pattern points downwards for a virtual source and upwards for a virtual receiver.

In the first article of this series [14], the FR matrix was built using a beamforming process in the temporal Fourier domain. Here, we show that this matrix of complex coefficients can also be directly (and more quickly) computed in the time domain from the recorded IQ signals via conventional DAS beamforming: (1) where τ in and τ out are the transmit and receive focusing laws such that τ in (θ in , r in ) = [x in sin(θ in ) + z in cos(θ in )]/c 0 , and (2a)

R(r out , r in ) = θin,uout A(u out , θ in , r in , r out ) R(u out , θ in , τ in (θ in , r in ) + τ out (u out , r out )),
τ out (u out , r out ) = |x out -u out | 2 + z 2 out /c 0 . (2b)
A is an apodization factor that limits the extent of the receive synthetic aperture, and (x in , z in ) and (x out , z out ) are the coordinates of the input and output focusing points r in and r out , respectively. In this paper, we will restrict our study to the x-projection of R rr , written R xx (z), in which only the responses between virtual transducers located at the same depth are considered (z = z in = z out ).

C. Manifestation of aberrations and multiple scattering

Each row of R xx (z) corresponds to the situation in which waves have been focused at r in = (x in , z) in transmission, and virtual detectors at r out = (x out , z) record the resulting spatial wave spreading across the focal plane [Fig. 1(a)]. Fig. 1(b) shows R xx (z) at z = 39 mm. Note that the coefficients R(x out , x in , z) associated with a transverse distance |x out -x in | larger than a superior bound ∆x max ∼ λ max /(2δθ) are not displayed because of spatial aliasing [14]. The angle increment δθ could have been decreased to remove this spatial aliasing; however, this would be at the cost of a longer measurement time which is critical for in-vivo imaging since the medium is moving. In the single scattering approximation, the FR matrix coefficients can be expressed theoretically as follows [14]:

R(x out , x in , z) = dxH out (x, x out , z)γ(x, z)H in (x, x in , z), (3)
where γ(x, z) is the medium reflectivity. H in (x, x in , z) and H out (x, x out , z) are the transmit and receive PSFs, i.e. the spatial amplitude distribution of the input and output focal spots at depth z. In the absence of aberration, their spatial extent δx 0 (r in/out ) scales as λ/[2 sin β(r in/out )], where β(r) is the maximum angle of illumination (in transmit mode) or collection (in receive mode) by the array from the associated focal point r.

The intensity distribution along the diagonal of R xx (z) yields a line of the confocal image I that would be obtained by plane wave synthetic beamforming [18]:

I(r) ≡ R (r, r) 2 . ( 4
)
The experimental ultrasound image is displayed in Figs. 2(a) and (b). In Fig. 2(a), the blurred appearance of the gallbladder internal wall indicates the presence of aberrations and/or multiple scattering. We now quantify these detrimental effects by investigating the off-diagonal coefficients of R xx .

In the accompanying paper [14], the intensity of the antidiagonals of R xx (z), referred to as the common midpoint (CMP) intensity profile, is shown to give access to

the local input-output incoherent PSF, |H in | 2 ∆x
⊛ |H out | 2 in the speckle regime (the symbol ⊛ here stands for a convolution product). Fig. 1(g) shows the CMP intensity profile in the area B displayed in Fig. 2(c). A confocal peak can be observed, originating from single scattering, which sits on a wider incoherent background resulting from high-order aberrations, multiple scattering and experimental noise. By applying the fitting procedure developed in the accompanying paper [14], maps of the focusing factor F (r) and of the incoherent background rate α(r) were extracted from local CMP intensity profiles (see Figs. 2(e) and (g), respectively). Blue areas in Fig. 2(e)[F (r) ∼ 1)] indicate high reliability; the ultrasound image accurately describes the medium reflectivity (see Fig. 2(a,c)). Green and yellow areas [F (r) < 0.5] indicate aberrated areas of the image . Gray areas correspond to the situation where the incoherent background is too large (α > 0.75) to obtain a reliable estimation of the focusing factor F . This occurs when the intensity level of the backscattered signal generated by the region of interest is threefold lower than the multiple scattering and/or electronic noise contributions (see yellow regions in Fig. 2(g)). Causes of this low intensity level could be weak reflectivity of the medium in the probed region, and/or strong fluctuations of the speed-of-sound upstream of (at shallower depths than) the focal plane, which would decrease the relative single scattering level close to the diagonal. Overall, the ultrasound image displayed in Fig. 2(c) suffers from: (i ) a focusing quality that is far from ideal in many areas (especially at shallow depths and on the right part of the image); (ii ) a weak contrast due to a predominant incoherent background (especially at large depths and on the left part of the image). Thus, prior to performing aberration correction, it is important to remove as much incoherent background as possible from the FR matrix.

D. Filtering multiple scattering and noise

The incoherent background can be partially suppressed using an adaptive confocal filter [12,13]. This process consists in weighting the coefficients R(x in , x out , z) of the FR matrix as a function of the distance |x out -x in | between the virtual transducers, such that:

R ′ (x out , x in ) = R(x out , x in ) exp - | x out -x in | 2 2l 2 c (r) . (5)
This filter has a Gaussian shape, with a width l c (r) that needs to be carefully set. If l c (r) is too large, the multiply-scattered echoes will prevent a correct estimation of the aberration phase law. If l c (r) is too small, the filter then acts as an apodization function that will smooth out the resulting aberration phase law. To be efficient, l c (r) is automatically chosen to scale with δx 0 (r) (hence adaptive): l c (r) = nδx 0 (r), with n a coefficient whose value is reported in Table . I. The choice of this value will be justified in Sec. III D.

Figs. 1(b,c) show the original and filtered FR matrices, R xx and R ′ xx , respectively, for z = 39 mm. It can be seen that the adaptive confocal filter has removed part of the incoherent background. However, R ′ xx still contains a residual multiple scattering and/or noise component which also pollutes matrix coefficients very close to the diagonal. Note that this filter has no impact on the raw ultrasound image since the confocal signals are unaffected. However, it constitutes a necessary step for an unbiased estimation of the aberration law.

III. MATRIX CORRECTION OF ABERRATIONS

The FR matrix R ′ xx (z) is now used to implement the distortion matrix concept [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. In this section, we will show how to estimate and correct for aberrations successively in the transmit and receive modes, both from the far-field and the surface of the transducer array. We describe all of the technical steps of the aberration correction process: (i ) the projection of the FR matrix at output or input into a correction basis (here either the Fourier or transducer basis) in order to investigate the reflected or incident wave-front associated with each virtual source or transducer, respectively, (ii ) the realignment of the transmitted or reflected wave-fronts to form the distortion matrix D, (iii ) the truncation of D into overlapping isoplanatic patches, (iv ) the singular value decomposition of D or of its normalized correlation matrix to extract an aberration phase law for each isoplanatic patch, and (v ) the application of the focusing law and back-projection of the reflection matrix into the focused basis. All of these steps are repeated by exchanging input and output bases, as well as the correction basis. Fig. 3 shows a workflow that sums up the different steps of the UMI procedure. The process is then be iterated while gradually reducing the size of isoplanatic patches in order to address higher order aberrations.

A. Projection of the reflection matrix into a dual basis

In adaptive focusing, the aberrating layer is often modeled as a phase screen. For an optimal correction, the ultrasonic data should be back-propagated to the plane containing this aberrating layer; indeed, from this plane, the aberration is spatially-invariant. By applying the phase conjugate of the aberration phase law, aberration can be fully compensated for at any point of the medium. However, in real life, speed-of-sound inhomogeneities are distributed over the entire medium and aberration can take place everywhere. To treat this case, the strategy here is to back-propagate ultrasound data into several planes from which the aberration phase law should be estimated and then compensated. The optimal correction plane is the one that maximizes the size of isoplanatic patches. For multi-layered media, the Fourier plane is the most adequate since plane waves are the propagation invariants in this geometry. For aberrations induced by superficial veins or skin nodules, the probe plane is a good choice. In this paper, the aberration correction will be performed in these two planes as they coincide also to the emission and reception bases used to record the reflection matrix. However, note that, in practice, other correction planes can be chosen according to the imaging problem.

Projection into the Fourier basis

To project the reflected wave-field into the Fourier plane, a spatial Fourier transform should be applied to Realign virtual sources (12) Spatial FT ( 6),( 8)

Time reversal analysis (SVD) (18)
Updating transmission matrices (36) Transducer basis Plane wave basis Spatial window filter around r p (14) Loop over

R uθ R rr R' rr R' / ur R' ru R' / kr R' rk D' ur (r p )/D' ru (r p ) D' kr (r p )/D' rk (r p ) R rr Q 1 T 1 U 1 (u,r p ) U 1 (k x ,r p )
Back-propagation ( 10),( 11)

D ur /D ru D kr /D rk feedback loop (c)
Matrix focusing (37)

r p u (a) (b) (c) T i T i † Q i Q i † x k x x FIG. 3. (a)
Workflow of the various steps of the aberration correction process. The blue lines refer to the estimation process of the aberration phase laws while the green lines refer to the aberration correction procedure. The purple line depicts the iteration of the UMI process. The gray rectangle symbolizes a loop over the spatial window centered on rp. (b,c) Schematic representation for the projection between the transducer/far-field and focused bases.

the output of R ′ xx (z) (Fig. 3(c)): R ′ kx (z) = T 0 × R ′ xx (z). ( 6
)
where T 0 is the Fourier transform operator

T 0 (k x , x) = exp (ik x x)/ N k , (7)
k x the transverse wave number, and N k the dimension of T 0 . On the one hand, the support ∆k of the wave-field in the Fourier plane is fixed by the resolution of the focal plane: ∆k ∼ 1/δx 0 . On the other hand, its sampling δk is the inverse of the field-of-view : δk ∼ 1/(N δx 0 ). This choice ensures a conservation of the information between the focused and Fourier bases, such that T 0 T † 0 = I, where I is the identity matrix and † stands for the transpose conjugate operation.

R ′ kx (z)≡ [R(k out , x in , z)] contains the set of aberrated wavefronts in the spatial Fourier domain generated by each virtual source r in . Fig. 4(c) shows the phase of R ′ kx (z) at z = 60 mm. Using the central frequency f c as a reference frequency, the transverse wave number k out can be associated with a plane wave of angle θ out , such that k out = k c sin(θ out ), with k c = 2πf c /c 0 . Expressing the far-field projection as a plane wave decomposition is useful to define the boundaries of this basis [white dashed lines in Fig. 4(a)]; the maximum transverse wave number is

k max ∼ k c sin[β(r)] [see Fig. 4(k)].
The matrix R ′ kx (z) will be used to tackle aberrations in the receive plane wave basis. To do the same in the transmit basis, a reciprocal projection to that of ( 6) can be performed at the input of R ′ xx (z):

R ′ xk (z) = R ′ xx (z) × T ⊤ 0 . ( 8
)
where the symbol ⊤ stands for matrix transpose. The coefficients R ′ (x out , k in , z) correspond to the wave-field probed by the virtual transducer at r out if a plane wave of transverse wave number k in illuminates the medium. This matrix will be used to investigate aberrations in the transmit plane wave basis.

Projection into the transducer basis

The strategy to treat aberrations in the transducer plane is similar to that described above. Free-space transmission matrix, Q 0 = [Q 0 (r, u)], is defined between the focused and transducer bases at central frequency f c :

Q 0 = T † 0 × (P • T 0 ) , (9)
where the symbol • stands for a Hadamard product and P = [P (k x , z)] is the plane wave propagator at the central frequency:

P (k x , z) = e i √ k 2 c -k 2
x z . The transmission matrix Q 0 can be considered here at the single central frequency because of the time gating operation performed in (1). If evanescent waves are neglected, the sampling δu in the transducer basis can be fixed to λ c /2, with λ c = c 0 /f c the wavelength at the central frequency f c , such that Q 0 Q † 0 = I. The operator Q 0 can be given a physical interpretation by reading the terms of ( 9) from right to left: (i) a spatial Fourier transform using the operator T 0 to project the wave-field in the plane wave basis; (ii) the plane wave propagation modeled by the propagator P between the focal and transducer planes over a distance z; (iii) an inverse Fourier transformation T † 0 that projects the wave-field into the transducer basis. Note that the coefficients of Q 0 actually correspond to the z-derivative of the free space Green's functions between transducers and focusing points at the central frequency f c .

Using this operator Q, the matrix R ′ xx can be projected into the transducer basis (Fig. 3(b)) either at input,

R ′ xu (z) = R ′ xx (z) × Q 0 ⊤ ( 10
)
or output,

R ′ ux (z) = Q 0 × R ′ xx (z). ( 11
) Each column of R ′ xu (z)= [R(x out , u in , z)
] holds the wavefield received by the virtual transducer at r out for an incident wave-field emitted from a transducer at u in . Reciprocally, each row of R ′ ux (z) = [R(u out , x in , z)] contains the wave-front recorded by the transducers for a virtual source in the focal plane at r in (Fig. 5(a)). Fig. 4(i) shows the phase of R ′ ux (z) obtained at z = 60 mm. At this relatively large depth, the spatial extension ∆u of the reflected wave-field in the transducer basis coincides with the physical aperture A of the array used to collect the echoes coming from a depth z. In contrast, Fig. 4(k) demonstrates that for shallower depths z < A tan[β(r)]/2, ∆u is limited by the numerical aperture of the probe such that ∆u ∼ 2z tan[β max ].

Discussion

We might expect to observe correlations between the columns of matrices R ′ kx (z) and R ′ ux (z) displayed in Figs. 4(c)-(i). Because neighboring virtual sources r in belong a priori to the same isoplanatic patch, the associated wave-fronts observed in the transducer plane or in the far-field should thus be, in principle, strongly correlated since they travel through the same area of the aberrating layer. Such correlations can be revealed by the spatial correlation matrix,

C (R) xx = exp(jarg{R ′ † vx }) × exp(jarg{R ′ vx }) (with v = k or u)
. However, whether it be in the Fourier (Fig. 4(e)) or transducer bases (Fig. 4

(g)), C (R)
xx exhibits a diagonal feature characteristic of uncorrelated wave fronts between the columns of the reflection matrices. In the following, we show how to reveal the hidden correlations in the reflection matrix by introducing the distortion matrix. We will consider mostly the transducer basis, as the far-field case has already been explored in a previous work [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF].

B. The distortion matrix

To reveal the isoplanaticity of the reflected wave-field, each aberrated wave-front contained in the reflection matrix R ′ ux (z) [Fig. 4(a)] should be decomposed into two components: (i ) a geometric component described by Q 0 (z), which contains the ideal wave-front induced by the virtual source r in that would be obtained in the homogeneous medium used to model the wave propagation [see dashed lines in Fig. 4(a)]; (ii ) a distorted component due to the mismatch between the propagation model and reality [Fig. 4(b)]. A key idea is to isolate the latter contribution by subtracting, from the phase of the experimentally measured wave-field, its ideal counterpart. Mathematically, this operation can be done by means of an Hadamard product between R ′ ux (z) and Q * 0 (z):

D ux (z) = R ′ ux (z) • Q * 0 (z). ( 12
)
where the symbol * stands for phase conjugate. We call

D ur = D ux (z) = [D(u out , {x in , z})] the distortion ma- trix.
The distortion matrix connects any input focal point r in to the distorted component of the reflected wavefield in the transducer basis.

To clarify the physical meaning of the distortion matrix, its coefficients can be expressed in the Fraunhoffer approximation as [START_REF] Goodman | Introduction to Fourier Optics[END_REF]:

D(u out , {x in , z}) = δx R ′ (x in + δx, x in , z)e i kc 2z uoutδx , (13)
with δx = x out -x in . Mathematically, each column of D ur is the Fourier transform of the focused wavefield re-centered around each focusing point r in . D ur can thus be seen as a dual reflection matrix for different realizations of virtual sources, all shifted at the origin of the focal plane (x in = 0, see Fig. The matrix operation in ( 12) is equivalent to the classic guide star principle of adaptive focusing [START_REF] Flax | [END_REF]. When considering the reflection matrix R ′ ux between focused and transducer bases in the transmit and receive modes, respectively, the focused transmission is the synthetically generated guide star. The phase conjugate of the propagation matrix Q 0 is the delay compensation one would conventionally use prior to determining the aberration law at the aperture data. The distortion matrix D ur is the conventional delay-corrected aperture data generated by each synthetic guide star.

Note that equivalent distortion matrices, D ru , D kr and D rk can be built from the other reflection matrices previously defined: R ′ xu (z), R ′ kx (z) and R ′ xk (z). For D ru , the same reasoning as above can be used by exchanging input and output. The far-field distortion matrices, D kr and D rk have already been investigated in a previous work [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. The comparison between the phase of R ′ kx (z) [Fig. high degree of correlation of the distorted wave-fronts in the far-field, resulting from the virtual shift of all the input focal spots to the origin. Again this is unambiguously quantified when comparing the corresponding correlation matrix C

C. Local distortion matrices

We have shown that virtual sources which belong to the same isoplanatic patch should give rise to strongly correlated distorted wave-fronts, even if the reflectivity of the medium is random [see Figs. 4(e,f)]. To correct for multiple isoplanatic patches in the field-of-view, recent works [15,[START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF] show that the distortion matrix can be analyzed over the whole field-of-view. Its effective rank is then equal to the number of isoplanatic modes contained in this field-of-view, while its singular vectors yield the corresponding aberration phase laws. The proof-ofconcept of this fundamental result was first demonstrated in optics for specular reflectors [15], then in ultrasound speckle for multi-layered media [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF] and lately in seismic imaging for sparse media [START_REF] Touma | [END_REF].

Here, we investigate the case of ultrasound in vivo imaging, in which fluctuations of the speed-of-sound occur both in the lateral and axial directions. This means that the spatial distribution of aberration effects can become more complex. In Fig. 2(e), strong fluctuations of the F -map illustrate the relatively small size of isoplanatic patches that can only be induced by a complex distribution of the speed-of-sound in the medium. Such complexity implies that any point in the medium will be associated with its own distinct aberration phase law. To construct an image in these conditions, the transmission matrix connecting the correction and focused bases should therefore be constructed to include all of these phase laws -an extremely difficult task. Here, this problem will be tackled using an analysis of a local distortion matrix. The idea is to take advantage of the local isoplanicity of the aberration phase law around each focusing point.

To begin, we divide the field-of-illumination into overlapping regions that are defined by their central midpoint r p and their spatial extension ∆r = {∆x, ∆z}. All of the distorted components associated with focusing points r in located within each region are extracted and stored in a local distortion matrix D ′ ur (r p ):

D ′ (u out , r in , r p ) = D(u out , r in ) W ∆r (r in -r p ), ( 14
)
where W ∆r (r) = 1 for |x| < ∆x and |z| < ∆z, and zero otherwise. Ideally, each sub-distortion matrix should contain a set of focusing points r in belonging to the same isoplanatic patch. In reality, the isoplanicity condition is never completely fulfilled. A delicate compromise thus has to be made on the size ∆r of the window function: it must be small enough to approach the isoplanatic condition, but large enough to encompass a sufficient number of independent realizations of disorder [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. This last point is discussed in Sec. III F.

D. Isoplanicity

The isoplanatic hypothesis can be ensured by means of the adaptive confocal filter described in Sec. II D. The parameter l c can actually control the extension of the aberrated PSF in the filtered matrix R ′

xx . If we model the spatially-distributed aberrations at each depth z by a thin aberrating layer located at z/2, l c should scale as l c ∼ λz/(2l p ), where l p is the coherence length of the aberrator. As l p also governs the isoplanatic length in the focal plane [21], the isoplanatic hypothesis will be fulfilled if ∆x ∼ l p . This, combined with the fact that δx 0 ∼ λz/(N p) in the far-field, gives l c ∼ nδx 0 , where n ∼ N p/(2∆x); this relation dictated the choice of extension l c of the adaptive confocal filter in Table I.

The isoplanatic condition can now be assumed to be fulfilled over each region of size ∆r. This hypothesis implies that the PSFs H in and H out are invariant by translation in each region. This leads us to define a local spatially-invariant PSF H (l) in/out around each central midpoint r p such that:

H in/out (x ′ , x, z) = H (l)
in/out (x ′x, z, r p ). Injecting (3) into (13) leads to the following expression for the D-matrix coefficients:

D ′ (u out , r in , r p ) = H(l) out (u out , r p ) dxγ(x + x in , z)H (l) in (x, r p )e i kc 2z uoutx . (15)
The physical meaning of this last equation is the following: Around each point r p , the aberrations can be modelled by a transmittance H(l) out (u out , r p ). This transmittance is the Fourier transform of the output PSF H (l) out (x, r p ):

H(l) out (u out , r p ) = x H (l) out (x, r p )e -i kc 2z uoutx . (16)
The aberration matrix H(l) out directly provides the true transmission matrix Q between the transducers and any point r p of the medium:

Q = H(l) out • Q 0 . ( 17
)
This transmission matrix Q, or equivalently the aberration matrix H(l) out , are the holy grail for ultrasound imaging since their phase conjugate directly provides the focusing laws that need to be applied on each transducer to optimally focus on each point r p of the medium.

E. Singular value decomposition

To extract the aberration phase law H(l) out (u out , r p ) from each local distortion matrix, we can notice from (15) that each line of D ′ ur (r p ) is the product between H(l) out (u out , r p ) and a random term associated with each virtual source. To unscramble the deterministic term Hout (u out , r p ) from the random virtual source term in (15), the singular value decomposition (SVD) of D ′ ur (r p ) can be applied. The SVD consists in writing D ′ ur (r p ) as

D ′ ur (r p ) = U(r p ) × Σ(r p ) × V † (r p ), (18)
where the symbol † stands for transpose conjugate. Σ is a diagonal matrix containing the singular values σ i (r p ) in descending order: σ 1 > σ 2 > .. > σ N . U(r p ) and V(r p ) are unitary matrices that contain the orthonormal set of output and input eigenvectors,

U i (r p ) = [U i (u out , r p )] and V i (r p ) = [V i (r in , r p )].
The physical meaning of this SVD can be intuitively understood by considering the asymptotic case of a point-like input focusing beam [H

(l) in (x, r p ) = δ(x)].
In this ideal case, (15) becomes D(u out , r in , r p ) = H(l) out (u out , r p )γ(r in ). Comparison with (18) shows that D ′ ur (r p ) is then of rank 1 -the first output singular vector U 1 (r p ) yields the aberration transmittance H(l) out (r p ) while the first input eigenvector V 1 (r p ) directly provides the medium reflectivity.

However, in reality, the input PSF H in is of course far from being point-like. As we will show below, the spectrum of D ′ ur (r p ) displays a continuum of singular values (Fig. 6(b)). The support of U 1 (r p ) is limited by a correlation length δu 1 related to the spatial extent δx in of the input PSF, and V 1 (r p ) only provides a low-resolution image of the medium reflectivity since its spatial frequency content is also limited by δu 1 . Interestingly, the normalized first output singular vector, Û1 (u out , r p ) = U 1 (u out , r p )/|U 1 (u out , r p )|, can still constitute an estimator of H(l) out (r p ) in this case. However, several effects will induce a bias on this estimator. First, the medium should ideally exhibit a random reflectivity; any local structure such as a plane reflector can alter the estimation of the aberration phase law. Second, in this speckle regime, the medium's random reflectivity should be smoothed out enough by considering a spatial window W ∆r containing a large number N in of independent input focal points r in . Finally, just like correlation techniques which rely on the synthetic guide star principle [START_REF] Flax | [END_REF], the spatial extension δx in of the transmitted focal spot degrades the quality of the estimator Û1 (u out , r p ). However, unlike correlation techniques, the SVD bias is of a different nature and can be mitigated for a sufficiently large number N in of disorder realizations. These aspects are now investigated in detail by considering the correlation matrix of D ′ ur (r p ) in the transducer basis.

F. Correlation matrix

To study the SVD of the distortion matrix D ′ ur (r p ) in the transducer basis, the correlation matrix (12), the resulting distortion matrix extracts the aberrated component of those wavefronts. From an other point of view, all the wavefronts are realigned as if they were generated by input focal spots that are virtually shifted at the origin (13). (C) The correlation matrix Cuu of Dur mimics the time reversal operator applied to a virtual reflector that results from a coherent average of all the shifted input focal spots. Applying an incident plane wave to Cuu is the first step of the time reversal process. (d) The resulting wave-field Xout is the estimator of the aberration phase law extracted by conventional correlation techniques. (e) This estimator can be improved by iterating the time reversal process, or equivalently performing the eigenvalue decomposition of Cuu [START_REF] Flax | [END_REF]. The phase conjugate of its first eigenvector U1 then yields the phase law to focus on the brightest part of this coherent reflector. needed:

C uu (r p ) is (a) (b) (e)
C uu (r p ) = N -1 in D ′ ur (r p ) × D ′ † ur (r p ), (19)
with N in the number of virtual sources contained in each spatial window ∆r. The SVD of D ′ (r p ) is indeed equivalent to the eigenvalue decomposition of C uu (r p ):

C uu (r p ) = U(r p ) × Σ 2 (r p ) × U † (r p ), (20)
or, in terms of matrix coefficients,

C uu (r p ) = i σ 2 i (r p )U i (u out , r p )U * i (u out , r p ). ( 21
)
The eigenvalues σ 2 i of C uu (r p ) are the square of the singular values of D ′ ur (r p ). The eigenvectors U i (r p ) of C uu (r p ) are the output singular vectors of D ′ ur (r p ). The study of C uu (r p ) should thus lead to the prediction of singular vectors U i (r p ).

The coefficients of C uu can be seen as an average over r in of the spatial correlation of each distorted wave-field:

C(u, u ′ , r p ) = 1 N in rin D ′ (u, r in , r p ) D ′ * (u ′ , r in , r p ). ( 22
)
In the absence of the adaptive confocal filter defined in (5), each element C(u out , u ′ out , r p ) of C uu (r p ) would correspond the correlation coefficient between the timedelayed complex (IQ) signals recorded at the ballistic time by the receiving elements, u out and u ′ out , averaged over the set of input focusing points r in contained in the spatial window centered on r p [5]. C uu can be decomposed as the sum of a covariance matrix ⟨C uu ⟩ and a perturbation term N:

C uu = ⟨C uu ⟩ + N. (23)
C uu will converge towards ⟨C uu ⟩ if the perturbation term N tends towards zero. In fact, the intensity of N scales as the inverse of the number N in of resolution cells in each sub-region [9,[START_REF] Robert | Evaluation of Green's functions in complex media by decomposition of the Time Reversal Operator: Application to Medical Imaging and aberration correction[END_REF]. In the following, we will thus assume a convergence of C towards its covariance matrix ⟨C⟩ due to disorder self-averaging.

The covariance matrix can be derived analytically in the speckle regime for which the medium reflectivity γ(r) is assumed to be random, meaning that ⟨γ(r)γ * (r ′ )⟩ = ⟨|γ| 2 ⟩δ(r-r ′ ). Under this assumption, injecting (13) into [START_REF] Robert | Evaluation of Green's functions in complex media by decomposition of the Time Reversal Operator: Application to Medical Imaging and aberration correction[END_REF] leads to:

C(u out , u ′ out , r p ) ∝ ⟨|γ| 2 ⟩ Hout (u out , r p ) H(l) out * (u ′ out , r p ) H(l) in * H(l) in (u out -u ′ out , r p ) , (24)
where the symbol * stands for a correlation product. The correlation term, H(l) in * H(l) in , results from the Fourier transform of the input PSF intensity H (l) in 2 . Equation ( 24) is reminiscent of the Van Cittert-Zernike theorem for an aberrating layer [START_REF] Mallart | [END_REF]. This theorem states that the spatial correlation of a random wavefield generated by an incoherent source is equal to the Fourier transform of the intensity distribution of this source (here the input aberrated focal spots).

Let us write the eigenvalue decomposition of the positive semidefinite correlation kernel,

H(l) in * H(l) in = H(l) in * H(l) in (u out -u ′ out ) , H(l) in * H(l) in = W × L × W † (25)
or, in terms of matrix coefficients, (24), one obtains the following expression for C uu :

H(l) in * H(l) in (u out -u ′ out ) = N i=1 l i W i (u out ) W * i (u ′ out ) (26
C uu ∝ H(l) out • W × L × W • H(l) out † (27)
Still under the assumption that the aberrations only induce phase retardation effects (| H(l) out (u out , r p )| = 1), the matrix H(l) out • W is unitary. Because of the uniqueness of the eigenvalue decomposition, Eqs. 20 and 27 show that the eigenvalue distribution of C uu is dictated by its correlation kernel (σ 2 i ∝ l i ) and that its eigenvectors U i are given by:

U i (r p ) ∝ H(l) out (r p ) • W i (r p ). (28)
Thus, the determination of whether the eigenvectors U 1 (r p ) can be satisfactory estimators of H(l) out rests on the properties of W i (r p ), the eigenvectors of H(l) in * H(l) in . W i (r p ) can be derived by solving a second order Fredholm equation with a Hermitian kernel [24,[START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]. An analytical solution can be found for certain analytical form of the correlation function

H(l) in * H(l) in : in the absence of aberration [ H(l) in (u in ) = 1], H(l) in * H(l)
in should be equal to a triangle function that spreads over the whole correlation matrix [9]. In presence of aberration, a significant drop of the correlation width δu in of H(l) in * H(l) in is expected. δu in is actually inversely proportional to the spatial extent δx in of the input PSF H in : δu in ∼ λz/δx in [START_REF] Mallart | [END_REF]. Fig. 6(a) illustrates that fact by showing the modulus of the correlation matrix C uu (r p ) computed over the area A in Fig. 2(c). If we assume that the aberrations only induce phase retardation effects (| H(l) out (u out , r p )| = 1), the modulus of C uu is actually a direct estimator of H(l)

in * H(l) in . As shown by Fig. 6(a), the correlation function H(l) in * H(l) in is far from having a triangular shape and it decreases rapidly with the distance |u out -u ′ out |. For such a bounded correlation function, the effective rank of C uu is shown to scale as the number of resolution cells contained in the input PSF H in [24]:

M δ ∼ (∆u/δu in ) ∼ (δx in /δx 0 ). ( 29
)
The shape of the corresponding eigenvectors W i (r p ) depends on the exact form of the correlation function, or equivalently on the shape of the virtual scatterer. For instance, a flat reflector (sinc correlation function) gives 3D prolate spheroidal eigenfunctions [24]; a cylindrical object leads to Hermite-Gaussian eigenmodes [27]. As the correlation function H(l) in * H(l) in is, in first approximation, real and positive (i.e associated with a symmetric PSF envelope | H(l) in |), a general trend is that the first eigenvector W 1 (r p ) shows a nearly constant phase. The phase of the first eigenvector U 1 (r p ) is then a direct estimator of H(l) out (r p ) [blue continuous line in Fig. 6(d)]. In practice, W 1 (r p ) can exhibit a linear phase ramp due to the asymmetry of the input PSF envelope H (l) in . This results in an additional phase ramp in U 1 (r p ) compared to the targeted aberration phase law H(l) out (u out , r p ). We will provide in Sec. III G a physical interpretation of this problem and a method to circumvent it.

A more critical issue comes from the fact that W 1 (r p ) (and thus |U 1 (r p )|) exhibits a bell curve shape whose characteristic width δu 1 is dictated by the spatial extent of the input PSF: [27]. The higher rank eigenvectors W i (r p ) are more complex and exhibit a number of lobes that scale with their rank i. Fig. 6(c) shows the modulus of the first two eigenvectors of the matrix C uu shown in Fig. 6(a). We recognize the typical signature of the two first eigenmodes with one and two lobes respectively. For aberration correction, it is thus important to consider the normalized vector Û1 = [U 1 (u out )/|U 1 (u out )] that only implies a phase shift, rather than the original vector U 1 as is done in Ref. [START_REF] Bendjador | [END_REF]. In the latter case, the bounded support of U 1 will limit the numerical aperture to δu 1 /z, and, for strong aberrations, deeply degrade the resolution of the corrected image. However, the shape of |U 1 | implies that the aberration phase law estimated from Û1 is inconsistent on the edges of its support.This bias is induced by the perturbation term in [START_REF] Mallart | [END_REF] whose variance

δu 1 ∼ λz 1 + z/δx in ∼ z λ/δx in
scales as |N (u out , u ′ out )| 2 = |C(u out , u out )| 2 /N in .
Taking a perturbational approach, U 1 can be written as the sum of its expected value, W 1 • H(l) out , and a first-order perturbation term δU 1 , with

|δU 1 (u out )| 2 ∼ M 2 δ /(N in N out )
. This decomposition can be seen as the sum of a constant phasor and a weak random phasor that implies the following phase error on the eigenvector 2 1 , with u 0 the center of the U 1 support. The phase error can then be expressed as:

U 1 : ⟨|arg{U 1 }| 2 ⟩ = |δU 1 | 2 /|W 1 | 2 [28]. Assum- ing a cylindrical virtual scatterer [27], |W 1 (u out )| 2 ∼ (M 1 /N out ) exp -(u out -u 0 ) 2 /δu
⟨|arg{U 1 (u out )}| 2 ⟩ ∼ M 2 δ N in M 1 exp u 2 out /δu 2 1 ( 30
)
with M 1 = (∆u/δu 1 ) ∼ √ λδx in /δx 0 . The exponential term in (30) confirms that the SVD bias mainly occurs on the edges of the array in the transducer basis. A correct estimation of the aberration phase law (28) over the whole probe aperture is obtained provided that :

N in > M 2 δ M 1 exp(M 2 1 ). ( 31
)
This condition is quite restrictive and justifies our initial choice for the area ∆r: N in = ∆x∆z/(δx 0 δz 0 ) ∼ 5000 (see Table I) with δz 0 = c/(2∆f ), the axial resolution of the ultrasound image. It also explains why the aberration correction process should then be iterated; at each iteration, the focal spot size, δx in or δx out , decreases and the spatial window ∆r can be reduced accordingly. In the end, the measurement of aberration matrices H(l)

out and H(l) in will thus have excellent spatial resolution.

G. Time reversal picture

In this section, we describe the SVD process using a time reversal picture. This picture allows for an intuitive physical interpretation of the theoretical concepts introduced in the previous section, and highlights the differences between our aberration correction approach and other, more conventional techniques [5,10,[START_REF] Flax | [END_REF]. We also show how the aforementioned additional phase ramp may arise in the estimated aberration phase law, and how to correct for this artifact.

In the previous section, we showed that the aberration phase law can be extracted from the SVD of the distortion matrix D ′ ur (r p ). This operation can be seen as a fictive time reversal experiment. Expressed in the form of ( 24), C uu is analogous to a reflection matrix R associated with a single scatterer of reflectivity |H in (x)| 2 [Fig. 5(c)]. For such an experimental configuration, it has been shown that an iterative time reversal process converges towards a wavefront that focuses perfectly through the heterogeneous medium onto this scatterer [START_REF] Prada | [END_REF]30]. Interestingly, this time-reversal invariant can also be deduced from the first eigenvector of the time-reversal operator RR † [START_REF] Prada | [END_REF][30][31]. The same decomposition could thus be applied to C uu in order to retrieve the wavefront that would perfectly compensate for aberrations and optimally focus on the virtual reflector. This effect is illustrated in Fig. 5(e).

This time reversal picture illuminates the difference between the approach in this article, and the correlation techniques that are widely used for aberration compensation. The latter basically build an estimator X out of the aberration phase law by a simple average of the correlation coefficients [32]:

X out (u out ) = N -1 out u ′ out C(u out , u ′ out )
. This operation is equivalent to the first iteration of the iterative time reversal procedure when a plane wave is first emitted from the transducers to initiate the process (Fig. 5(c)). X out then corresponds to the wave-field reflected by the virtual scatterer (Fig. 5(d). The input PSF intensity |H (l) in (x)| 2 is blurred due to aberration, meaning that the virtual scatterer (or, equivalently, the guide star) has an enlarged spatial extent. An aberration phase law estimated using such a enlarged/distorted guide star will have an inherent bias. Indeed, using (24), the cross-correlation estimator X out can be expressed as follows:

X out (u out ) ∝ H(l) out (u out , r p ) H(l) in * H(l) in ⊛ H(l) out (u out , r p ) , (32)
The estimator X out thus yields the aberration phase law (left term in ( 32)) modulated by a correlated random phase term (right term in ( 32)) with a coherence length scaling with δu in ∼ λz/δx in . The estimator thus exhibits a strong bias due to the input PSF blurring [START_REF] Flax | [END_REF]. This bias can be circumvented by iterating the time reversal process, or equivalently by performing a SVD of the distortion matrix. Indeed, the iteration of the time reversal process tends to maximize the energy backscattered by the virtual reflector [24,27,33]. It yields a wave-front U 1 that tends to focus on the center of the virtual reflector over a resolution length δx 1 ∼ λz/δu . The phase of U 1 can then exhibit an additional linear phase ramp compared to the true aberration phase law [33]: Û1 (u out ) = Hout (u out )e -i kc 2z uoutx0 (33) where x 0 corresponds to the lateral shift of the corresponding PSF H 1 , such that

H 1 (x) = uout Û1 (u out )e i kc 2z uoutx = H 1 (x -x 0 ). ( 34
)
If no effort is made to remove this shift, each selected area defined by the spatial window function W ∆r could suffer from arbitrary lateral shifts x 0 compared to the original image. This artifact can be suppressed by removing the linear phase ramp in (33). One way to do this is to estimate x 0 using the auto-convolution product of the incoherent PSF |H 1 | 2 :

|H 1 | 2 x ⊛ |H 1 | 2 (x) = |H (l) out | 2 x ⊛ |H (l) out | 2 (x -2x 0 ) (35) If a Gaussian covariance model is assumed for aberra- tions, the auto-convolution |H (l) out | 2 x ⊛ |H (l)
out | 2 should be maximum at x = 0. In that case, the maximum position of

|H 1 | 2 x
⊛ |H 1 | 2 (x) leads to an estimation of 2x 0 . Once the latter parameter is known, the undesired linear ramp can be removed from Û1 (33), as illustrated in Fig. 6(d).

In the present case, this linear phase ramp compensation is applied to each estimated aberration phase law. Nevertheless, it should be noted that some specific shapes of aberrating layer -such as a wedge -can manifest as a linear phase ramp in the aberration phase law. For such a particular case, the lateral shift observed on the corresponding PSF is physical and should not be compensated for. The removal of the linear phase ramp should thus be used with caution as it could cancel, in some specific cases, the benefits of the aberration correction process.

H. Transmission matrix estimator

Aberration correction in the receive transducer basis

Now that an estimator Û1 of the aberration matrix H(l) out has been derived, its phase conjugate can be used as a focusing law to compensate for aberrations. We start by correcting in the transducer basis in receive (output). This means building an estimator Q 1 of the transmission matrix from the Hadamard product of the free space transmission matrix Q 0 and the phase conjugate of Û1 :

Q 1 = Q 0 • Û * 1 . ( 36
)
Q 1 is then used to recalculate the broadband FR matrix R xx (z) (1), leading to a corrected FR matrix R (c) xx (z) (Fig. 3(b)):

R (c) xx (z) = Q † 1 (z) × Q 0 (z) × R xx (z). ( 37
)
Note that the correction is applied to the raw FR matrix R xx (z), and not to the filtered FR matrix R ′ xx (z). This is to make sure that no singly-scattered echo is removed during the aberration correction process; as highly-aberrated singly-scattered echoes can extend quite far from the diagonal of the FR matrix, some of this information might be lost after a confocal filter is applied.

Aberration correction in the transmit plane wave basis

Next, aberration correction is performed for the transmit mode (input). It is important to note that it is possible to change the correction basis if need be. Here, correcting in the plane-wave basis in the transmit mode is more efficient since the ultrasound emission sequence was performed in this basis. Indeed, it can also help to compensate for axial movements of the medium that may have occurred during the acquisition. Such movements give rise to additional phase shifts in the illumination basis. A set of dual reflection matrices R ′ xk (z) (8) is built from the updated FR matrix and an aberration phase law is extracted from the corresponding distortion matrices D ′ rk (r p ) [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. The aberration correction process in the plane wave basis is similar to that described here for the transducer basis (36)- (37), replacing the matrix Q 0 by T 0 (7) (Fig. 3(c)). The result is an updated FR matrix R (c) xx (z), an example of which is displayed in Fig. 1(e) at z = 39 mm. The corresponding CMP intensity profile is displayed in Fig. 1(g) for the area B shown in Fig. 2(c). Compared to the initial CMP intensity profile, the result of this first step of the aberration correction process seems quite modest (see Table . II). This is explained by the relatively large size of the spatial window ∆r chosen at the first step of the UMI process (Table. I). While the central part of the FR matrix seems thinner in Fig. 1(e) compared to its initial counterpart in Fig. 1(c), the focusing quality is not drastically improved on the eccentric parts of the field-of-view. These residual aberrations will be tackled in the next iteration of the aberration correction process by: (i ) projecting the ultrasound data in the transducer basis at input in order to compensate for aberrations induced by superficial layers at large depths; (ii ) projecting the ultrasound data in the plane wave basis at output in order to compensate for aberrations at shallow depths; (ii ) reducing the area ∆r in order to address higher-order aberrations associated with smaller isoplanatic patches. The aberration correction process can now be iterated over smaller areas (see Table . I). Note that the correction bases are also exchanged between input and output to minimize any redundancy in the algorithm and optimize the efficiency of this second step. Compared to the previous step, the CMP intensity profile displayed in Fig. 1(g) illustrates the gain both in terms of resolution and contrast of the input and output PSFs. After this second step, the transverse resolution is actually enhanced by a factor of three in the area B and the contrast shows an improvement of almost 3 dB [see Table . II].The transverse resolution is estimated from the full width at half maximum, w -3dB , of the CMP intensity profile. The contrast is here computed as the ratio between the single scattering energy and the incoherent background at focus.

Convergence of the UMI process

The iteration of the aberration correction process can then be pursued since the quality of focus is improved at each step. This results in more highly-resolved virtual transducers which, in return, provides a better estimation of the aberration phase law (in particular at large angles in the plane wave basis or near the edge of the array in the transducer basis). As before, this process can be repeated at input and output in both correction bases, while again reducing the size of the isoplanatic patches and opening the confocal numerical pinhole (Table. I). The final corrected FR matrix R (c) xx (z) is displayed in Fig. 1(e). Comparing this result with the initial R xx (z) [Fig. 1(b)] illustrates the benefit of UMI. Whereas the single scattering contribution originally spread over multiple resolution cells, it now lies along the diagonal of the FR matrix. The part of the back-scattered energy that remains off-diagonal is mainly due to multiple scattering events taking place ahead of the focal plane [11,14]. The corresponding CMP intensity profile is shown in Fig. 1

(f).

A comparison with the initial profile illustrates both the gain in terms of contrast (> 5 dB) and resolution (×3.5) provided by UMI (see Table. II).

Figure 7 displays the depth evolution of the focusing factor F obtained with the UMI process described above (red curve). It confirms the clear improvement of the focusing quality compared to its original value (black curve). It also shows the importance of alternating the correction bases (Fourier and transducer planes) both at input and output. If the aberration correction is only done in the measurement bases (plane wave basis at input and transducer basis at output), the focusing quality is much lower especially at shallow and large depths. Finally, Fig. 7 highlights the effect of the adaptive confocal filter that improves the focusing quality, especially at large depths where the incoherent background dominates.

Nevertheless, note that the final focusing criterion reached by UMI (F ∼ 0.65) remains far from its ideal value (F = 1). In the next section, we will show through numerical simulations that this limit of the UMI process can be accounted for by the continued presence of highorder aberrations. These aberration effects are associated with an isoplanatic length smaller than the final size ∆r of the spatial windows (Table. I). Reducing ∆r even further is not possible, since the condition (31) would be no longer fulfilled, thereby leading to a biased estimation of the aberration phase laws.

IV. RESULTS

A. Numerical validation

To validate our aberration correction method, and to investigate its limits, we now apply it to the numerical simulation described in Ref. [14] (low speckle regime). For this simulation, the ultrasound transmission sequence used 161 steering angles spanning from -40 o to 40 o . The corresponding ultrasound image is displayed in Fig. 8(a). The black, round inclusion near z = 20 mm displays unexpected speckle noise -this noise is caused by aberrations induced by the inhomogeneous speed-of-sound and density distributions in superficial layers of tissues [14]. The F-map displayed in Fig. 8(c) quantifies the aberration level in the ultrasound image. Unsurprisingly, the quality of focus is more degraded at short depths while aberrations are smoothed out at larger depths due to a smaller numerical aperture. This is confirmed by Fig. 8(g) that shows several aberration laws across the field-of-view. The previously described aberration correction process is applied using the parameters given in Tab. I. The resulting ultrasound matrix image is displayed in Fig. 8 ] that is directly extracted from the local input and output PSFs, H in and H out [14]. A slight discrepancy can be observed on top right part of the image and can be explained by the short-scale fluctuations of aberrations at shallow depths (see Fig. 8(e)). As the focusing factor F is not able to grasp the short-scale variations of focusing quality in this area [14], the spatial resolution of the transmission matrix estimator is probably not sufficient to capture the local aberration phase laws in this area.

The limits of UMI can also be highlighted by comparing the aberration laws estimated by the UMI process with the true time delays. In the transducer basis, the latter can actually be measured numerically by placing a source at a point of interest in the medium and recording the transmitted wave-front on the array of transducers. The comparison between the estimated and true time delay laws is shown for different locations in Fig. 8(i). While UMI succeeds in retrieving the low spatial frequency components of the wave-front, its short-scale fluctuations are not captured by the aberration correction process. The main reason for this is that the local isoplanatic condition required by the aberration correction procedure is not satisfied by the high-order aberrations. Indeed, in the transducer basis, the isoplanatic length is directly related the coherence length of the aberrating layer. The spatial resolution of the estimated aberration phase law is thus roughly given by the transverse size ∆ x (here 3 mm, see Tab. I) of the spatial window considered for each local distortion matrix. The aberration law is then a low-pass filtered version of the true time delay law. This explains why F does not reach a uniform and optimal value of 1 over the whole field-of-view after the aberration correction process [Fig. 8(e)]. Note also that the poor estimation of the focusing law on the top right part of the image in Fig. 8(i) is in agreement with the poor focusing quality shown by Fig. 8(f) in the same area.

To provide a more systematic study of UMI performance, we now consider a set of four numerical layers introduced by Mast et al. in their seminal paper [34]. The corresponding speed-of-sound distributions are provided The UMI process is then applied using the parameters given in Tab. I. The final Strehl ratio S F can be assessed from our estimator φ(u, r) of the aberration phase law:

S F (r) = exp i[ϕ(u, r) -φ(u, r)] u 2 .
Its averaged value ⟨S F (r)⟩ r is reported in Tab. III. For each realization, a clear increase of the Strehl ratio is observed compared to its initial value S, the gain being TABLE III. Result of the UMI process on four realizations of aberrating layers extracted from Ref. [34].

Realization #1 #2 #3 #4 ⟨S⟩ 0.49 0.8 0.79 0.76 ⟨S F ⟩ 0.64 0.87 0.89 0.83 higher when S is smaller. For instance, this gain is drastic for realization #1: ⟨S F ⟩/⟨S⟩ ∼ 130 %. Note also that the final Strehl ratio S F never exceeds a value of 0.9. As already highlighted by Fig. 8(i), this is because UMI fails in capturing the short-scale fluctuations of the aberration phase laws.

B. Ultrasound matrix image

The UMI process now having been validated numerically, we can turn our attention to the experimental ultrasound image I M which is built from the diagonal of the corrected FR matrices R (c,d) displays the gallbladder in close-up; UMI reveals a clear view of the internal wall that was completely blurred in the original image. This improvement is a valuable one for in vivo imaging, as gallbladder inflammation commonly manifests as a thickening of the gallbladder wall.

To quantify this improvement, maps of the focusing criterion are shown in Figs. 2(e,f), before and after the UMI process [14]. The focusing parameter F is improved over a large part of the image [Fig. 2(f)]. The gain in focusing quality is particularly spectacular at shallow depths, as is evidenced by the increase in F from F ∼ 0.2 [see Fig. 2(e)] to an optimal value [F ∼ 1, see Fig. 2(e)]. At larger depths, the F-factor is improved on both the left and right sides of the gallbladder but does not reach an optimal value, saturating around F ∼ 0.7. The subsistence of high-order aberrations associated with small isoplanatic lengths can account, at least partially, for this non-ideal focusing quality, even after the UMI process. An inhomogeneous attenuation of the wave-field across its angular spectrum can also explain this saturation of the F -factor, as will be discussed below.

To estimate the benefit of the UMI process in terms of image contrast, the background rate α was estimated before and after the UMI process [14]. The corresponding α-maps are shown in Figs. 2(g,h). The background rate α is drastically decreased over a large part of the image [Fig. 2(h)]. This is particularly spectacular around the gallbladder which explains the much better quality of its close-up image (Fig. 2(b)) compared to the orig-inal one (Fig. 2(b)). The α-map also shows that the contrast is less improved in a triangular region to the left of the image, and at larger depths (z > 40 mm) behind the gall bladder. This effect could be explained by an spatially-varying attenuation across the angular spectrum of the incident wave-field, induced by upstream echogeneous structures. It would also account for the saturation of the F -factor at large depths. Overcoming this issue would require an extension to the UMI process in order to estimate both the amplitude and phase of the transmission matrix. This estimator would enable an inversion procedure of both the attenuation and aberrations induced by the medium rather than just a phase conjugation of wave-front distortions. Fig. 10 shows the spatial distribution of estimated aberration phase laws at the conclusion of the matrix imaging process. The truncated aspect of some of the aberration laws results from the maximal angles of illumination or collection imposed by the finite size of the ultrasonic array. At small depths [z < A tan[β(r)]/2, see Fig. Fig. 10 highlights the spatial variations of the aberration phase across the field-of-view. This spatial resolution was made possible by gradually reducing the size of the spatial window W ∆r at each step of the matrix imaging process (see Table. I). The four colored straight rectangles in Fig. 2(c) depict the size of the spatial windows used at each step. A 75% overlap is applied between subsequent spatial windows in order to retrieve the aberration matrix Hin/out at a high resolution ∆r. The spatial correlations in the aberration matrices can now be quantitatively investigated (Sec. III F).

C. Estimated aberration phase laws

D. Aberration matrices and isoplanatic modes

Matrix imaging also enables the mapping of isoplanatic modes across the field-of-view. As already shown for the distortion matrix with specular reflectors [15], a SVD can be highly useful in extracting the characteristic spatial variations of aberration. To that aim, we define a full aberration matrix from the input and output aberration matrices defined either in the plane wave or transducer bases:

H = Hin Hout . ( 38
)
The input and output aberration phase laws, estimated for each point r p , are stored along each column of H = [ H({v in , v out }, r p )], with v = u or k depending on the considered basis. The SVD of H can be written as follows:

H = p s p A p × I † p . (39
C (H) (r ′ p , r p ) = {vin,vout} H * ({v in , v out }, r ′ p ) H({v in , v out }, r p ).
The first eigenvector I 1 is thus the spatial domain where the degree of correlation between aberration phase laws is maximized. This degree of correlation is quantified by (c) corresponding isoplanatic modes Ip.

the normalized eigenvalue ŝ2

1 , such that

ŝ2 p = s 2 p i s 2 i = I † p × C (H) × I p Tr[C (H) ] . (40)
The corresponding singular vector A 1 = s -1 1 H × I 1 contains the associated input and output aberration phase laws,

A in,1 = [A in,1 (v in )] and A out,1 = [A out,1 (v out )], such that A p = A in,p A out,p .
The same process can be iterated on the matrix Hs 1 A 1 × I † 1 to retrieve the second eigenstate and so on. A set of orthogonal isoplanatic modes I p is finally obtained with a degree of correlation ŝp that decreases with their rank.

Fig. 11(a) displays the corresponding eigenvalues in the plane wave and transducer bases. A few predominant eigenvalues associated with the main isoplanatic modes seem to emerge from a continuum of lower eigenvalues associated with a noise background in each case. To determine the effective number of isoplanatic modes supported by the field-of-view in each basis, one can consider the entropy H(ŝ 2 p ) of the normalized eigenvalues ŝ2

p that yields the effective rank of the matrix C (H) : H(ŝ p ) =i ŝi log 2 (ŝ i ). Here, the entropy is 3.1 and 3.9 in the Fourier and transducer bases, respectively. Fig. 11 shows the three first eigenstates of H in the plane wave and transducer bases. We first remark that the retrieved aberration phase laws, A in,p and A out,p , are not strictly equal at input and output; this is despite the fact that the transmit and back-scattered waves travel through the same heterogeneities, and so the manifestation of aberrations should be identical in input and output. The partial non-reciprocity stems from the different input and output bases used to: (i ) originally record the reflection matrix; (ii ) correct aberrations at each step of the UMI process [see Table . I]. For this latter reason, some distorted components can emerge more clearly in one basis rather than the other, as aberrations in each basis are not fully independent (especially in the far-field).

Interestingly, the measured features of an isoplanatic mode differ according to the correction basis. In particular, the first isoplanatic mode I 1 [Figs. 11(c 1 ,d 1 )] confirms the fact that the plane wave and transducer bases are more effective at small and large depths, respectively [see Sec. IV C]. Moreover, each correction basis addresses aberrations of a different nature. Indeed, the spatial extension of the isoplanatic patches is deeply affected by the correction basis, thereby impacting the result of the aberration correction process. Depending on the location of the aberrating layer and/or its spatial dimension, one basis is therefore more suitable than another to extract the aberration law.

On the one hand, in most in-vivo applications, the organs under study are generally separated from the probe by layers of skin, adipose and/or muscle tissues. In such layered media, aberrations are invariant by translation from the plane wave basis. The plane wave eigenstates displayed in Fig. 11 Interestingly, the overall observation is a correlation between the calculated isoplanatic modes and the tissue architecture revealed by UMI [Fig. 2(d)]. The SVD of the aberration matrices thus provides a segmentation of ultrasound images that could be related to the actual distribution of tissues in the field-of-view. This information would be valuable for the quantitative mapping of important bio-markers such as the speed-of-sound [11,36,37].

V. DISCUSSION AND PERSPECTIVES

In this paper, we have presented a validation of the UMI process for local aberration correction by means of a numerical simulation of aberrations through the abdominal wall, and an in vivo experiment performed on a human gallbladder. Ultrasound imaging is actually often used for diagnosis of cholecystitis [38] (inflammation of the gallbladder), either by detecting gallstones or identifying pericholic fluid and thickened gallbladder wall. This application of in-vivo ultrasound imaging is one of the most challenging, involving various types of tissues and a strongly heterogeneous distribution of the speed-of-sound. Such a medium includes areas of both strong and weak scattering, which is a difficult problem for standard adaptive focusing techniques. Backscattered echoes are generated either by unresolved scatterers (ultrasound speckle), bright point-like scatterers, anechoic regions (gallbladder) or specular structures that are larger than the image resolution (for example, the skin-muscle interface around z = 5 mm). Previous works have shown that the distortion matrix concept can be applied to both specular objects [15], random scattering media [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF], or sparse media made of a few isolated scatterers [START_REF] Touma | [END_REF]. This article details a global strategy for aberration correction that is essential for applying UMI to invivo configurations where all of these scattering regimes can be found.

While our method is inspired by previous works in ultrasound imaging [6,9,[START_REF] Mallart | [END_REF]33,39], it features several distinct and important differences. The first one is its primary building block: The broadband FR matrix that precisely selects all of the singly-scattered echoes originating from each focal point. This is a decisive step since it greatly reduces the contribution of: (i ) out-offocus echoes thanks to ballistic time gating; (i ) multiply-scattered echoes by means of an adaptive confocal filter. Secondly, the matrix approach provides a generalization of the virtual transducer interpretation [9]. By decoupling the location of the input and output focal spots, this approach becomes flexible enough for a local estimation and compensation of aberrations at both input and output and in any correction basis, thereby unifying the different methods developed in the past for aberration correction from a single basis [4-6, 9, 10].

On the one hand, one can easily project the FR matrix at input (or output) to build a D-matrix in the different bases of interest by simply considering monochromatic propagators at the central frequency. On the other hand, the output (or input) of the D-matrix remains in a focused basis which enables a local estimation of aberration phase laws by truncating the field-of-view into limited spatial windows. It thus enables an iterative procedure of aberration correction with alternative changes of correction basis between input and output, while gradually improving the spatial resolution of the transmission matrix estimator. Last but not least, local aberration phase laws are estimated in the Fourier domain, which is much more efficient in terms of computational speed [5] than the temporal cross-correlation of time-delayed signals usually performed for adaptive focusing. It also enables the estimation of the aberration phase laws via a SVD of the distortion matrix [START_REF] Bendjador | [END_REF][START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. Equivalent to an iterative time reversal process [6,9,33], this estimator is more robust than usual correlation techniques [5,10,[START_REF] Flax | [END_REF] with respect to guide-star blurring.

UMI thus constitutes a powerful tool for imaging a heterogeneous medium when little to no previous knowledge on the spatial variations of the speed-of-sound is available. It can be applied whatever the array geometry (curved probes, phased arrays, matrix arrays, etc.) or acquisition scheme (focused excitations, plane wave, diverging waves, etc.). Optimized contrast and resolution can be recovered for any pixel of the ultrasound image. This approach also paves the way towards a mapping of the speed-of-sound distribution inside the medium [11,36,37] by revealing the different isoplanatic modes in the ultrasound image. Such an optimized focusing process is also critical for accurate characterization measurements such as local measurements of ultrasound attenuation [40], scattering anisotropy [41] or the microarchitecture of soft tissues [42].

3 FIG. 1 .

 31 FIG. 1. Focused reflection matrix. (a) UMI consists in splitting the locations of the transmit (rin) and receive (rout) focusing points. (b) Conventional ultrasound image of the gallbladder. (c)-(f) Evolution of Rxx(z) at depth z = 39 mm (white dashed line in (b) during the aberration correction process. Rxx is shown (b) prior to correction, (c) after the application of the adaptive confocal filter (5), (d,e) after the first and fourth iteration of the aberration correction process, respectively. (g) Evolution of the corresponding CMP intensity profiles I(r, ∆x) spatially averaged over the area B [Fig. 2(c)].

FIG. 2 .

 2 FIG. 2. Results of the aberration correction process applied to in-vivo imaging of the gallbladder. (a,b) Zooms on the gallbladder of the initial and corrected UMI images displayed in (c,d), respectively. (c) Conventional dynamically focused (i.e confocal) image. The continous white square A shows the area where the correlation matrix displayed in Fig. 6(a) is computed. The dashed white square B defines the area used to estimate the average CMP intensity profiles shown in Fig. 1(g). The four solid, colored rectangles correspond to the spatial windows W∆r used in the four steps of the correction process (from the largest to smallest, see Table.I). (b) Corrected dynamically focused UMI image. (e,f) Focusing criterion superimposed onto the (e) conventional and (f) corrected images. (g,h) Incoherent background rate superimposed onto the (g) conventional and (h) corrected images.

  FIG. 2. Results of the aberration correction process applied to in-vivo imaging of the gallbladder. (a,b) Zooms on the gallbladder of the initial and corrected UMI images displayed in (c,d), respectively. (c) Conventional dynamically focused (i.e confocal) image. The continous white square A shows the area where the correlation matrix displayed in Fig. 6(a) is computed. The dashed white square B defines the area used to estimate the average CMP intensity profiles shown in Fig. 1(g). The four solid, colored rectangles correspond to the spatial windows W∆r used in the four steps of the correction process (from the largest to smallest, see Table.I). (b) Corrected dynamically focused UMI image. (e,f) Focusing criterion superimposed onto the (e) conventional and (f) corrected images. (g,h) Incoherent background rate superimposed onto the (g) conventional and (h) corrected images.

Focusing ( 1 ),( 2 )

 12 Confocal filter(5)

  5(b)). The co-location of the virtual sources at the same point is the reason why the phase laws along each column of D ux (z) [Figs. 4(j)] are much more correlated than those of R ′ ux (z) [Fig. 4(g)]. This observation is quantitatively confirmed by the corresponding correlation matrix, C (D) xx = exp(jarg{D † ux }) × exp(jarg{D ux }), displayed in Fig. 4(h). It shows much larger off-diagonal correlation coefficients than the original reflection matrix R ′ ux (Fig. 4(g)).

FIG. 4 .

 4 FIG. 4. Revealing the spatial correlations between reflected wave-fields. By subtracting (a) the phase of the dual reflection matrix (sketched by continous lines) by its geometrical counterpart (dashed line), (b) a distortion matrix is computed and reveals the long-range spatial correlations exhibited by the reflected wave-fronts over each isoplanatic patch. (c,d) Phase of the dual reflection matrix, R ′ kx , and distortion matrix, D kx , in the plane wave basis. (e,f) Corresponding spatial correlation matrices of their phase, C (R) xx and C (D) xx , respectively. (g,h) Phase of the dual reflection matrix, R ′ ux , and distortion matrix, Dux, in the transducer basis. (i,j) Corresponding spatial correlation matrices of their phase, C (R) xx and C (D) xx , respectively. The matrices shown in panels (c)-(j) correspond to a depth of z = 60 mm. (k) The white dashed lines in panels (c)-(d) and (g)-(h) account for the finite angular extent β(r) and spatial support ∆u(r) of the reflected waves in the plane wave and transducer bases, respectively.

  (D) xx (Fig. 4(f)) with its original counterpart C (R) xx (Fig. 4(e)).

1 FIG. 5 .

 15 FIG.5. Sketch of the time reversal analysis of the distortion matrix. (a) Each column of the reflection matrix Rur corresponds to the reflected wavefield induced by the associated virtual source r in . (b) By removing the geometrical curvature of each reflected wavefront(12), the resulting distortion matrix extracts the aberrated component of those wavefronts. From an other point of view, all the wavefronts are realigned as if they were generated by input focal spots that are virtually shifted at the origin(13). (C) The correlation matrix Cuu of Dur mimics the time reversal operator applied to a virtual reflector that results from a coherent average of all the shifted input focal spots. Applying an incident plane wave to Cuu is the first step of the time reversal process. (d) The resulting wave-field Xout is the estimator of the aberration phase law extracted by conventional correlation techniques. (e) This estimator can be improved by iterating the time reversal process, or equivalently performing the eigenvalue decomposition of Cuu[START_REF] Flax | [END_REF]. The phase conjugate of its first eigenvector U1 then yields the phase law to focus on the brightest part of this coherent reflector.

FIG. 6 .

 6 FIG. 6. Extraction of the aberration phase law from the correlation matrix in the output transducer basis computed over the area A in Fig. 2(c). (a) Modulus of Cuu. (b) Ten first normalized eigenvalues of the correlation matrices Cuu. (c) Modulus of the two first eigenvectors of Cuu, U1 (blue line) and U2 (red line). (d) The aberration phase law Û1, before (continuous blue line) and after the linear phase ramp correction (dashed blue line).

FIG. 7 .

 7 FIG.7. Depth evolution of the focusing factor F for several aberration correction schemes applied to the gallbladder experiment.

FIG. 8 .

 8 FIG. 8. Results of the aberration correction process applied to the numerical experiment . (a) Conventional dynamicallyfocused image. (b) Corrected UMI image.Estimated (c,d) and true (e,f) F -maps before and after matrix correction of aberrations, respectively. (g,h) Incoherent background rates before and after matrix correction of aberrations, respectively. (i) Examples of aberration laws (red lines) in the transducer basis resulting from the matrix imaging process and compared to their true value (dashed black line).

  (b). Compared to the original image [Fig.8(a)], the corrected image shows a clear contrast improvement in the anechoic inclusion. It can be quantified by the incoherent background rate [Fig.8(g,h)] that is drastically decreased by the UMI process. A contrast improvement of 6 dB is observed in the vicinity of the anechoic inclusion. The transverse resolution is also vastly improved; the F -factor now shows a nearly optimal value (F ∼ 1) throughout the field-of-view [compare Figs.8(c) and (d)]. This gain in focusing quality is confirmed by the ground-truth value of the F -factor [Fig.8(f)

FIG. 9 .

 9 FIG. 9. Speed-of-sound distributions for the different numerical models of the abdominal wall[34] simulated with k-wave.

  (c) xx : I M (r) = |R (c) (r, r)| 2 . Fig. 2 compares the resulting image I M of the gallbladder [Fig. 2(d)] with the original one [Fig. 2(c)]. The two images are normalized by their mean intensity and are displayed over the same dynamic range. A significant improvement of the image quality is observed especially at shallow depth and on the left of the gallbladder. Figs. 2

2 FIG. 10 .

 210 FIG. 10. Examples of estimated aberration phase laws resulting from the matrix imaging process, computed in transmit (blue curves) and in receive (red curves) modes, in the plane wave (1) and transducer (2) bases. The dimension of the selected areas is defined by the size of the spatial window used at the fourth step of the whole process (see Table.1).

  FIG. 10. Examples of estimated aberration phase laws resulting from the matrix imaging process, computed in transmit (blue curves) and in receive (red curves) modes, in the plane wave (1) and transducer (2) bases. The dimension of the selected areas is defined by the size of the spatial window used at the fourth step of the whole process (see Table.1).

  Fig. 10 shows the spatial distribution of estimated aberration phase laws at the conclusion of the matrix imaging process. The truncated aspect of some of the aberration laws results from the maximal angles of illumination or collection imposed by the finite size of the ultrasonic array. At small depths [z < A tan[β(r)]/2, see Fig. 4(k)], the spatial extent of the reflected wave-field is limited by the numerical aperture of the probe [see e.g Fig. 10(b 2 )]. A plane wave basis is thus preferable since its angular range is almost invariant over the focal plane [see Fig. 10(b 1 )]. On the contrary, in the far-field, the angular range of plane waves reaching each focal point is

  ) where s p are the singular values arranged in decreasing order. A p = [A p ({v in , v out })] are the singular vectors in the transducer or Fourier bases, and I p = [I p (r p )] are the singular vectors in the focused basis. For a physical interpretation of these vectors, we take advantage of the equivalence between the SVD of H and the eigenvalue decomposition of the spatial correlation matrix, C (H) = H † × H. The elements of C (H) correspond to the correlation coefficients between aberration phase laws obtained for each image pixel r p and r ′ p :

FIG. 11

 11 FIG. 11. (a) First eigenvalues of the correlation matrix C (H) in the plane wave and transducer bases. (b) Phase of the first three singular vectors Ap (from top to bottom) in the plane wave basis (blue line: Ain, red line: Aout) and (b) corresponding isoplanatic modes Ip. (d) Phase of the first three singular vectors Ap (from top to bottom) in the transducer basis (blue line: Ain, red line: Aout).(c) corresponding isoplanatic modes Ip.

  confirm this statement. While I 2 focuses on the superficial layers of skin and fat [Fig. 11(c 2 )], I 1 mainly spreads between these superficial layers and the gallbladder [Fig. 11(c 1 )]. The corresponding aberration phase laws exhibit a parabolic shape with different curvatures, which is characteristic of a layered medium [Figs. 11(b 1 ,b 2 )]. On other hand, a local perturbation of the medium speed-of-sound located at shallow depth such as superficial veins will have a strong impact on the signals that are measured by transducers located directly above. The transducer basis is then the most adequate for those local variations of the speed-of-sound. I 2 and I 3 focus on the top right and left parts of the ultrasound image, respectively [see Figs. 11(d 2 ,d 3 )]. The corresponding aberration phase laws exhibit fluctuations of high spatial frequencies; such phase laws are characteristic of those which would be induced by local speed-of-sound variations at shallow depths [Fig. 11(e 2 ,e 3 )].

TABLE I .

 I Parameters used for the UMI process. Transmit basis k in u in k in u in Receive basis u out k out u out k out

	Correction steps 1	2	3	4
		l c /δx 0 7.5 10 12 12
	exp.	∆x (mm) 12 9.6 7.7 6.2
		∆z (mm) 16 12 10.2 8.1
		∆x (mm) 10 7	5	3
	num.	∆z (mm) 10 7	5	3

  1 ∼ √ λδx in < δx in . The SVD estimator of the aberration phase law is thus less impacted by the blurring of the input focal spot and the presence of multiple scattering and/or noise than standard correlation techniques based on the guide star principle.However, another issue may occur if the scattering dis-

	tribution H in (x) (l)	2	is too complex. U 1 then focuses on
	the brightest spot exhibited by the input PSF H	in (x) (l)	2

TABLE II .

 II Result of the UMI process in the area B [Fig.2(a)].

	Correction step 0	1	2	3	4
	F	0.31 0.36 0.42 0.46 0.48
	w -3dB (mm) 1.69 1.65 0.59 0.51 0.50
	Contrast (dB) -2.7 -2.4 -0.1 1.7 2.55
	3. Iteration of the UMI process

ACKNOWLEDGMENT

The authors wish to thank Paul Balondrade and Ulysse Najar whose own research works in optics inspired this study, as well as anonymous reviewers whose feedback has helped us to improve the quality of the manuscript.