Reflection matrix approach for quantitative imaging of scattering media

We present a physically intuitive matrix approach for wave imaging and characterization in scattering media. The experimental proof-of-concept is performed with ultrasonic waves, but this approach can be applied to any field of wave physics for which multi-element technology is available. The concept is that focused beamforming enables the synthesis, in transmit and receive, of an array of virtual transducers which map the entire medium to be imaged. The inter-element responses of this virtual array form a focused reflection matrix from which spatial maps of various characteristics of the propagating wave can be retrieved. Here we demonstrate: (i) a local focusing criterion that enables the imaging quality to be evaluated everywhere inside the medium, including in random speckle; (ii) a tomographic measurement of wave velocity, which allows for aberration corrections in the original image; (iii) an highly resolved spatial mapping of the prevalence of multiple scattering, which constitutes a new and unique contrast for ultrasonic imaging. More generally, this matrix approach opens an original and powerful route for quantitative imaging in wave physics.


I. INTRODUCTION
In wave imaging, we aim to characterize an unknown environment by actively illuminating a region and recording the reflected waves. Inhomogeneities generate backscattered echoes that can be used to image the local reflectivity of the medium. This is the principle of, for example, ultrasound imaging [1], optical coherence tomography for light [2], radar for electromagnetic waves [3] or reflection seismology in geophysics [4]. This approach, however, rests on the assumption of a homogeneous medium between the probe and target. Large-scale fluctuations of the wave velocity in the medium can result in wavefront distortion (aberration) and a loss of resolution in the subsequent reflectivity image. Smaller-scale inhomogeneities with high concentration and/or scattering strength can induce multiple scattering events which can strongly degrade image contrast. In the past, numerous methods such as adaptive focusing have been implemented to correct for these fundamental issues in reflectivity imaging [5][6][7][8]. However, such methods are largely ineffective in situations where the focus quality inside the medium can not be determined. An extremely common example of this situation for ultrasound imaging is the presence of speckle, the signal resulting from an incoherent sum of echoes due to randomly distributed unresolved scatterers. Speckle often dominates medical ultrasound images, making adaptive focusing difficult.
Alternately, one can try to exploit effects which are detrimental to reflectivity imaging (such as distortion and scattering) to create different imaging modalities. In the ballistic regime (where single scattering dominates), the refractive index can be estimated by analyzing the distortion undergone by the wave as it passes * alexandre.aubry@espci.fr through the medium. This is the principle of quantitative phase imaging [9] in optics and computed tomography [10] in ultrasound imaging. These approaches, however, require a transmission configuration, which is not practical for thick scattering media, and which is impossible for most in-vivo or in-situ applications in which only one side of the medium is accessible. In the multiple scattering regime, optical diffuse tomography [11] is a well-established technique to build a map of transport parameters. However, the spatial resolution of the resulting image is poor as it scales with imaging depth. Moreover, this approach assumes a purely multiple scattering medium, which does not exist in practice. A single scattering contribution always exists, and is furthermore typically predominant in ultrasound imaging. To evaluate the validity of such images, a local (spatially-resolved) multiple scattering rate would be a valuable observable, but is not accessible with state-of-the-art methods.
Recently, a reflection matrix approach to wave imaging was developed with the goals of: (i ) processing the huge amount of data that can now be recorded with multielement arrays [12][13][14], and (ii ) optimizing aberration correction [15][16][17][18] and multiple scattering removal [19][20][21][22] in post-processing. Such matrix approaches provide access to much more information than is available with conventional imaging techniques. Their recent successes suggest that access to detailed information on aberration and multiple scattering could be capitalized upon for more accurate characterization of strongly heterogeneous media. In this paper, we introduce a universal and non-invasive matrix approach for new quantitative imaging modes in reflection.
Our method is based on the projection of the reflection matrix into a focused basis [22,23]. This focused reflection matrix can be thought of as a matrix of impulse responses between virtual transducers located inside the medium [16]. These virtual transducers are created via numerical simulation of wave focusing, i.e. combining all of the backscattered echoes in such as way as to mimic focusing at a set of focal points that spans the entire medium, in both transmit and receive. While each pixel of a confocal image is associated with the same virtual transducer at emission and reception, the focused reflection matrix also contains the cross-talks between each pixel of the image, and thus holds much more information than a conventional image. Importantly, this matrix allows to probe the input-output point spread function (PSF) in the vicinity of each pixel even in speckle. A local PSF in reflection is a particularly relevant observable since (i ) it allows the quantification of the contribution of aberration and multiple scattering to the image, (ii ) it can be used as a guide star for wave velocity tomography in the medium, thereby correcting the aberrations of the confocal image, and (iii ) it paves the way towards spatially resolved measurements of wave transport parameters, and hence a local characterization of the medium micro-architecture. In this paper, we concentrate on two of the fundamental quantities which can be accessed with our method: a local focusing quality, and a spatially resolved multiple scattering rate. These quantities are crucial since they can lead to estimation of the wave velocity [24] and of transport parameters [20,25,26] such as the absorption length and the scattering mean free path (the mean distance between two successive scattering events). Not only are these parameters quantitative markers for biomedical diagnosis in ultrasound imaging [27][28][29][30][31][32] and optical microscopy [11], but they are also important observables for non-destructive evaluation [33][34][35] and geophysics [36][37][38].
In this paper, we present the principle and first experimental proofs of concept of our approach in the context of medical ultrasound imaging. However, the concept can be extended to any field of wave physics for which multi-element technology is available.

