William Lambert

Justine Robin

Laura A Cobus

Mathias Fink

Alexandre Aubry

Ultrasound Matrix Imaging-Part I: The focused reflection matrix, the F -factor and the role of multiple scattering

This is the first article in a series of two dealing with a matrix approach for aberration quantification and correction in ultrasound imaging. Advanced synthetic beamforming relies on a double focusing operation at transmission and reception on each point of the medium. Ultrasound matrix imaging (UMI) consists in decoupling the location of these transmitted and received focal spots. The response between those virtual transducers form the so-called focused reflection matrix that actually contains much more information than a confocal ultrasound image. In this paper, a time-frequency analysis of this matrix is performed, which highlights the single and multiple scattering contributions as well as the impact of aberrations in the monochromatic and broadband regimes. Interestingly, this analysis enables the measurement of the incoherent input-output point spread function at any pixel of this image. A fitting process enables the quantification of the single scattering, multiple scattering and noise components in the image. From the single scattering contribution, a focusing criterion is defined, and its evolution used to quantify the amount of aberration throughout the ultrasound image. In contrast to the state-of-the-art coherence factor, this new indicator is robust to multiple scattering and electronic noise, thereby providing a contrasted map of the focusing quality at a much better transverse resolution. After a validation of the proof-of-concept based on time-domain simulations, UMI is applied to the in-vivo study of a human calf. Beyond this specific example, UMI opens a new route for speed-of-sound and scattering quantification in ultrasound imaging.

I. INTRODUCTION

To investigate soft tissues in ultrasound imaging, a sequence of incident waves is used to insonify the medium. Inside the medium, the waves encounter short-scale fluctuations of acoustic impedance, generating backscattered echoes that are used to build an ultrasound image. Conventionally, this estimation of the medium reflectivity is performed using the process of delay-and-sum (DAS) beamforming, which relies on a coherent summation of the signals associated with each echo generated by scatterers in the medium. Signals from a particular echo are selected by computing the time-of-flight associated with the forward and return travel paths of the ultrasonic wave between each transducer and the image voxel. From a physical point of view, time delays in transmission are used to concentrate the ultrasound wave on a focal area whose size is ideally only limited by diffraction. Time delays at reception select echoes coming from this excited area. This process falls into the so-called confocal imaging techniques, meaning that, for each point of the image, a double focusing operation is performed.

The critical step of computing the time-of-flight for each insonification and each focal point is achieved in any clinical device by assuming the medium as homogeneous with a constant speed of sound. This assump- * alexandre.aubry@espci.fr tion is necessary in order to achieve the rapidity required for real-time imaging; however, it may not be valid for some configurations in which long-scale fluctuations of the medium speed of sound impact wave propagation [1]. In soft tissues, such fluctuations are around 5%, as the speed of sound typically ranges from 1400 m/s (e.g. fat tissues) to 1650 m/s (e.g. skin, muscle tissues) [START_REF] Duck | Physical properties of tissue: A comprehensive reference book[END_REF]. In such situations, the incident focal spot spreads beyond the diffraction-limited area, the exciting pressure field at this focusing point is reduced, and undesired echoes are generated by surrounding areas. In reception, the application of an incorrect time delay profile mixes echoes which originate from neighboring points in the medium, resulting in a distorted point spread function (PSF) at the output. These aberrating effects can strongly degrade the image resolution and contrast. For highly heterogeneous media, such aberrations may impact the diagnosis of the medical exam or limit the capability to image some organs. A classic example of this effect is liver imaging of difficult-to-image patients. Because the ultrasonic waves must travel through successive layers of skin, fat, and muscle tissue before reaching the liver, both the incident and reflected wave-fronts undergo strong aberrations [START_REF] Hinkelman | [END_REF].

To assess the quality of focus and monitor the convergence of aberration correction techniques, Mallart and Fink [4] introduced a focusing parameter C, known as the coherence factor in the ultrasound community [START_REF] Hollman | IEEE Ultrasonics Symposium. Proceedings. International Symposium[END_REF]. Based on the van Cittert-Zernike theorem, this indicator is linked to the spatial coherence of the back-scattered field measured by the probe for a focused transmit beam. This method assumes that images are dominated by ultrasonic speckle -a grainy, noise-like texture which is generated by a medium of random reflectivity. This assumption is in fact often valid in medical ultrasound, where soft tissues are composed of randomly distributed scatterers that are unresolved at ultrasonic frequencies. In that case, each transmitted focusing wave excites an ensemble of scatterers which are randomly distributed within the focal spot. The back-scattered wavefront results from the random superposition of the echoes generated by each of those scatterers. If the focus is perfect, the time delays used to focus (in transmission and in reception) are exactly the same as the actual roundtrip time-of-flight of the waves which arrive at, and are backscattered by, scatterers located within the resolution cell. The spatial coherence of the received signal is then maximal; the associated coherence factor C tends towards 2/3 in the speckle regime and in the absence of clutter (contributions from multiple scattering to the detected backscattered signal). A decrease in the quality of focus is matched by a decrease in C; thus, this indicator is popular in the literature for the evaluation of ultrasonic focusing quality.

Very recently, ultrasound matrix imaging (UMI) has been proposed for a quantitative mapping of aberrations in ultrasound imaging [START_REF] Lambert | [END_REF]. Experimentally, the first step in this approach is to record the reflection matrix associated with the imaged medium. This matrix contains the response of the medium recorded by each transducer of the probe, for a set of illuminations. Depending on the problem one is facing, this matrix can be investigated in different bases. Here, for imaging purposes, the reflection matrix will be projected into a focused basis. In contrast with conventional ultrasound imaging that relies on confocal beamforming, the idea here is to apply independent focused beamforming procedures at the input and output of the reflection matrix. This process yields a focused reflection (FR) matrix that contains the responses between virtual transducers synthesized directly in the medium from the transmitted and received focal spots. The FR matrix holds much more information on the medium than a conventional ultrasound image. Importantly, it can be leveraged to discriminate between the single scattering, multiple scattering and noise contributions. From the former, a new focusing criterion F is introduced to assess the local focusing quality, notably in the speckle regime, for any pixel of the image. While this indicator shows some similarities with the coherence factor introduced by Mallart and Fink [4], the parameter F constitutes a much more sensitive and spatially-resolved probe of the focusing quality. In the quest for local correction of distributed aberrations [START_REF] Ali | IEEE International Ultrasonics Symposium (IUS)[END_REF][START_REF] Jaeger | [END_REF][START_REF] Rau | IEEE Int. Ultrason. Symp[END_REF][START_REF] Bendjador | [END_REF][START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF], this parameter F can be of particular interest as it can play the role of a guide star for any pixel of the ultrasound image. More generally, the FR matrix is a key operator for imaging applications [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Badon | [END_REF][13][14][15]. The accompanying article [16] will present in details a matrix method of aberration correc-tion in which the FR matrix plays a pivotal role.

The current paper describes the mathematical construction and physical meaning of the FR matrix. In particular, it is shown how the FR matrix can be exploited to probe the local ultrasound focusing quality, independently of the medium reflectivity. The FR matrix was originally introduced in Ref. [START_REF] Lambert | [END_REF]; here, it is investigated in more detail, and its properties are demonstrated by means of time domain simulations and applied to invivo ultrasound imaging of the calf of a human subject. First, a time-frequency analysis of the FR matrix is performed, and its mathematical expression is derived rigorously for the monochromatic and broadband regimes. In the broadband regime, the analytical expression of the FR matrix goes beyond the result of Ref. [START_REF] Lambert | [END_REF] that consisted in a simple expansion of the monochromatic result at the central frequency. In this article, we discuss the relative weight of single and multiple scattering contributions and the impact of aberrations for each regime. While Ref. [START_REF] Lambert | [END_REF] only considered virtual transducers belonging to the same transverse plane, here we investigate the impulse responses between focusing points located at any depth [see Fig. 1(a)]. A local incoherent inputoutput PSF can thus be extracted from the FR matrix along both the axial and transverse directions. Going beyond Ref. [START_REF] Lambert | [END_REF], a novel fitting process of this PSF is also proposed to discriminate between the single and multiple scattering components of the reflected wave-field. This process leads to a new definition of the focusing factor F as a scaling factor in the transverse direction of the inputoutput PSF compared to its diffraction-limited (ideal) counterpart. (Note that in the latter quantity, the local frequency spectrum of ultrasonic data is taken into account). This entire process is first validated by means of numerical simulations. The comparison with the state-ofthe-art C-map shows a drastic gain in terms of contrast and resolution for the local evaluation of the focusing quality. The gain in contrast is due to the fact that, unlike the coherence factor C, F is robust to both multiple scattering and electronic noise. The gain in transverse resolution is provided by the input-output focusing operation inherent to UMI while the C-factor is estimated from a single focusing operation. We then illustrate this property of F by applying UMI to in-vivo ultrasound data. A predominant incoherent background due to multiple scattering and electronic noise yields a very weak and badly contrasted coherence factor C, while the focusing factor F highlights different regions in which the focus is almost perfect (F ∼ 1), and others in which the focusing process is degraded by irregular layers of muscle, fat and skin at shallow depths.

The paper is structured as follows: Section II describes the experimental procedure to acquire the reflection matrix and the numerical simulations that will be used in the paper. In Sec. III, the FR matrix is investigated in the monochromatic regime. A 2D local input-output PSF is extracted from the FR matrix. The manifestation of aberrations, out-of-focus echoes and multiple scatter-ing in this monochromatic PSF is discussed by means of numerical simulations and experimental in-vivo data. In Sec. V, a coherent sum of the monochromatic FR matrices leads to a broadband FR matrix. This operation, which is equivalent to a time-gating process, suppresses contributions from out-of-focus reflectors. In addition, as shown by numerical simulations, this summing process enables a clear discrimination between the single and multiple scattering components of the reflected wave-field. In Sec. VI, a time-frequency analysis of the FR matrix is performed and shows, in particular, the impact of absorption and scattering on the experimental FR matrix as a function of depth and frequency. In Sec. VII, a fitting procedure is proposed to extract the relative weights of single scattering, multiple scattering and electronic noise. From the single scattering component, a focusing parameter F is obtained by rescaling the transverse evolution of the broadband input-output PSF with the ideal, diffraction-limited one. The overall fitting procedure is validated by means of the time domain simulations and applied to the in-vivo ultrasound imaging of a human calf. Finally, Sec. VIII outlines the advantages of F compared to the coherence factor C widely used in the literature [4,[START_REF] Hollman | IEEE Ultrasonics Symposium. Proceedings. International Symposium[END_REF]. Section IX presents conclusions and general perspectives.

II. METHODS

A. Experimental Procedure

UMI begins with an experimental recording of the reflection matrix, R. In principle, this measurement can be achieved using any type of illumination (element-byelement [17], focused beams [18], etc.); here, a planewave acquisition sequence [19] has been arbitrarily chosen. The probe was placed in direct contact with the calf of a healthy human volunteer, orthogonally to the muscular fibers. (This study is in conformation with the declaration of Helsinki). The acquisition was performed using a medical ultrafast ultrasound scanner (Aixplorer Mach-30, Supersonic Imagine, Aix-en-Provence, France) driving a 5 -18 MHz linear transducer array containing 192 transducers with a pitch p = 0.2 mm (SL18-5, Supersonic Imagine). The acquisition sequence consisted of transmission of 101 steering angles spanning from -25 o to 25 o , calculated for the hypothesis of a tissue speed of sound of c 0 = 1580 m/s [START_REF] Duck | Physical properties of tissue: A comprehensive reference book[END_REF]. The pulse repetition frequency was set at 1000 Hz. The emitted signal was a sinusoidal burst lasting for three half periods of the central frequency f c = 7.5 MHz. For each excitation, the backscattered signal was recorded by the 192 transducers of the probe over a time length ∆t = 80 µs at sampling frequency f s = 40 MHz. Acquired in this way, the reflection matrix is denoted R uθ (t) ≡ R(u out , θ in , t), where u out defines the transverse position of the receiving transducer, θ in , the incident angle of the emitted plane wave and t, the time-of-flight. In the following, the subscripts in and In a monochromatic regime, the synthesized virtual transducer displays an elongated shape in the z-direction because of diffraction. In the broadband domain, the axial resolution is inversely proportional to the signal bandwidth, giving a much thinner virtual transducer in the z-direction. (b)-(c) Monochromatic (f = 7.5 MHz) and broadband common mid-point intensity profile (5) averaged over a set of common mid-points r contained in the white rectangle displayed in Fig. 3(a1) and (a2). Both profiles have been normalized by their maximum.

out will denote the transmitting and receiving steps of the reflection matrix recording process.

B. Numerical Simulations

To validate our aproach, a multistatic synthetic aperture dataset has also been computed with k-Wave [20], a time domain simulation software based on the k-space pseudospectral method. The density and speed-of-sound distributions of the simulated medium [Figs. 2(a) and(b)] mimic the aberrations undergone by ultrasound through the abdominal wall [21]. They result from a spatial lowpass filtering of the numerical tissue layers introduced by Mast et al. [22] [Figs. 2(c) and (d)] on which speckle has been superimposed. Those short-scale fluctuations of these mechanical properties generate a random wavefield characteristic of ultrasound imaging in soft tissues. A 2.6 mm-diameter anechoic disk is also included at a depth of 20 mm to highlight the impact of aberrations and multiple scattering on the imaging process. To investigate the effect of multiple scattering, two k-wave simulations have been performed to mimic ultrasound imaging in low and high speckle scattering regimes. The speckle strength is controlled by means of the standard deviation of the speed-of-sound and of the density, σ c and σ ρ , respectively. The values considered in our simulations are reported in Tab. I. In the Rayleigh regime, the impact of density fluctuations is negligible compared to speed-of-sound variations [23]. Under a scalar, Gaussian and bi-dimensional model of disorder, the scattering mean free path ℓ s can be estimated through the following relation [START_REF] Aubry | Approche matricielle de l'opérateur de propagation des ondes ultrasonores en milieu diffusant aléatoire[END_REF]:

ℓ s ∼ 2(σ c /c 0 ) -2 k -3 0 ℓ -2 c
, with k 0 the wave number and ℓ c the coherence length of speckle. The latter parameter is here roughly equal to the spatial sampling used for the simulation. It yields an order of magnitude for ℓ s that is reported in Tab. I for each simulation.

Each k-Wave simulation is performed over a twodimensional grid. All simulation parameters such as the transducer configuration, transmit pulse, computational grid, and sampling are described in Tab. I. The recording of the reflection matrix is performed as described above for the experiment. The acquisition sequence consisted of transmission of 101 steering angles spanning from -25 o to 25 o . The emitted signal was a 2.25 MHz two-cycle sinusoidal burst. For each excitation, the back-scattered signal was recorded by the 155 transducers of the probe at sampling frequency f s = 80 MHz.

III. MONOCHROMATIC FOCUSED REFLECTION MATRIX

The first post-processing step of UMI consists in projecting the reflection matrix into a focused basis. This step is greatly simplified by performing beamforming operations in the frequency domain, i.e. applying appropriate phase shifts to all frequency components of the received signals in order to realign them at each focal point. A matrix formalism is particularly suitable for this operation, because, in the frequency domain, the projection of data from the plane-wave or transducer bases to any focal plane can be achieved with a simple matrix product previously described in Ref. [START_REF] Lambert | [END_REF]. Note that the focused basis corresponds to a set of vectors that are not linearly independent; this is therefore not, strictly speaking, a basis in a mathematical sense. In contrast, the incident plane waves used to record the reflection matrix are linearly independent, but they do not form a complete basis. Nevertheless, for convenience, we will refer to these two sets of wave-fronts as focused and plane wave bases in the following.

The result of this beamforming process is a set of focused reflection matrices R rr (f ) obtained at each frequency f of the bandwidth. In contrast to standard synthetic ultrasound imaging, in which input and output focusing points coincide, the approach presented here decouples these points. In emission, the incident energy is concentrated at the focusing point r in ; this point can thus be seen as a virtual source [START_REF] Passmann | [END_REF]26]. Similarly, in reception, a virtual sensor is synthesized by selecting the part of the back-scattered wave-field originating from the vicinity of point r out . Therefore, each coefficient R(r out , r in , f ) of R rr (f ) contains the monochromatic responses of the medium between a set of virtual transducers. The position of each virtual transducer maps onto a pixel of the ultrasound image. Fig. 1(a) illustrates this matrix focusing process. Note that, for clarity, the input focusing operation in Fig. 1(a) is represented by a cylindrical wave-front instead of a combination of plane waves. This is justified by the fact that a plane wave synthetic beamforming numerically mimics a focused excitation [19]. Note that the concept of virtual transducers is merely didactic and that, of course, they do not act as real source or sink of energy. Moreover, they are strongly directive: in the downward direction for the virtual source, in the upward direction for the virtual receiver. At last, in the monochromatic regime, the coupling between the transmit and receive beams can also occur above and below the focal depth due to spurious echos from out-of-focus scatterers.

Fig. 3(a 2 ) shows the

x-projection R xx (z, f ) of R rr (f ) at depth z = 18 mm and frequency f = 6 MHz in the experiment. R xx (z, f ) = [R(x out , x in , f, z)]
contains the responses between virtual transducers, r in = (x in , z) and r out = (x out , z), located at the same depth. Note that ∆x has to be limited to avoid the spatial aliasing induced by the incompleteness of the plane wave illumination basis; its maximal value ∆x max is inversely proportional to the angular step δθ of the plane wave illumination basis: ∆x max ∼ λ max /(2δθ), with λ the wavelength. Thus, the coefficients R(x out , x in , f, z) associated with a transverse distance |x out -x in | larger than a superior bound ∆x max are not displayed.

Fig. 3(a 2 ) shows that most of the signal in R xx (f ) tends to concentrate around its diagonal. This feature is characteristic of a predominant single scattering contribution [START_REF] Lambert | [END_REF]. Indeed, singly-scattered echoes mainly originate from a virtual detector r out which is located near to virtual source r in . This is confirmed by Fig. 4 that compares two monochromatic FR matrices R xx (f ) (z = 13 mm, f = 2.25 MHz) obtained numerically in the low and high speckle scattering regimes [Fig. 4(d 1 ) and (d 3 ), respectively]. While the multiple scattering rate is supposed to differ in each case, the FR matrices display a similar feature, which corroborates that single scattering is here predominant.

The diagonal elements of R rr (f ), which obey r in = r out , directly provide the image I(r, f ) which would be obtained via multi-focus (a.k.a. confocal) imaging at frequency f :

I (r, f ) ≡ |R (r, r, f )| 2 . ( 1
)
An example of such an image is displayed in Fig. 3(a 1 ).

Compared to a standard ultrasound image built from broadband signals (Fig. 3b 1 ), this monochromatic image displays poor axial resolution and is therefore difficult to exploit. This observation is confirmed numerically whether it be in the low (Fig. 4(b 1 )) or high (Fig. 4(b 3 ) speckle regime. However, as will be shown in the following, the off-diagonal elements of R rr provide valuable information on the physical properties of the medium, as well as on the wave focusing quality.

Under the Born approximation and its matrix formalism [27], R rr (f ) can be expressed theoretically as follows [START_REF] Lambert | [END_REF]:

R rr (f ) = H ⊤ out (f ) × Γ(f ) × H in (f ), (2)
where the matrix Γ describes the scattering process in the sample; in the single scattering regime, Γ is diagonal and its elements correspond to the medium reflectivity γ(r, f ) at frequency f . H in (f ) and H out (f ) are the input and output focusing matrices, respectively [START_REF] Lambert | [END_REF]. Each column of H in (f ) = [H in (r, r in , f )] corresponds to one focused illumination, and contains the monochromatic PSF at emission -the spatial amplitude distribution of the input focal spot resulting from that particular illumination. Similarly, each column of

H out (f ) = [H out (r, r out , f )]
contains the spatial amplitude distribution of an output focal spot. As will be shown in the following, the spatial extension of these focal spots is a direct manifestation of the mismatch between the actual velocity distribution in the medium and the spatially-invariant velocity model used in the matrix beamforming process. For lighter notation in the rest of this section, the frequency dependence f of each physical quantity is made implicit.

Equation ( 2) can be rewritten in terms of matrix coefficients as follows: R(r out , r in ) = dr H out (r, r out ) γ(r) H in (r, r in ). [START_REF] Hinkelman | [END_REF] This equation shows that each pixel of the ultrasound image (diagonal elements of R rr ) results from a convolution between the sample reflectivity γ(r) and an imaging PSF, H(r, r ′ ), which is itself a product of the input and output PSFs: H(r, r ′ ) = H in (r, r ′ ) × H out (r, r ′ ). As H in (r, r in ) and H out (r, r out ) define the characteristic size of the input focal spot at r in and the output focal spot at r out , these matrices also define the resolution of the resulting confocal image. In the absence of aberration, the transverse and axial dimensions of these focal spots, δx 0 (r) and δz 0 (r), are only limited by diffraction [START_REF] Born | Principles of optics[END_REF]:

δx 0 (r) = λ 2 sin[β(r)] , δz 0 (r) = 2λ sin 2 [β(r)] , (4)
with β(r) the maximum angle by which each focal point is illuminated or seen by the array of transducers. In ultrasound imaging, the radiation pattern of the transducers usually limits the numerical aperture to ∼25 • . The focal spots thus typically display a characteristic elongated shape in the z-direction [δx 0 << δz 0 , see Fig.

IV. COMMON MID-POINT INTENSITY

The off-diagonal points in R rr can be exploited for a quantification of the focusing quality at any pixel of the ultrasound image. To that aim, the relevant observable is the intensity profile along each anti-diagonal of R rr [START_REF] Lambert | [END_REF]:

I(r m , ∆r) = |R(r m + ∆r/2, r m -∆r/2)| 2 . ( 5
)
All pairs of points on a given anti-diagonal have the same midpoint r m = (r out + r in )/2 , but varying spacing ∆r = (r out -r in ). In the following, I(r m , ∆r) is thus referred to as the common-midpoint (CMP) intensity profile. To express this quantity theoretically, we first make an isoplanatic approximation in the vicinity of each CMP r m . This means that waves which focus in this region are assumed to have travelled through approximately the same areas of the medium, thereby undergoing identical phase distortions [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]16]. The input and/or output PSFs can then be considered to be spatially invariant within this local region. Mathematically, this means that, in the vicinity of each common mid-point r m , the spatial distribution of the input or output PSFs, H in/out (r, r in/out ), only depends on the relative distance between the point r and the focusing point r in/out . This leads us to define a local spatially-invariant PSF H x in [mm] -10 0 10

x (mm) common mid-point r m such that:

H in/out (r, r in/out ) = H (l)
in/out (r -r in/out , r m ). ( 6)

The range of validity of this approximation is given by the size of each isoplanatic patch. In our numerical simulation, azimuthal and axial isoplanatic patch sizes are of the order of 3 mm and 5 mm, respectively, at 80% correlation. For the experimental data (acquired on a human calf), we obtain similar values for isoplanatic patch size -these were extracted from the aberration phase laws given by UMI [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF] We next make the assumption that, as is often the case in ultrasound imaging, scattering is due to a random distribution of unresolved scatterers. Such a speckle scattering regime can be modelled by a random reflectivity:

⟨γ(r 1 )γ * (r 2 )⟩ = ⟨|γ| 2 ⟩δ(r 2 -r 1 ), (7)
where ⟨...⟩ denotes an ensemble average and δ is the Dirac distribution. By injecting (3) into (5), then reducing the implicit double integral over r by the δ-function in [START_REF] Ali | IEEE International Ultrasonics Symposium (IUS)[END_REF], and finally using [START_REF] Lambert | [END_REF] and the change of variables (r-r m → r), the following expression can be found for the CMP intensity:

I(r m , ∆r) = dr |H (l) out (r -∆r/2, r m ) | 2 |H (l) in (r + ∆r/2, r m ) | 2 |γ(r + r m )| 2 . ( 8
)
To smooth the intensity fluctuations due to the random reflectivity, a spatial average over a few resolution cells is required while keeping a satisfactory spatial resolution.

To do so, a normalized spatially averaged intensity profile I av (r m , ∆r) is computed in the vicinity of each point r m , such that

I av (r m , ∆r) = ⟨W L (r -r m )I(r, ∆r)⟩ r ⟨W L (r -r m )I(r, 0)⟩ r (9)
where the symbol ⟨...⟩ r denotes the spatial average and W L (r) is a spatial window function, such that

W L (r) = 1 for ∥r∥ < L/2 0 otherwise. ( 10
)
This spatial averaging process leads to the replacement of |γ(r)| 2 in the last equation by its ensemble average ⟨|γ| 2 ⟩. I av (r m , ∆r) then directly provides the convolution between the incoherent input and output PSFs,

|H (l) in | 2 and |H (l) out | 2 : I av (r m , ∆r) ∝ |H (l) in | 2 ∆r ⊛ |H (l) out | 2 (∆r, r m ), (11)
where the symbol ∆r ⊛ denotes a spatial convolution over ∆r. Note that this formula only holds in the speckle regime; for a specular reflector, the CMP intensity profile is equivalent to the intensity of the coherent input-output [START_REF] Lambert | [END_REF]. In either case, this quantity gives two interesting pieces of information. Firstly, the CMP intensity at ∆r = 0 is proportional to the confocal energy at focus. Secondly, the spatial extension of the CMP is linked to the lateral and axial dimensions of the input and output PSFs. Fig. 1(b) displays an example of a two dimensional CMP intensity profile. The corresponding cross-section of this profile at ∆z = 0 is displayed in Fig. 3(a 3 ). This cross-section has been averaged over a set of CMPs contained in the white rectangle of Fig. 3(a 1 ). Similarly to the input/output PSFs H out | 2 displays a cigarlike shape. However, while its axial FWHM is close to the diffraction limit (δz 0 ∼ 2.2 mm, with β = 25 • ), its transverse FWHM (∼ 1 mm) is far from being ideal (δx 0 ∼ 0.25 mm).

PSF, |H (l) in ∆r ⊛ H (l) out | 2
Ideally, the CMP intensity profile would be given by

I 0 (∆x, r m ) = |H (l) 0 | 2 ∆x ⊛ |H (l) 0 | 2 (∆x, r m ), (12)
where H

(l) 0 is the diffraction-limited PSF, such that

H (l) 0 (∆x, r m ) = sinc (2π∆x sin [β(r m )] /λ) (13)
with sinc(x) = sin(x)/x. The CMP intensity profile should thus ideally display a transverse full width at half maximum (FWHM) roughly equal to the diffractionlimited resolution δx 0 .

The poor transverse resolution highlighted by Fig. 1(b) can be explained by several potential effects. First, longscale variations of the speed-of-sound can give rise to aberrations that distort the input and output PSFs [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF]. Second, if there are spurious echoes from out-of-focus scatterers (located above or below the focal plane), then the expansion of the resulting input/output beams will be greater than those originating exactly at the focal plane [Fig. 4(a)]. Note also that short-scale heterogeneities may induce multiple scattering events that give rise to an incoherent background in the CMP intensity profile [START_REF] Lambert | [END_REF]. As we will see, the prevalence of multiple scattering events can be estimated from the amount of signal at off-diagonal elements -an incoherent background which is a combination of multiple scattering contributions and, in the experiment, of electronic noise.

All these contributions constitute a problem for imaging; in the monochromatic regime under examination here, it is extremely difficult to discriminate between the effects of aberration, multiple scattering, and singlyscattered echoes taking place out-of-focus. In the next section, we show that the contributions from out-of-focus echos can be greatly reduced via a time-gating operation. This process enables the recovery of the standard axial resolution exhibited by ultrasound images.

V. BROADBAND REFLECTION MATRIX

Under the matrix formalism, time gating can be performed by building a broadband FR matrix R rr . In the following, we show that besides improving the axial resolution and contrast of the ultrasound image, R rr allows a clear distinction between the contributions from single and multiple scattering.

In the frequency domain, the FR matrix is built by dephasing each RF signal such that scattering paths whose first and last scattering events take place at r in and r out constructively interfere. A coherent sum over the overall bandwidth ∆f can then be performed to build a broadband FR matrix:

R rr (∆f ) = 1 ∆f f+ f- df R rr (f ) (14)
with f ± = f c ± ∆f /2 and f c the central frequency of the RF signal bandwidth. One row of the broadband FR matrix corresponds to the situation in which the transmitted waves are focused at r in , creating a virtual source while the virtual detector probes the spatial spreading of this virtual source at the expected ballistic time, i.e at time t = 0 in the focused basis [Fig. 1(a)]. Thus, the sum of monochromatic FR matrices over the entire bandwidth can be interpreted as a time-gating operation in which echos originating from a certain range of times-of-flight are extracted.

With the time-gating applied ( 14), the axial resolution of the virtual transducers should be drastically improved [see Fig. 1(a)]. To prove this assertion, we now derive an expression for the broadband FR matrix within the formalism of this work. For sake of simplicity and analytical tractability, paraxial and isoplanatic (6) approximations are made. The monochromatic PSFs can be decomposed as follows:

H in/out (r, r ′ , f ) = H in/out (r -r ′ , r m , f ) e j2πf (z-z ′ )/c
(15) where H in/out represents the envelope of the PSF. Injecting (3) and ( 15) into ( 14) leads to the following expression for the coefficients of R rr (∆f ):

R(r out , r in , ∆f ) = dr e j2πfc(2z-zin-zout)/c sinc π∆f c (2z -z in -z out ) H out (r -r out , r m ) γ(r) H in (r -r in , r m ). ( 16
)
where we have assumed, in first approximation, that H in/out is constant over the frequency bandwidth. The occurrence of the sinc factor in the integrand of the last equation accounts for the time gating operation of (14). Its spatial extent yields the expected axial resolution δz 0 ∼ c/2∆f in a broadband regime. Not surprisingly, the coherent sum of ( 14) drastically improves the axial resolution and contrast of the image I(r, ∆f ) built from the diagonals of R rr (∆f ) (14), revealing the micro-architecture of the calf tissues in the experiment (Fig. 3 ], the single scattering contribution along the diagonal of R rr is enhanced with respect to the off-diagonal coefficients. This enhancement is due to the time-gating procedure; contributions from scatterers which sit above and below the focal plane have been eliminated, so the remaining singly-scattered echoes are located near the diagonal of R rr .

Interestingly, numerical simulations show that the offdiagonal energy is now mainly due to multiple scattering events taking place at depths which are shallower than the focal plane. While the monochromatic FR matrices did show a similar off-diagonal energy whatever the speckle level [see Figs. 4(d 1 ,d 3 )], the broadband one now exhibits a larger off-diagonal energy in the high speckle regime [see the comparison between Figs. 4(d 2 ) and (d 4 )]. This observation is quantitatively confirmed by investigating the broadband CMP intensity profile in Fig. 4(c). It now displays a confocal, steep peak mainly due to single scattering on top of a lower but wider pedestal linked to multiple scattering. The relative amplitude of the multiple scattering background with respect to single scattering is larger in the high speckle regime. In Sec. VII, we will show how to map this relative amplitude quantitatively, and demonstrate its sensitivity to speckle level.

Experimentally, the improvement in transverse resolution of the single scattering contribution can also be seen in the broadband CMP intensity profile displayed in Fig. 1(c), and in its cross-section shown in Fig. 3(b 3 ). Nevertheless, although we are in a broadband regime, the 2D focal spot in Fig. 1(c) still exhibits a cigar-like shape. To illustrate why this is expected,, we now express the ensemble average of the CMP intensity profile (15) in the broadband regime under the paraxial approximation [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF]:

I av (r m , ∆r, ∆f ) = A ∆f f+ f- df |H in | 2 ∆x ⊛ |H out | 2 (∆r, r m , f ), ( 17
)
where A is a constant. Equation (17) shows that, unlike traditional speckle statistics [START_REF] Wagner | [END_REF], the CMP intensity profile essentially eliminates the phase components of the wave-field and yields a sum over the frequency bandwidth of the incoherent input-output PSF

|H in | 2 ∆x ⊛ |H out | 2 .
This explains why the axial resolution in Fig. 1(c) has not been improved by the time gating process, as was the case for the broadband FR matrix (Eq. 16). Nevertheless, the CMP intensity evolution along ∆x still offers a way to estimate the transverse resolution of the imaging PSF in the broadband regime. Indeed, Fig. 3(b3) shows that the transverse FWHM of I av (r m , ∆r, ∆f ) (∼ 0.5 mm) remains far from optimal (δx 0 (f c ) ∼ 0.25 mm).

In the next section, we aim to develop a better way to evaluate the aberration and multiple scattering levels in ultrasound imaging. We will define quantitative parameters to measure the focusing quality and multiple scattering rate at any pixel of the ultrasound image. The focusing criterion will correspond to the scaling factor that maps the transverse evolution of the CMP intensity profile, I av (r m , ∆r, ∆f ), onto its ideal diffraction-limited counterpart, Ī0 (r m , ∆r, ∆f ). A local multiple scattering rate will be also deduced from this fitting process. To make our measurement quantitative, the theoretical prediction of Ī0 (r m , ∆r, ∆f ) should be as accurate as possible; in the following, we work towards this accuracy by developing a theoretical time-frequency analysis of the FR matrix.

VI. TIME-FREQUENCY ANALYSIS OF THE FOCUSED REFLECTION MATRIX

A time-frequency analysis of the FR matrix is required to investigate the evolution of absorption and scattering as a function of frequency. To do so, the coherent sum of the monochromatic FR matrices [Eq. 14] can be performed over a smaller bandwidth δf centered on a given frequency f :

R rr (f, δf ) = 1 δf f +δf f -δf df ′ R rr (f ′ ) ( 18
)
We have shown that the axial dimension δz 0 of the virtual transducers is inversely proportional to the frequency bandwidth δf . Thus, a compromise must be made between the spectral and axial resolutions. Here, the following choice has been made: δf = 0.5 MHz and δz 0 = 3 mm. ; nevertheless, the time-frequency analysis of the FR matrices yields the evolution of the SNR versus depth and frequency. At z = 18 mm, for instance, the FR matrix at f = 9.75 MHz exhibits a tiny confocal intensity enhancement on top of a predominant noise background (SNR∼3dB). Conversely, the FR matrices at f = 5.25 and 7.5 MHz exhibit a CMP intensity profile which more closely resembles its broadband counterpart. This weak SNR at 9.75 MHz can be partially explained by the finite bandwidth of the transducers (5 -10 MHz). Absorption and scattering losses undergone by ultrasonic waves in soft tissues also have a strong impact on the ultrasound image. Fig. 3(d) illustrates this effect by displaying the spectrum of the confocal signal, I(r, f, δf ) x , as a func-tion of depth. This spectrum shifts towards low frequencies as a function of depth. This frequency shift seems characteristic of absorption losses in soft tissues [START_REF] Duck | Physical properties of tissue: A comprehensive reference book[END_REF]. However it is difficult to discriminate between absorption and scattering losses [31] since the scattering mean free path ℓ s can also vary across the frequency bandwidth. Anyway, whether due to absorption or scattering losses, the variation of the temporal frequency spectrum of backscattered echoes has a strong impact on the local resolution of the ultrasound images. In the next section, we show how to incorporate this frequency dependence in the theoretical expression of the ideal CMP intensity profile, in order to establish a precise and quantitative focusing factor.

VII. THE LOCAL FOCUSING FACTOR AND MULTIPLE SCATTERING RATE

In this section, a local focusing criterion is established for the broadband ultrasound image. For the sake of lighter notation, the dependence of each physical quantity with respect to ∆f is omitted. Aberrations caused by medium heterogeneities degrade the resolution of the ultrasound image and induce a spreading of singlyscattered echoes over the off-diagonal coefficients of R rr . In the speckle regime, it is difficult to determine by eye whether the image is aberrated, and if so, which areas are the most impacted. Interestingly, the CMP intensity profile can yield an unambiguous answer to this question. In the speckle regime, this profile yields the convolution between the incoherent input-output PSF averaged over the frequency bandwidth (Eq. 17). While the incoherent input-output PSF is not exactly equal to the confocal imaging PSF (Eq. 16) it nevertheless fully captures the impact of transverse aberrations. It thus constitutes a relevant observable for assessing focusing quality.

To measure a local focusing factor and multiple scattering rate, the normalized CMP intensity profile should be first decomposed as the sum of three contributions:

I av (r m , ∆x) = α S (r m )I 0 (r m , F (r m )∆x) + α M (r m )I M (r m , ∆x) + α N (r m ), (19)
where α S (r m ), α M (r m ), and α N (r m ), are the respective rates of local single scattering, multiple scattering, and noise. These three quantities obey the following relation:

α S (r m ) + α M (r m ) + α N (r m ) = 1.
The first term accounts for the single scattering contribution and assumes that, on average, aberrations reduce the effective numerical aperture by a focusing factor F . The spatial dependence of F is thus a re-scaled version of the diffraction-limited CMP profile I 0 (r m , ∆x) that would be obtained in absence of aberrations. The second term in Eq. 19 accounts for the multiple scattering contribution. In first approximation, its spatial profile is assumed to be Gaussian [32], such that I M (r m , ∆x) = exp -∆x 2 /(2σ M (r m ) 2 ) , with σ M the spatial extent of the diffuse halo. At last, the third term accounts for the noise background that is constant with respect to ∆x on average.

In the following, the free parameters F , α S , α M and σ M will be determined through a fitting procedure, but we must first accurately determine the reference CMP profile I 0 (r m , ∆x). To do so, the frequency spectrum of the ultrasound image should be taken into account. For each CMP r m , I(r m , f ) is an estimation of the desired frequency spectrum [Fig. 3(d)]. I 0 (r m , ∆x) can then be computed by performing a frequency average of the theoretical CMI intensity profile, I 0 (r m , ∆x, f ) [START_REF] Badon | [END_REF], weighted by the spectrum I(r m , f ) of the ultrasound image. Now that Ī0 (r m , ∆x) is properly estimated, the CMP intensity profile is fit with Eq. 19 using the Matlab nonlinear programming solver fiminsearchcon. This fitting procedure is first validated by the numerical simulations. Fig. 5(a) shows the result of this fit in the low speckle regime at an arbitrary depth z = 33 mm. While the single scattering peak seems to be nicely fit by a rescaled version of I 0 (|∆x| < 2 mm), a larger fit error is observed for the multiple scattering background (|∆x| > 2 mm). Two main reasons can account for this discrepancy: (i ) the Gaussian model used to describe the transverse evolution of multiple scattering contribution that is only strictly valid in the diffusive limit [32], and (ii ) a lack of statistical average for smoothing the fluctuations induced by the random reflectivity of the medium. Indeed, a compromise has to be made between the size L of the spatial window W L (10) used to average the CMP intensity profile and the spatial resolution of the fitting parameters in (19). This size L should be small enough to preserve the spatial variations of aberrations across the field-of-view, but large enough to smooth out the ultrasound speckle. Indeed, the speckle fluctuations in Īav (r m , ∆x) decrease as 1/N in with N in ∼ L 2 /δx 0 δz 0 ), the number of resolution cells contained in each spatial window. Here we chose L = 3.5λ c (with λ c the wavelength at the central frequency) for the numerical data which yields N ∼ 45. Fig. 5(b) shows the spatial evolution of the F -factor across the field-of-view. To evaluate the performance of the fitting procedure, the F -map is compared to its ground truth value (Fig. 5(c)). This ground truth map is built using a local measure of the input and output PSFs, which is possible in numerical simulations. To do so, each PSF H in/out (r, r in/out , t) is computed by simulating the propagation of an incident wave-field designed to focus at r in/out from the plane wave (input) or transducer (output) basis. The PSF corresponds to the transmitted wave-field recorded in the focal plane with a time origin at the expected ballistic time. A temporal Fourier transform and a change a variable yields the corresponding de-scanned PSFs, H

in and H (l) out , in the frequency domain. The ground truth CMP intensity profile is then provided by the following expression of I av [START_REF] Lambert | Matrix approach for ultrasound imaging and quantification[END_REF]:

I av (r m , ∆x, ∆f ) = B ∆f 2 f+ f- df f+ f- df ′ dx H (l) out (x -∆x/2, r m , f ) H (l) * out (x -∆x/2, r m , f ′ ) H (l) in (x + ∆x/2, r m , f ) H (l) * in (x + ∆x/2, r m , f ′ ) . ( 20
)
where B is a constant. This last expression is more general than (17) since the frequency-dependence of the PSF is here considered. Fig. 5 displays a satisfactory agreement between the F -factor and its ground truth value. Their slight discrepancy can be explained by several reasons. First, as with any ultrasound image, the F is shown as a function of an effective depth z = c 0 t/2 that scales with the time-of-flight t and the reference speed-of-sound c 0 . On the contrary, the reference map displayed in Fig. 5c evaluates the focusing factor at true depths in the medium. Second, the F -factor varies quite rapidly across the field-of-view for z < 20 mm. This limited isoplanicity also contributes to the difference observed between the measured F-factor and its ground-truth value, which is especially evident at shallow depths (z < 20 mm). At larger depths, the slight discrepancy observed around x = -10 mm between the measured F-factor and its true value might be explained by slight imprecision at certain regions of the discrimination between single and multiple scattering components. Again, the statistical fluctuations of the multiple scattering component and our hypothesis that the single scattering component can be described by a re-scaled version of I 0 (r m , F (r m )∆x) may explain, at least partially, the discrepancy observed between the measured F -map and its ground truth.

At shallow depths (z < 20mm), the multiple scattering rate extracted by our fitting process shows a clear difference between the low (Fig. 6(c)) and high (Fig. 6(d This property highlights the fact that our approach enables discrimination between aberration effects (quantified by F ) and multiple scattering phenomena (quantified by α M ). Note that, at larger depths, the invariance of F between the two scattering regimes (Eq. 19) is no longer strictly obeyed. This effect (the difference between Figs. 6(a) and (b) at z > 25 mm) may be explained by a different effective wave velocity of the medium in the high speckle regime.

Beyond short-scale fluctuations of speed-of-sound and density in the medium, multiple scattering can also be caused by the strong impedance mismatch between tissues and connecting tissues in the abdominal wall. To illustrate this effect, a last numerical simulation considering the density and speed-of-sound distributions of Ref. [22] (Figs. 2(c) and (d)) has been performed. The resulting maps of the focusing factor F and multiple scattering rate α M are displayed in Figs. 6(c) and (f). Compared with the initial simulation that considers smoother boundaries between tissues (Figs. 6(a) and (b)), F and α M are both higher in the speckle area behind the abdominal layer. This result can be understood as follows: Stronger impedance mismatch in the abdominal layer gives rise to more intense multiple reflection events, while smoothed variations of the speed-of-sound only deviate the trajectory of the incident and reflected waves, giving rise to more aberrations and less multiple scattering. Fig. 7 shows the application of the fitting process to the experimental data. In Fig. 7(a), the original ultrasound image is shown. The extension L of the spatial window W L has been set to 1.5 mm. Figs. 7(b-d) show separate maps of the factors which affect the image quality. Unlike in numerical simulations, electronic noise is present and a noise rate α N should be quantified. The evolution of α N is shown in Fig. 7(b). Unsurprisingly, a predominant noise occurs at large depth (z > 50 mm) and in areas where the medium reflectivity is weak. In any case, if a medical diagnosis is desired, the features of the ultrasound image in these areas should be more carefully inter- preted. However, this is not the only factor which affects the quality of the ultrasound image (Fig. 7(a)). Clutter is also far from being negligible at shallow depths due to multiple reflection events between superficial layers and at larger depths because of multiple scattering processes in speckle (Fig. 7(c)). With regards to the single scattering contribution (Fig. 7(d)), the focusing quality shows strong variations across the field-of-view, particularly between the left and right parts of image. While high values of F (blue areas) indicate good image reliability, low values of F (green/yellow areas) indicate a poor quality of focus. As expected, yellow areas seem to correspond to blurring of the ultrasound image [Fig. 7(a)]. Gray areas correspond to the situation in which the estimation of the image resolution has failed because of a low SNR (α M + α N > 0.75). In these areas, the single scattering contribution is drowned out by a predominant incoherent background either caused by multiple scattering processes and/or electronic noise.

The ultrasound image shows different structures that are associated with their own speed of sound: (i) mus-cles tissues with three different fiber orientations [areas I, II, III on Fig. 7(a)]; (ii) two veins located at {x, z} = 12, 5 mm and -5, 33 mm; and (iii) the fibula, located at the bottom left of the figure -12, 45 mm. Fig. 7(d) reveals a poor focusing quality at shallow depth that could be explained by the near-field effects induced by the discrete sampling of the array. At larger depths, the F -map exhibits a strong contrast between the right part of the image where the focusing factor F is close to 1 and its left part where the focusing is degraded by the thicker and curved arrangement of superficial layers of skin, fat and muscle in area III.

However, note there cannot be a direct correlation between the B-mode image and the F -and α M -maps. Aberrations and multiple scattering are actually induced by wave velocity inhomogeneities and scatterers along the trajectories of the incident and reflected waves. For a given focusing point, the values of F and α M are thus not related to the tissue architecture revealed by the Bmode image at this point but by speed-of-sound fluctuations and scattering properties of the medium above that point. Retrieving the spatial distribution of the speed-ofsound throughout the medium from the quality of focus provided by the F -factor is a difficult task that is beyond the scope of this paper.

VIII. DISCUSSION

The study presented here provides new insights into the construction of the FR matrix and a refined definition of the focusing criterion and multiple scattering rate. Our results are representative of in-vivo ultrasound imaging in which the medium under investigation, a human calf, is composed of different kinds of tissues, which can themselves be heterogeneous. As the medium is composed of a mix of unresolved scatterers and specular reflectors, the scattering of ultrasound varies considerably in space, ranging from areas of strong to weak scattering. Moreover, ultrasound imaging of muscles is a relevant tool for monitoring neuro-muscular diseases [33]. In that respect, a sharp quantification of aberrations and multiple scattering is crucial for mapping the speed-of-sound [34] and scattering anisotropy [35] in muscles. In shear wave elastography [36], these phenomena can also have a detrimental impact on the inversion procedure that consists in extracting the muscle stiffness from ultrasound movies of shear wave propagation [37].

A direct and important application of this work would thus be to aid aberration correction by employing the parameter F as a virtual guide star for adaptive focusing techniques. Currently, in the literature in this area, this role is performed by the coherence factor C [38][39][START_REF] Dahl | IEEE International Ultrasonics Symposium (IUS)[END_REF]. Here, to highlight the benefit of our matrix approach with respect to the state-of-the-art, we build a map of the standard coherence factor C for the ultrasound image of the human calf. C is equal to the ratio of the coherent intensity to the incoherent intensity of the realigned reflected wave-fronts recorded by the probe for each input focusing beam [START_REF] Robert | [END_REF]. Just as with the calculation of the CMP intensity profile ( 9), the raw coherence factor C(r in ) is then spatially averaged over overlapping spatial windows to smooth the fluctuations due to the random reflectivity. The result is mapped onto the ultrasound image in Fig. 7e. Compared to the F -map (Fig. 7d), the coherence factor C provides a weakly contrasted image of the focusing quality. To understand this difference, the analytical expression of C can be derived under the isoplanatic hypothesis [START_REF] Lambert | [END_REF]. While the CMP intensity enables a full retrieval of the spatial convolution between the incoherent input and output PSFs [START_REF] Lambert | Proc. Natl. Acad. Sci. USA[END_REF], the coherence factor C only probes this quantity at ∆x = 0 [4,[START_REF] Robert | [END_REF][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]:

C(r) = |H in | 2 ∆x ⊛ |H out | 2 (∆x = 0, r). (21)
In the speckle regime, the coherence factor C ranges from 0, for strong aberrations, to 2/3 in the ideal case [4]. However, the presence of multiple scattering and noise can strongly hamper its measurement. As highlighted by Figs. 7(b,c), the incoherent input-output PSF, |H in | 2 ∆x ⊛ |H out | 2 , exhibits an incoherent background that is far from being negligible. Unlike the C-factor, the CMP intensity enables a clear discrimination between the single scattering contribution and the noise components. The F -factor is sensitive only to aberrations, since it is estimated after removing the incoherent background. This crucial feature accounts for the difference in behavior between F and C in Fig. 7, especially at large depths. While F shows a close to ideal value of 1 on the right part of the image beyond z = 10 mm (Fig. 7 Beyond the C-factor, the CMP intensity profile is an interesting analogue of the speckle autocorrelation (SAC) function, where the autocorrelation of the ultrasound image yields the pulse-echo PSF, Hin × Hout , convolved against itself in the single scattering regime [START_REF] Wagner | [END_REF]. Together, the CMP profile and SAC function seem to describe pieces of an even more general convolution over both ∆r and ∆r m , with the CMP intensity profile being the case where ∆r = 0 and SAC function being the case where ∆r m = 0 (i.e. common input and output focusing points). Both are related measures of focusing quality in speckle: the CMP profile yields the proposed F criterion, whereas the SAC FWHM is the traditional spatial resolution [START_REF] Wagner | [END_REF]. Nevertheless, multiple scattering and/or electronic noise give rise to a peaked SAC function, potentially leading to the false impression of a diffraction-limited focus. In contrast, those contributions manifest themselves as an incoherent background in the CMP intensity profile which allows them to be discriminated from the single scattering component. The CMP approach is thus much more robust to multiple scattering and noise; moreover, it is not restricted to the speckle regime but can also be used to probe the focusing quality in presence of specular reflectors [START_REF] Lambert | [END_REF].

In the second paper of the series [16], we present a matrix approach to aberration correction for in vivo ultrasound imaging of a gallbladder. An important aspect of that study is the use of F to map the the image at each step of the aberration correction process. Again, numerical simulations are used to validate the whole process and the measured F -factor is shown to be in agreement with its ground truth value.

IX. CONCLUSION

In summary, we successfully applied the concept of the FR matrix to the case of in vivo ultrasound imaging of a human calf. This approach has been priorly optimized and validated by means of numerical simulations. Thanks to the intuitive concept of virtual transducers, the FR matrix provides a wealth of information on the medium that goes well beyond a single confocal image.

By splitting the location of the transmitted and received focal spots, the local resolution of the ultrasound image can be assessed at any pixel. By performing a timefrequency analysis of the reflection matrix, the contributions of single and multiple scattering and their impact on the resolution and contrast were carefully investigated. This time-frequency study of the FR matrix paves the way towards a quantitative characterization of soft tissues by measuring parameters such as the attenuation coefficient or the scattering mean free path. In the accompanying paper [16], the FR matrix will be used as a key building block of UMI for a local aberration correction. Relatedly, a focusing criterion was defined from the FR matrix in order to quantify the impact of aberrations in the vicinity of each pixel of the ultrasound image. Compared to the coherence factor generally used in the literature [4,[START_REF] Hollman | IEEE Ultrasonics Symposium. Proceedings. International Symposium[END_REF], our focusing parameter is much more robust to noise and multiple scattering and locally maps the focusing quality with a much better transverse resolution. Our focusing parameter is thus promising for use in medical imaging as a reliability index of the ultrasound image. Such an index would be of great use in the development and evaluation of advanced ultrasound imaging modes, and more importantly could help clinicians in their diagnostic. A quantitative and local assessment of image quality would indeed add strength to any ultrasound-based diagnostic, potentially avoiding any further and more invasive examination to the patient. Last but not least, the focusing factor can also be used as a guide star for adaptive focusing techniques, or as a local aberration indicator for UMI [16].

FIG. 1 .

 1 FIG.1. Principle of the FR matrix. (a) UMI consists in splitting the location of the transmitted (rin) and received focusing (rout) points both in the axial and transverse directions, thereby synthesizing virtual transducers that can act as source and detector, respectively at any point in the medium. In a monochromatic regime, the synthesized virtual transducer displays an elongated shape in the z-direction because of diffraction. In the broadband domain, the axial resolution is inversely proportional to the signal bandwidth, giving a much thinner virtual transducer in the z-direction. (b)-(c) Monochromatic (f = 7.5 MHz) and broadband common mid-point intensity profile (5) averaged over a set of common mid-points r contained in the white rectangle displayed in Fig.3(a1) and (a2). Both profiles have been normalized by their maximum.

FIG. 2 .

 2 FIG. 2. Acoustic properties of low speckle media simulated with k-wave. (a,b) Density and speed-of sound distributions resulting from the low-pass filtering of original maps (c,d) introduced by Mast et al. [22]. The speckle statistics is provided in Tab. I.

  coefficient 0.15 -0.4 dB.MHz -1.5 .cm -1

  4(a)], which accounts for the poor axial resolution exhibited by the monochromatic images obtained experimentally [Fig. 3(a 1 )] and numerically [Figs. 4(b 1 ,b 3 )].

FIG. 3 .

 3 FIG. 3. UMI in a monochromatic, broadband and multi-frequency regime: Experimental result. (a) Monochromatic regime (f =7.5 MHz): (a1) Multi-focus image (dB-scale), (a2) reflection matrix Rxx(z) at z = 18 mm [see white dashed line in (a)], and (a3) common mid-point intensity profile I(z, ∆x) ((5)) averaged at the same depth. (b) Broadband regime (5-10 MHz): (b1) Multi-focus image (dB-scale), (b2) Broadband reflection matrix Rxx(z) at z = 18 mm, and (b3) common mid-point intensity profile I(r, ∆x) averaged at the same depth. (c) Multi frequency regime (5-5.5 MHz, 7.25-7.75 MHz, 9.5-10 MHz, from left to right): (c1) Multi-focus images (dB-scale), (c2) multi-frequency reflection matrices Rxx(z) at z = 18 mm, and (c3) common mid-point intensity profiles I(r, ∆x) averaged at the same depth. (d) Spectrogram of the confocal signals versus depth z extracted from the diagonal elements of Rxx(z) (linear scale). The confocal spectrum is normalized by its maximum at each depth. The white rectangle in (a1) and (a2) accounts for the position of CMPs considered for the computation of the CMP intensity profiles displayed in Fig. 1(b) and (c), respectively.

FIG. 4 .

 4 FIG. 4. UMI in low and high speckle scattering regimes: Numerical results. (a) x-cross section of the FR matrix: The transmitted (rin) and received focusing (rout) points are at the same depth. (b) Multi-focus images (dB-scale) in the low (b1, b2) and high (b3, b4) speckle regimes built over a restricted (2.245-2.255 MHz) and broad (0.5-4 MHz) bandwidth, respectively. (c) Correspond CMP intensity profile I(z, ∆x) ((5)) at z = 13 mm [see white dashed line in (b)] in the monochromatic (dashed lines) and broadband (continuous lines) regimes for a weak (blue lines) and high (red lines) speckle. (d) FR matrix Rxx(z) at the same depth in the low (d1, d2) and high (d3, d4) speckle regimes built over a restricted (2.245-2.255 MHz) and broad (0.5-4 MHz) bandwidth, respectively.

(b 1 )

 1 ) and the anechoic region in the numerical simulations (Figs. 4(b 2 ,b 4 ). Fig. 3(b 2 ) shows the cross-section R xx (∆f, z) of the experimental FR matrix R rr (∆f ) at depth z = 18 mm [dotted white line of Fig. 3(b 1 )]. Compared to its monochromatic counterpart [Fig. 3(a 2 )

Fig. 3 (

 3 c) shows the experimental ultrasound images [Fig. 3(c 1 )], FR matrices [Fig. 3(c2)] and CMP profiles [Fig. 3(c3)] for three different frequency bandwidths: 5-5.5 MHz, 7.25-7.75 MHz, 9.5-10 MHz. The axial resolution in each ultrasound image is of course deteriorated compared to the broadband image [Fig. 3(b 1 )]

FIG. 5 .

 5 FIG. 5. Numerical validation of the fitting process on the CMP intensity profile by Eq. 19. (a) Result of the fit at z = 33 mm in the low speckle regime: the fit parameters are F = 0.96, αM = 5.5% and σM = 3.4 mm. (b) Map of the focusing factor F (rm) extracted from the CMP intensity profile averaged over spatial windows of size L = 4 mm. (c) Ground-truth map of F (rm) computed from the measured PSFs Hin and Hout (17).

  )) speckle scattering regimes. On the contrary, the focusing factor F remains remarkably constant in the same depth range [see the comparison between Figs. 6(a) and (b)].

FIG. 6 .

 6 FIG. 6. Mapping of local focusing quality and multiple scattering in numerical simulations. Maps of the (a,b,c) F -factor, (d,e,f) Multiple scattering rate and (g,h,i) C-factor in the low and high speckle regimes of the first simulation (Figs. 2(a) and (b)) and in the low speckle regime of the second simulation (Figs. 2(c) and (d)), respectively.

FIG. 7 .

 7 FIG. 7. Mapping of local focusing quality and multiple scattering in the calf imaging experiment. (a) Ultrasound image of the calf. (b) Noise rate αN (rm). (c) Multiple scattering rate αM (rm). (d) F -factor. (d) C-factor. In (b,c) yellow areas correspond to a weak SNR and high multiple scattering rate, respectively. In(d,e) blue and yellow areas correspond to a high and low quality of focus respectively. In (d), gray areas highlights location where the estimation of the focusing criterion is not reliable because the incoherent background is too high (αN + αM > 0.75)

  Figs.7(b,c), the incoherent input-output PSF, |H in | 2 ∆x ⊛ |H out | 2 , exhibits an incoherent background that is far from being negligible. Unlike the C-factor, the CMP intensity enables a clear discrimination between the single scattering contribution and the noise components. The F -factor is sensitive only to aberrations, since it is estimated after removing the incoherent background. This crucial feature accounts for the difference in behavior between F and C in Fig.7, especially at large depths. While F shows a close to ideal value of 1 on the right part of the image beyond z = 10 mm (Fig.7(d)), C exhibits a low value everywhere (Fig.7(e)) because of a

TABLE I .

 I Simulation parameters.

ACKNOWLEDGMENT

The authors wish to thank Victor Barolle, Amaury Badon and Thibaud Blondel whose own research works in optics and seismology inspired this study, as well as anonymous reviewers whose feedback has helped us to improve the quality of the manuscript.