A. Experimental measurement
The sample under investigation is a tissue-mimicking phantom which generates an ultrasonic speckle characteristic of human tissue, with point-like targets placed at regular intervals, and two bright 5 mm-diameter inclusions at a depth of 40 mm [ Fig. 1(a)]. A 20 mm-thick layer of bovine tissue is placed on top of the phantom and acts as both an aberrating and scattering layer. This experiment mimics the situation of in-vivo liver imaging in which layers of fat and muscle tissues generate strong aberration and scattering at shallow depths. We acquire the acoustic reflection matrix experimentally using a linear ultrasonic transducer array placed in direct contact with the sample [ Fig. 1(a)]. The simplest acquisition sequence is to emit with one element at a time, and for each emission record with all elements the time-dependent field reflected back from the medium. This canonical basis was first used to describe the so-called time-reversal operator [39], and is now commonly used in non destructive testing where it is referred to as the full matrix capture sequence [40]. A matrix acquired in this way can be written mathematically as R uu (t) ≡ R(u out , u in , t), where u is the position of elements along the array, 'in' denotes transmission, and 'out' denotes reception. Alternately, the response matrix can be acquired using beamforming (emitting and/or receiving with all elements in concert with appropriate time delays applied to each element) to form, for example, focused beams as in the conventional B-mode [1] or plane waves for high frame rate imaging [12]. To demonstrate the compatibility of our method with state-of-the-art medical technology, our data was acquired using plane-wave beamforming in emission and recording with individual elements in reception. The experimental procedure is described in detail in Appendix A. A set of plane waves is used to probe the medium of interest. For each plane wave emitted with an incident angle θ in , the time-dependent reflected wavefield is recorded by the transducers. The corresponding signals are stored in a reflection matrix R uθ (t) = [R(u out , θ in , t)] [ Fig. 1(a)]. An ultrasound image can be formed by coherently summing the recorded echoes coming from each focal point r, which then acts as a virtual detector inside the medium. In practice, this is done by applying appropriate time delays to the recorded signals [12]. The images obtained for each incident plane wave are then summed up coherently and result in a final compounded image with upgraded contrast. This last operation generates a posteriori a synthetic focusing (i.e. a virtual source) on each focal point. The compounded image is thus equivalent to a confocal image that would be obtained by focusing waves on the same point in both the transmit and receive modes.

B. Monochromatic focused reflection matrix
We now show how all of the aforementioned imaging steps can be rewritten under a matrix formalism. The reflection matrix can actually be defined in general as containing responses between one or two mathematical bases. The bases implicated in this work are: (i ) the recording basis which here corresponds to the transducer array, (ii ) the illumination basis which is composed of the incident plane waves, and (iii ) the focused basis in which the ultrasound image is built. In the frequency domain, simple matrix products allow ultrasonic data to be easily projected from the illumination and recording bases to the focused basis where local information on the medium proprieties can be extracted.
Consequently, a temporal Fourier transform should be first applied to the experimentally acquired reflection matrix to obtain R uθ (ω), where ω = 2πf is the angular frequency of the waves. The matrix R uθ (ω) can be expressed as follows: where the matrix Γ, defined in the focused basis, describes the scattering process inside the medium. In the single scattering regime, Γ is diagonal and its elements correspond to the medium reflectivity γ(r). T = [T (r, θ)] is the transmission matrix between the plane wave and focused bases. Each column of this matrix describes the incident wavefield induced inside the sample by a plane wave of angle θ. G = [G(u, r)] is the Green's matrix between the transducer and focused bases. Each line of this matrix corresponds to the wavefront that would be recorded by the array of transducers along vector u if a point source was introduced at a point r = (x, z) inside the sample. The holy grail for imaging is to have access to these transmission and Green's matrices. Their inversion, pseudo-inversion, or more simply their phase conjugation can enable the reconstruction of a reliable image of the scattering medium, thereby overcoming the aberration and multiple scattering effects induced by the medium itself. However, direct measurement of the transmission and Green's matrices T and G would require the introduction of sensors inside the medium, and therefore these matrices are not accessible in most imaging configurations. Instead, sound propagation from the plane-wave or transducer bases to the focal points is usually modelled assuming a homogeneous speed of sound c. In this case, the elements of the corresponding free-space transmission matrix T 0 (ω) are given by where x and z describe the coordinates of r in the lateral and axial directions, respectively [ Fig. 1(a)], and k = ω/c is the wave number. The elements of the freespace Green's matrix G 0 (ω) are the 2D Green's functions between the transducers and the focal points [41] G 0 (r out , u out , ω) where H (1) 0 is the Hankel function of the first kind. T 0 (ω) and G 0 (ω) can be used to project the reflection matrix R uθ (ω) into the focused basis. Based on Kirchhoff's diffraction theory [42], such a double focusing operation can be written as the following matrix product where the symbol * and † stand for phase conjugate and transpose conjugate, respectively. The matrices T † 0 and G * 0 contain the phase-conjugated wavefronts that should be applied at emission and reception in order to project the reflection matrix into the focused basis both at input (r in ) and output (r out ). Equation (4) thus mimics focused beamforming in post-processing in both emission and reception. Each coefficient of R rr = [R(r out , r in )] is the impulse response between a virtual source at point r in and a virtual detector at r out [ Fig. 1(c)].
The aberration issue in imaging can be investigated by expressing the matrix R rr mathematically using Eqs. (1) and (4): where and are the input and output focusing matrices, respectively [ Fig. 1 and H out = [H out (r, r out )] corresponds to the transmit and receive PSFs, i.e. the spatial amplitude distribution of the input and output focal spots. Their support defines the characteristic size of each virtual source at r in and detector at r out [ Fig. 1(b)]. In the absence of aberration, the transverse and axial dimension of these focal spots, δx 0 and δz 0 , are only limited by diffraction [43]: where β is the maximum angle of wave illumination or collection by the array and λ the wavelength. In the presence of aberration, i.e. if the velocity model is inaccurate, there is a mismatch between the transmission and Green's matrices, T(ω) and G(ω), and their free-space counterparts, T 0 (ω) and G 0 (ω). The focusing matrices, H in and H out , are far from being diagonal [Eqs.(6)- (7)]. The corresponding PSFs are strongly degraded and the virtual transducers can overlap significantly. A better model of wave propagation is thus needed to overcome aberrations and restore diffraction-limited PSFs. In Sec. III A, we will show how a multilayer model can be used to reach a better estimate of T and G in the experimental configuration depicted in Fig. 1(a).

C. Broadband focused reflection matrix
For broadband imaging, we can restrict our study to pairs of virtual transducers, r in = (x in , z) and r out = (x out , z), located at the same depth z. Furthermore, an inverse Fourier transform of the corresponding submatrices R xx (z, ω) = [R(x out , x in , z, ω)] should be performed in order to recover the excellent axial resolution of ultrasound images. For direct imaging, only echoes at the ballistic time (t = 0 in the focused basis) are of interest. This ballistic time gating can be performed via a coherent sum of R xx (z, ω) over the frequency bandwidth δω. A broadband focused reflection matrix R rr (z) is thus obtained at each depth z: where ω ± = ω c ± δω/2 and ω c is the central frequency.
Each element of R xx (z) contains the signal that would be recorded by a virtual transducer located at r out = (x out , z) just after a virtual source at r in = (x in , z) emits a pulse of length ∆t = δω −1 at the central frequency ω c . The broadband focusing operation of Eq. (9) gives virtual transducers which now have a greatly reduced axial dimension δz ∼ c∆t [ Fig. 1(c)]. Figures 2(a) and (b) display R xx at the bovine tissue/phantom interface (z = 18 mm) and in the phantom (z = 30 mm), respectively. In both cases, most of the signal is concentrated around the diagonal. This indicates that single scattering dominates at these depths [22], since a singly-scattered wavefield can only originate from the point which was illuminated by the incident focal spot. In fact, the elements of R rr which obey r in = r out hold the information which would be obtained via multifocus (or confocal) imaging, in which transmit and receive focusing are performed at the same location for each point in the medium. A line of the ultrasound image can thus be directly deduced from the diagonal elements of R xx (z), computed at each depth: The corresponding image is displayed in Fig. 2(c). It is equivalent to the coherent compounding image computed via delay-and-sum beamforming of the same data set [12], constituting a validation of our matrix approach for imaging. Interestingly, the matrix R rr contains much more information than a single ultrasound image. In particular, focusing quality can be assessed by means of the off-diagonal elements of R xx . To understand why this is, R xx shall be expressed theoretically. To that aim, a time-gated version of Eq. (5) can be derived: where H in (z), Γ(z) and H out (z) are the time-gated sub-matrices of H in [Eq. (6)], Gamma [Eq. (1)] and H out [Eq. (7)] at depth z = ct/2 and central frequency ω c . In the single scattering regime and for spatiallyand frequency-invariant aberration, the previous equation can be rewritten in terms of matrix coefficients as follows: This last equation confirms that the diagonal coefficients of R xx (z), i.e. a line of the ultrasound image, result from a convolution between the sample reflectivity γ and the confocal PSF H in ×H out . As we will see, access to the offdiagonal elements of R xx will allow our analysis of the experimental data to go far beyond a simple image of the reflectivity. In particular, off-diagonal elements can be used to extract the input-output PSF in the vicinity of each focal point, which will lead to a local quantification of the focusing quality.

III. LOCAL FOCUSING CRITERION
In this section, we detail how an investigation of the off-diagonal points in R xx can directly provide a focusing quality criterion for any pixel of the ultrasound image. To that aim, the relevant observable is the mean intensity profile along each antidiagonal of R xx : where · · · denotes an average over the pairs of points r in = (x in , z) and r out = (x out , z) which share the same midpoint r = (r out + r in )/2, and ∆x = (x out − x in ) is the relative position between those two points. We term I(r, ∆x) the common-midpoint intensity profile.
Whereas I (r) [Eq. (10)] only contains the confocal intensity response from an impulse at point r, I(r, ∆x) is a measure of the spatially-dependent intensity response to an impulse at r. This means that whatever the scattering properties of the sample, I(r, ∆x) allows an estimation of the input-output PSFs. However, its theoretical expression differs slightly depending on the characteristic length Homogeneous model Optimized two-layer model scale l γ of the reflectivity γ(r) at the ballistic depth and the typical width δx of the input and output focal spots.
In the specular scattering regime [l γ >> δx, see Fig.  2(a)], the common-midpoint intensity is directly proportional to the convolution between the coherent input and output PSFs, H in and H out (see Appendix B): where the symbol * stands for convolution. However, in ultrasound imaging, scattering is more often due to a random distribution of unresolved scatterers. In this speckle regime [l γ < δx, see Fig. 2(b)], the common midpoint intensity I(r, ∆x) is directly proportional to the convolution between the incoherent input and output PSFs, |H in | 2 and |H out | 2 (see Appendix B): The ensemble average in Eq. 15 implies access to several realizations of disorder for each image, which is often not possible for most applications. In the absence of multiple realizations, a spatial average over a few resolution cells is required to smooth intensity fluctuations due to the random reflectivity of the sample while keeping a satisfactory spatial resolution. To do so, a spatially averaged intensity profile I av (r, ∆x) is computed at each point r of the field of view, such that I av (r, ∆x) = I(r , ∆x) (r −r)∈A (16) where the symbol · · · denotes an average over the set of focusing points r contained in an area A centered on r. The compromise between intensity fluctuations and spatial resolution guided our choice of a 7.5 mm-diameter disk for A. Whatever the scattering regime, the averaged common-midpoint intensity profile I av (r, ∆x) is a direct indicator of the focusing resolution at each point r of the medium. One can then build a local focusing parameter where the input-output focusing resolution w(r) is defined as the full width at half maximum (FWHM) of I av (r, ∆x), and w 0 (r) is a reference value based on the theoretical diffraction limit for a homogeneous medium. This parameter is bounded between 0 (w >> w 0 , bad focusing) and 1 (w = w 0 , perfect focusing). F (r) is the equivalent in the focused basis of the coherence or focusing factor originally introduced by Mallart and Fink [8] in the transducer basis. The definition of a focusing parameter in the focused basis offers an important advantage in that the wave focusing quality and the image resolution can now be probed locally. Figure 3(a) shows the focusing criterion F (r) calculated for the bovine tissue/phantom system [corresponding to the reflectivity image of Fig. 2(c)]. The reference resolution w 0 (r) has been computed under a speckle scattering hypothesis [Eq. (B2)]. The poor quality of focus over a large part of the image can be attributed to the fact that the presence of the bovine tissue layer was not taken into account in our (homogeneous) model of the system [Eqs.
(2) -(4)]. Fortunately, as discussed in the following section, the focused reflection matrix approach enables the determination of a more accurate model for this a priori unknown medium.

A. Wave velocity mapping
The ability to locally probe the focus quality offers enormous advantages for local characterization of heterogeneous media, in particular for a quantitative mapping of their refractive index, or more specifically, as shown in the following, the speed of sound.
We begin by observing the focusing criterion F (r) at the bovine tissue-phantom interface as a function of the wave velocity c t assumed in the bovine tissue. The result is displayed in Fig. 2(e). The corresponding focusing criterion is optimized for c t = 1573 m/s. Fig. 2(g) shows the reflection matrix R xx obtained by considering this optimized wave velocity. The comparison with the orig- inal matrix displayed in Fig. 2(a) shows a narrowing of the input-output PSFs along the antidiagonal of R xx ; the focusing resolution w(r) is now much closer to the diffraction limit w 0 due to the use of the optimized c t . This approach works for a reasonably homogeneous medium (the tissue layer). However, our assumption of a homogeneous wave velocity model does not conform to the bi-layer system under experimental investigation [ Fig. 1(a)]. To probe more deeply into the system, we extend our approach to model a multilayer medium. Using the ultrasound image [ Fig. 2(c)] as an approximate guide, we define two layers: one at z = 0−18 mm with our measured c t , and a second for depths below z = 18 mm with unknown wave velocity c p . New transmission and Green's matrices, T 1 and G 1 , are computed using this two-layer wave velocity model. A new reflection matrix R xx is then built via Figure 2(h) displays R xx at depth z = 30 mm. The corresponding focusing criterion F (r), averaged over the full width of the image and a 2 mm range of depths, is shown as a function of the phantom speed of sound hypothesis c p in Fig. 2(e). The optimization of F yields a quantitative measurement of the speed of sound in the phantom: c p = 1546 m/s. To build an entire profile of wave velocity throughout the medium, the F −optimization is repeated for each depth. The resulting depth-dependent velocity estimate is shown in Fig. 2(d). The presence of two layers can be clearly seen, corresponding to the bovine tissue with mean wave velocity c t = 1570 m/s (z < 18 mm), and the phantom with c p = 1547 m/s (z > 30 mm). These values are in excellent agreement with the manufacturer's specification for the phantom (c p = 1542 ± 10 m/s) and the speed of sound estimated from the travel time of the pulse reflected off of the tissue/phantom interface (c t ≈ 1573 m/s). Near the interface between the two layers, the measurement of c p is less precise because the wave has not spent enough time in the phantom to have had a significant influence on F . More quantitatively, the measurement error ∆c p /c p on the wave velocity scales as the inverse of z p , the depth of the focal plane from the phantom surface [see Appendix C, Eq. C9]: with k p = ω c /c p . A measurement precision ∆F/F ∼ 5 × 10 −4 on the focusing criterion [see Fig. 2(e)] means that the wave has to travel through a medium thickness z p ∼ 10 mm to reach a precision ∆c p ∼ 5 m/s on the wave velocity. This value is in qualitative agreement with the axial resolution of the wave velocity profile displayed in Fig. 2(d). Figure 2(f) shows the ultrasound image deduced from the confocal elements of R xx . Compared with the homogeneous model [ Fig. 2(c)], it can be seen by eye that the two-layer model slightly improves the imaging of bright targets, although there is no clear difference in areas of speckle. However, the result for F (r) after optimization with the two-layer model shows that a significant improvement in quality of focus is obtainable with the twolayer model (Fig. 3). Thus, F (r) constitutes a promising new metric for speed of sound measurement, offering a sensitivity in speckle which is lacking in state of the art methods based on image brightness [44][45][46].

B. Discussion
The study of the focused reflection matrix yields a quantitative and local focusing criterion that constitutes a sensitive guide star to map the wave velocity of an inhomogeneous scattering medium. The perspective of this work will be to go beyond a depth profile of the wave velocity to map its variations in 2D (1D probe) or 3D (2D array). In this respect, we shall mention the recent work of Jaeger et al. [47,48] that investigated the local phase change at a point when changing the transmit beam steering angle. Looking at this local phase change under a matrix formalism would be a way to make the best of the two approaches. An inverse problem would then have be solved to retrieve a map of the local phase velocity [47,48]. Again, a matrix formalism could be relevant to optimize this inversion.
In the same context, we would also like to mention the work of Imbault et al. [24] that investigated the correlations of the reflected wavefield in the transducer basis for a set of focused illuminations. Combined with a time reversal process consisting in iteratively synthesizing a virtual reflector in speckle [49], the wave speed was measured by maximizing a focusing criterion based on the spatial correlations of the reflected wavefield [8]. The downside of this approach is that the construction of the virtual reflector requires the focusing algorithm to be iterated several times before the guide-star becomes point-like. Moreover, the whole process should be both averaged and repeated over each isoplanatic patch of the image [50], which limits the spatial resolution and the practicability of such a measurement. Nevertheless, inspired by previous works [8,16,24,49], novel potential applications can be imagined for the local focusing parameter F (r). It could, for instance, be used as a guide star for a matrix correction of aberration. Based on its maximization, the goal would be to converge towards the best estimators of the transmission matrices T and G. An inversion or pseudo-inversion of these matrices would then lead to an optimized image whose resolution would be only limited by diffraction.
The developments presented thus far have been based on a single scattering assumption. However, multiple scattering is far from being negligible in real-life ultrasound imaging, whether it be for example in soft human tissues [20] or coarse-grain materials [34]. In the following, we show how the reflection matrix approach suggests a solution for the multiple scattering problem, and how it can furthermore be exploited to create a new contrast for ultrasound imaging.

IV. MULTIPLE SCATTERING
In the previous sections, we developed a local focusing criterion by considering the near-confocal elements of R xx . We now turn our attention to the points which are farther from the confocal elements, i.e. for which |x out − x in | > w. Figure 4(a) shows R xx at a depth of z = 60 mm. Signal can clearly be seen at points far from the diagonal. Because each matrix R xx is investigated at the ballistic time (t = 2z/c), the only possible physical origin of echoes between distant virtual transducers is the existence of multiple scattering paths occurring at depths shallower than the focal depth, as sketched in Fig. 4(b). In this section, we will see that a significant amount of multiple scattering takes place in our bovine tissue/phantom system. Signal from such multiple scattering processes has traditionally been seen as a nightmare for classical wave imaging, as it presents as an incoherent background which can greatly degrade image contrast. However, because they are extremely sensitive to the micro-architecture of the medium, multiply scattered waves can be a valuable tool for the characterization of scattering media [20,25,26]. In the following, we show how our matrix approach enables the measurement of a local multiple scattering rate for each pixel of the ultrasound image. This multiple scattering rate can be directly related to the concentration of scatterers, paving the way towards a novel contrast for ultrasound imaging.
A. Multiple scattering in the focused basis For our experimental configuration, multiple scattering can be investigated by examining the averaged spatial intensity profiles I av (r, ∆x) [Eq. (16)]. Each intensity profile is composed of three contributions: 1. The single scattering component, I S . Signals from single scattering mainly lie along the near-confocal elements of R xx [∆r < w(r)]. This is the contribution that has been investigated in the previous sections.
2. The multiple scattering component, I M . This contribution can be split into two terms: An incoherent part which corresponds to interferences between waves taking different paths through the medium, and a coherent part which corresponds to the interference of waves with their reciprocal counterparts [see the blue and red paths in Fig. 4(b). Referred to as coherent backscattering (CBS), this interference phenomenon results in an enhancement (of around two) in intensity at exact backscattering. Originally discovered in the plane wave basis [51][52][53][54], this phenomenon also occurs in a point-to-point basis, whether the points be real sensors [55][56][57] or created via focused beamforming [25,58]. In the point-to-point basis, contributions from multiple scattering give to the backscattered intensity profile the following shape: a narrow, steep peak (the CBS peak) in the vicinity of the source location [∆x < w(r)], which sits on top of a wider pedestal (the incoherent contribution).
3. Electronic noise, I N . These contribution can decrease the contrast of an ultrasound image in the same way as I M . Noise contributes to a roughly constant background level to the backscattered intensity profiles I av (r, ∆x).
To estimate the level of each contribution, the relevant observables are the mean confocal intensity (I on ) and offdiagonal intensity (I off ) of R xx . The confocal intensity I on is given by I on (r) = I av (r, ∆x = 0) = I S (r) + 2I M (r) + I N , (20) where the factor of 2 accounts for the CBS enhancement of the multiple scattering intensity at the source location. I off is the sum of the multiple scattering incoherent background and of the additive noise component: where · · · ∆x>w(r) indicates an average over off-diagonal elements which obey ∆x > w(r). This average constitutes an average over several realizations of disorder, which is necessary to suppress the fluctuations from constructive and destructive interference between the various possible multiple scattering paths through the sample. Figure 5(a) shows three examples of normalized intensity profiles I av (r, ∆x)/I av (r, ∆x = 0). Each profile has been averaged over a different zone of the ultrasound image [ Fig. 5(c)]: green and blue curves (solid and dotted rectangles) correspond to zones situated respectively above and below the bright speckle disk. It is clear that the incoherent background I off is higher in the deeper (blue) zone, suggesting that multiple scattering is greatly enhanced behind the reflective object. Surprisingly, the incoherent background I off is far from being negligible in the red zone at shallower depths (dashed line rectangle).
To investigate these phenomena further, we define two new observables: (1) the multiple-to-single scattering ratio and (2) the multiple scattering-to-noise ratio, To calculate these quantities, it is necessary to be able to distinguish between I M , I N , and I S .

B. Coherent backscattering as a direct probe of spatial reciprocity
Discrimination between I M and I N can be achieved by exploiting the spatial reciprocity of propagating waves in a linear medium. While the multiple scattering contribution gives rise to a random but symmetric reflection matrix R xx [22,23], additive noise is fully random. Thus, the symmetry of R xx gives us a tool to determine the relative weight between noise and multiple scattering in the incoherent background of R xx .
An elegant approach to probe spatial reciprocity is the measurement of the CBS effect in the plane-wave basis (the far-field). The CBS effect can be observed by measuring the average backscattered intensity as a function of the angle ∆θ ≡ |θ in − θ out | between the incident and reflected waves. In the presence of multiple scattering, this profile displays a flat plateau (the incoherent background), on top of which sits a CBS cone centered around the exact backscattering angle ∆θ = 0. The cone is solely due to constructive interference from waves following reciprocal paths inside the sample [ Fig. 5(b), inset]. Thus, CBS in the far-field is a direct probe of spatial reciprocity in the focused basis [25,58].
Appendix D describes how to eliminate the single scattering contribution and extract a far-field intensity profile I(r, ∆θ) for the area A surrounding each focusing point r. In Fig. 5(b), normalized intensity profiles I av (r, ∆θ)/I av (r, ∆θ = 0) are shown for the three areas A highlighted in Fig. 5(c). For each area, a CBS cone is clearly visible, showing that the experimental data contain contributions from multiple scattering. Just as with the CBS peak in the focused basis [ Fig. 5(a)], the highest amount of multiple scattering is observed for the red zone at shallow depths.
To estimate the relative weight of the noise and multiple scattering contributions, we examine the mean intensity for two cases: (1) at exact backscattering I av (r, ∆θ = 0) = 2I M (r) + I N (r), (24) and (2) where θ c is the width of the CBS peak and · · · ∆θ>θc indicates an average over all angles ∆θ which obey ∆θ > θ c . The enhancement factor of the CBS peak is given by χ(r ) can have values ranging from 1 to 2; it is at a minimum when I M = 0 and a maximum when all backscattered signal originates from multiple scattering.

C. Maps of multiple scattering rates
The multiple scattering-to-noise ratio (r) [Eq. 23] can be expressed as a function of the enhancement factor χ(r) by injecting Eqs. 24 and 25 into Eq. 26 : The multiple-to-single scattering ratio ρ(r) [Eq. 22] can be derived by injecting the last equation into Eqs. 20 and 21 Figures 5 (d) and (e) show experimental results for ρ(r) and (r), respectively, superimposed onto the original ultrasound image. Both of these maps constitute new contrasts which are complementary to the reflectivity maps produced by conventional ultrasound imaging. To begin with, ρ(r) can be used as an indicator of the reliability of a reflectivity image [such as that in Fig. 5(c)]. Because the single-to-multiple scattering ratio is a direct indicator of the validity of the single-scattering (Born) approximation, the reliability of the ultrasound image should scale as the inverse of ρ(r ). An interesting example is displayed in Fig. 5(d) at a depth of 37 mm (white solid rounded rectangle), where a high multiple-to-single scattering rate is observed. This abrupt increase of ρ(r) can be accounted for by double reflection events between the probe and the tissue-phantom interface. We can thus conclude that the structures that seem to emerge in Fig.  5(c) at the same depth are in fact artifacts due to multiple reflections.
With respect to the quantification of multiple scattering, the parameter (r) seems to be particularly relevant. The areas in Fig. 5 bubbles that generate resonant scattering, thereby inducing a strong multiple scattering 'tail' behind them; (ii) a set of bright targets close to each other give rise to strong multiply-scattered echoes at depth z =50 mm. The bright speckle area in the phantom (z ∼ 60 mm, white dashed-dotted circle) also contains a sufficient concentration of scatterers to generate multiple scattering events. Figure 5(e) thus demonstrates how the parameter (r) can provide a highly contrasted map of the multiple scattering rate.
An important technical note is that the study of spatial reciprocity in the reflection matrix requires in principle that the bases of reception and emission be identical. Because this is not the case for our experimental measurements, this tends to slightly underestimate I M . Relatedly, here we have employed the CBS effect to probe spatial reciprocity, but an equivalent measurement can be performed directly in the focused basis by computing correlations between symmetric elements of R xx .

D. Discussion
The maps of ρ(r) and (r) help to provide an overall assessment of the factors impacting image quality. For instance, they can be used to explain the apparent poor focus quality in some areas of Fig. 3(b), which appear even after correction for wavefront distortion. The compensation for an incorrect hypothesis for c does not compensate for the effect of multiple scattering or reflections [e.g. in the highlighted areas of Fig. 5(e)]. Thus, the multiple scattering rates combined with the measurement of F (r) provide a sensitive local mapping of the heterogeneities in the medium which includes both small-and large-scale variations of the refractive index.
Beyond image reliability, the maps displayed in Fig.  5(d,e) provide a great deal of quantitative information about the system under investigation. Because (r) is calculated from the off-diagonal elements of R xx , it contains only negligible contributions from single scattering, and thus constitutes an interesting new contrast for imaging which is much more sensitive to the microstructure of the medium than it is to its reflectivity. Conversely, because ρ(r) is independent of noise, it can constitute a useful biomarker for medical imaging deep inside tissue, and a potentially valuable tool for future research seeking to characterize multiple scattering media, even at large depths where I N becomes important. The perspective of this work will be to extract from ρ(r) quantitative maps of scattering parameters such as the elastic mean free path or the absorption length [20,34], and transport parameters such as the transport mean free path [55,59,60] or the diffusion coefficient [55,56,61]. While diffuse tomography in transmission only provides a macroscopic measurement of such parameters, preliminary studies have demonstrated how a reflection matrix recorded at the surface can provide transverse measurements of transport parameters [25,26,58,62]. The focused reflection matrix we have introduced here connects each point inside the medium to all other points. Hence, a 2D or 3D map of transport parameters can now be built by solving the radiative transfer inverse problem.

V. CONCLUSION AND PERSPECTIVES
In summary, a powerful and elegant matrix approach of wave imaging has been presented. By focusing at distinct points in emission and reception, one can build a focused reflection matrix that contains the impulse responses between a set of virtual transducers mapping the entire medium. From this focused reflection matrix, a local focusing parameter can be estimated at any point of the inspected medium. Because it can be applied to any type of media, including in the speckle regime or in the presence of specular reflectors, this focusing criterion is suitable to any situation encountered in medical ultrasound, and enables wave velocity mapping of the medium. Knowledge of the spatial variation of velocity can then in turn be used with the focused reflection matrix to overcome wavefront distortions. The contrast and resolution of the image could then be restored almost as if the inhomogeneities had disappeared. This physically intuitive approach does not depend on arbitrary parameters such as image quality, and does not require guide stars or complex iterative adaptive focusing schemes.
Furthermore, we have shown that the focused reflection matrix enables a local examination of multiple scattering processes deep inside the medium. We have demonstrated the effectiveness of using fundamental interference phenomena such as coherent backscattering -a hallmark of multiple scattering processes -to discriminate between multiple scattering and measurement noise. A novel imaging method is proposed based on the multiple scattering contrast. Unexplored but promising perspectives for this work include the quantitative imaging of parameters such as the scattering, absorption or transport mean free paths. Finally, we emphasize that we have only concentrated here on the relationship between virtual transducers located in the same focal planes at the ballistic time; there is enormous further potential for the analysis of the entire focused matrix R rr across the whole medium and beyond the ballistic time, which will be explored in future works. Last but not least, our matrix approach of wave imaging is very general and could be applied to any kind of system for which emission and detection of waves can be varied in a controllable way.
Thus, the potential of this work goes far beyond ultrasound imaging, with immediate foreseeable impacts in a range of wave physics including optical microscopy, radar and seismology. The experimental set up consisted in an ultrasound phased-array probe (SuperLinear TM SL15-4) connected to an ultrafast scanner (Aixplorer R , SuperSonic Imagine, Aix-en-provence, France). This 1D array of 256 transducers with a pitch p = 0.2 mm was used to emit 41 plane waves with an angle of incidence θ in spanning from −20 o to 20 o [ Fig. 1(a)]. The emitted signal was a sinusoidal burst of central frequency f c = 7.5 MHz, with a frequency bandwidth spanning from 2.5 to 10 MHz. In reception, all elements were used to record the reflected wavefield over a time length t = 124 µs at a sampling frequency of 30 MHz. The matrix acquired in this way is denoted R uθ (t) ≡ R(u out , θ in , t).
Appendix B: Derivation of the common-midpoint intensity profile The theoretical expression of the common-midpoint intensity profile I(r, ∆x) is derived in the specular and speckle scattering regimes.
In the specular scattering regime, the characteristic size l γ of reflectors is much larger than the width of the focal spot δx. γ(r) can thus be assumed as invariant over the input and output focal spots. Equation 12 then becomes R(r, ∆x) = γ(r) × (H in * H out ) (∆r). (B1) The injection of Eq. B1 into Eq. 13 yields the expression of I(r, ∆x) given in Eq. 14.
In the speckle scattering regime, l γ << δx. This is the most common regime in ultrasound imaging, as scattering is more often due to a random distribution of unresolved scatterers. To a first approximation, such a random medium has the property that γ(r 1 )γ * (r 2 ) = |γ| 2 δ(r 2 − r 1 ), where δ is the Dirac distribution. The combination of Eqs. 12, 13 and B2 gives directly I(r, ∆x) as expressed in Eq. 15.
Appendix C: Measurement errors on the focusing criterion and the speed of sound In this work, we have defined the focusing parameter F as the ratio between the width w 0 of the ideal diffraction-limited PSF and the width w of the experimentally measured PSF. In optics, the Strehl S ratio is generally used to quantify aberration [63]. It is defined as the ratio between the maximum of the PSF intensity, I, and that in the ideal diffraction-limited case, I 0 . Due to energy conservation, we have I 0 × w 0 = I × w. The focusing criterion and Strehl ratio, as well as their relative measurement errors, are thus equivalent: S can also be expressed as the square magnitude of the averaged aberration transmittance e iφ(sin θ) [63]: where φ(sin θ) is the far-field phase delay induced by the mismatch between the propagation model and the real medium in the θ-direction. In Fig. 2, a two-layer medium is used to model the bovine tissue/phantom system. Assuming that the wave velocity c t is properly estimated in the first layer (bovine tissue), the phase φ(sin θ) accumulated in the phantom is given by φ(sin θ) = k p z p cos(θ p ), where k p = ω/c p is the wavenumber in the phantom and θ p is the refraction angle in the phantom, obeying sin θ p /c p = sin θ/c t . If a wrong value of c p is used to model sound propagation in the phantom, the resulting phase distortion is given by where η = ∆c p /c p is the relative error of the speed of sound hypothesis in the phantom. For the sake of simplicity, we will assume in the following that cos θ ∼ cos θ p . This approximation is justified by the small relative difference between c p and c t . Assuming relatively weak aberrations (∆φ(sin θ) << π), the transmittance aberration function e i∆φ(sin θ) can be expanded as e i∆φ(sin θ) ∼ 1 − i k p z p cos θ η − 1 2 The angular average of e iφ(sin θ) is then deduced e i∆φ(sin θ) sin θ = 1 sin β sin β 0 e i∆φ(sin θ) d(sin θ) Injecting the last expression into Eq. C3 leads to the following expression of the Strehl ratio : For weak aberrations (F, S ∼ 1), the relative error ∆F/F [Eq. C2] of the focusing criterion can then be directly deduced from the previous expansion of the Strehl ratio: