Method for local compensation of aberrations in a dynamic medium in ultrasound imaging.

Abstract

The invention relates to an ultrasonic method for constructing a confocal image of a dynamic medium, the method comprising the following steps: a) acquisition of a series of canonical reflection matrices Rui(t, μ), b) focused reflection matrices Rxx(z, μ), c) dynamic focused reflection matrices, d) correction laws, e) corrected focused reflection matrices, f) dynamic confocal signal Sc(x,z) from any point at spatial position (x, z), g) construction of an image from the dynamic confocal signals. Figure for the abstract: Fig. 3

Description

Title of the invention: Method for local compensation of aberrations in a dynamic medium in ultrasound imaging. technical field

[0001] The present invention relates to a method and a system for constructing a confocal image of a dynamic medium with local compensation of aberrations.

[0002] The invention is advantageously applicable to the field of medical imaging, but can be applied to any field of ultrasound imaging. State of the art

[0003] In the field of acoustic imaging, the aim is to characterize an unknown medium by actively probing it with ultrasonic waves. This is notably the principle of the ultrasound scanner in medical imaging.

[0004] However, due to inhomogeneities in the speed of sound between the different tissues of the human body, ultrasound images suffer from both transverse and axial aberrations which degrade their contrast and alter their resolution.

[0005] Furthermore, in ultrasonic localization microscopy, it is possible to image the vascular network of an organ with high resolution by detecting, locating, and tracking isolated bubbles. However, this bubble detection and tracking process relies on convolution with a Gaussian spreading function, which is only valid in an ideal case without aberrations.

[0006] Figure 1 illustrates a conventional focusing process for producing an ultrasound image of a medium. Unfortunately, the medium here contains an aberrant layer with a different speed of sound than the speed of sound observed in the rest of the medium. This results in spatial distortion and temporal dispersion (reverberation) of the acoustic wavefront, leading to transverse and axial aberrations in the resulting ultrasound image. These phenomena lead to a degradation of its resolution and contrast, as well as the appearance of reverberation artifacts, which are particularly problematic during a medical examination.

[0007] In the diagram on the left of this [Fig. 1], the transducer array, positioned opposite a medium, allows for the insonification and imaging of the medium. The conventional method consists of insonifying the medium using focused emissions by a technique known as beamforming. A set of appropriate delays r, based on a homogeneous velocity model c0, is applied to the signals emitted by each transducer in order to constructively interfere the waves produced by each transducer at the targeted focal point with spatial position rin. = (x in, z). Due to the Due to the physical limitations of diffraction, the ultrasound waves emitted through the aperture of the ultrasound probe are concentrated in an area often called the "focal spot." Furthermore, the waves passing through the aberrant layer are distorted, causing distortion and broadening of the focal spot around the focal point. This undesirable effect is illustrated in the diagram on the left of [Fig. 1].

[0008] Waves reflected at the focal point are returned to the transducer array and pass through the aberrant layer again, further distorting the reflected wavefront measured by the transducer array. A path-forming process applied to this type of signal results in an ultrasound image exhibiting significant lateral distortions due to the presence of the aberrant layer. If this layer is reverberant, multiple reflection echoes can lead to axial distortion of the ultrasound image. These various effects result in a loss of resolution and contrast in the ultrasound image. A heterogeneous distribution of sound velocity in the traversed tissues therefore impacts the quality of the reconstructed image.

[0009] The diagram on the right of [Fig. 1] illustrates the object of the invention: determining the wavefront to be emitted in order to optimally focus the ultrasonic waves both spatially and temporally towards each point in the medium. Adaptive focusing techniques, or more recently matrix imaging, have been developed for this purpose. However, they rely on the invariance of the focal spot over a sufficient area to intelligently combine the waves reflected by different contiguous points. This effectively eliminates disorder and allows access to the aberration law associated with the area under consideration, commonly called the isoplanetary zone. These adaptive focusing techniques are well known but remain limited because they only allow compensation for relatively low-order aberrations associated with sufficiently large isoplanetary patches.Higher-order aberrations and reverberations vary too rapidly to be detected by these state-of-the-art techniques. This results in a spatiotemporal distortion of the acoustic wavefront, leading to significant aberrations in the ultrasound image, and therefore a degradation of its resolution and contrast. These aberrations can be so severe that they compromise the ultrasound characterization, particularly in the case of a medical examination.

[0010] The present invention relates to a new ultrasound imaging method in which higher-order aberrations and reverberations are optimally compensated for each focal point. The objective is to obtain an ultrasound image with the highest possible resolution and optimal contrast. Description of the invention

[0011]

[0012]

[0013]

[0014]

[0015]

[0016]

[0017]

[0018]

[0019]

[0020] At least one of the aforementioned objectives is achieved with an ultrasonic method for constructing a confocal image of a dynamic medium, the method comprising the following steps: a) acquisition, by means of a network of transducers, of a series of canonical reflection matrices R ui(t, #m)=[Æ(u out,i mt,# ,„)] at different times, each canonical reflection matrix is ​​defined between an ultrasonic wave emission basis i at the input and a reception basis u at the output; the coefficients of this canonical reflection matrix correspond to the signals received by the transducers and induced by the ultrasonic waves reflected in the medium; t denoting the echo time and #m denoting the mth canonical reflection matrix; b) determination of a focused reflection matrix R xx(z, #m) for each canonical reflection matrix by input-output focusing for any point of at least one region of the medium, the coefficients of this focused reflection matrix are obtained by calculating an acoustic pressure field between all points of the region with lateral positions x in and x out, located at an expected depth z for a speed of sound assumed to be c0; c) determination of a dynamic component for each coefficient of each focused reflection matrix so as to constitute dynamic focused reflection matrices # ) vV.4- lit / d) determination of a correction law Q(x. z) for each point x and depth z of the medium from the dynamic focused reflection matrices, e) Determination of corrected focused reflection matrices #J By applying the correction law 0(x, z) at every point in the medium, f) determination of a dynamic confocal signal S c (x^ # m) of any point of spatial position (x, z) from the diagonal coefficients of the corrected focused reflection matrix R^. / z, # m), g) construction of an image from dynamic confocal signals. With the present invention, the method advantageously allows for local probing of the medium to obtain a local estimate of an aberration correction law adapted to correct the ultrasonic channel formation process. This correction makes it possible to reduce or eliminate aberrations, for example, due to variations in the speed of sound in the medium or to multiple reflections of waves generated by one or more aberration zones in the medium. The acquisition step can consist of performing measurements to process the reflection matrices or of a posteriori processing the matrices from data stored in a memory space. The variable t represents the echo time associated with the signals recorded during the measurement and f the ultrasonic wave frequency during the measurement. For each measurement, an amplitude and a phase are acquired for each pixel.

[0021] Thus the various calculations of the invention can be carried out independently of the measurement acquisition phase, in particular by modifying various calculation parameters, which makes it possible to carry out various ultrasonic construction analyses either in real time or a posteriori.

[0022] These correction calculations benefit from local information extracted from the dynamic part of the echoes reflected by the tissues.

[0023] While existing aberration correction methods (adaptive focusing, matrix imaging) rely on an assumption of local isoplanetism (spatial invariance of the focal spot) in order to average the correlations between ultrasonic signals from several contiguous focal points, the method according to the invention uses the dynamics of the signals to obtain an independent focusing law for each pixel. Spatial averaging is not required, and it is therefore possible to access very high-order aberrations exhibiting little or no isoplanetism.

[0024] The invention is particularly remarkable in that a series, and therefore several canonical reflection matrices, are acquired at different times. Each acquisition can comprise several measurements with the application of several different plane waves to constitute a matrix. An acquisition can be called a frame, symbolized by #m.

[0025] Having several canonical reflection matrices of the same medium acquired at different times makes it possible to take into account the dynamic aspect of the medium. This dynamic aspect, which could be considered a disadvantage, is used to improve the accuracy in determining the correction laws.

[0026]

[0027] According to an advantageous embodiment of the invention, the step of acquiring a series of canonical reflection matrices R ui(t, #m) can include the emission of an ultrasonic pulse from each transducer of the network whose position is located by the coordinate u in; this pulse gives rise to a divergent cylindrical or spherical incident wave which is reflected by diffusers of the medium; these reflected echoes form a backscattered field which is recorded by each of the transducers as a function of time; the canonical reflection matrix RUu(t># m) expressed in the basis of the transducers being composed of a set of impulse responses A*(u out,u imt,# m ) between transducers.

[0029] According to one embodiment, the step of acquiring a series of canonical reflection matrices Ru(t, #m) may include insonifying the medium with a series of plane waves with a delay r' applied to each signal at emission to form a wavefront inclined at an angle θin with respect to the transducer array. A backscattered field by the medium, θ(uout, θia, t, #m), is measured by all the position transducers uout for each incident plane wave θin. The set of responses forms a canonical reflection matrix Ru0(t, #, θ) = [θ(uout, in - t, #m)] -

[0030]

[0031] According to yet another variant, the step of acquiring a series of canonical reflection matrices R ui(t, #m) may include an insonification of the medium with a series of divergent waves.

[0032]

[0033] According to the invention, the step of determining the focused reflection matrix R xx(z, # m) may include:

[0034] - an input focusing process from each reflection matrix canonical R ui(t,# m) which uses a forward time of flight of the waves between the ultrasonic wave emission base i and a virtual input transducer TVin and which creates a so-called input focal spot around a first point PI of spatial position r in =(x in,z), said input focal spot corresponding to the virtual input transducer TVin, and

[0035] - an output focusing process from the canonical reflection matrix R m) which uses a return time of the waves between a virtual output transducer TVout and the receiving base transducers u and which creates a focal spot called output around a second point P2 of spatial position r out =(x out,z), said focal spot output corresponding to the virtual output transducer TVout.

[0036]

[0037] By way of example, the xout coordinate of the second point P2 is located at a distance from the xin coordinate of the first point PI which is less than or equal to a maximum distance Axmax which is a function of a number of ultrasonic waves generated during the acquisition of the reflection matrix.

[0038]

[0039] In particular, the focused reflection matrix can be determined in the time domain or in the frequency domain.

[0040] If we are dealing with aberrations that resemble simple time shifts of relatively small amplitude (typically less than the temporal resolution of ultrasonic signals), then it will be advantageous to directly calculate the reflection matrix in the time domain and only at ballistic time.

[0041]

[0042]

[0043]

[0044]

[0045]

[0046]

[0047]

[0048] If we are dealing with larger amplitude aberrations and / or reverberations and / or frequency dispersion, then it is preferable to adopt a polychromatic approach and calculate the focused reflection matrix in the frequency domain. According to an advantageous feature of the invention, the dynamic component can be determined by subtracting from each coefficient of each focused reflection matrix a moving average over N focused reflection matrices. Furthermore, the dynamic component can be determined by applying a high-pass or band-pass filter along a dimension of the numbers #m of the reflection matrices. We are taking advantage here of having acquired several frames of canonical reflection matrices. The dynamic component can also be determined by performing a singular value decomposition of the focused reflection matrices rearranged into a two-dimensional matrix, one dimension of which is that of the numbers #m in the reflection matrices. In other words, the set of reflection matrices is concatenated into a global two-dimensional matrix, which is then analyzed as a singular value decomposition.

[0049]

[0050]

[0051] According to the invention, the step of determining a correction law 0(x, z) can include, for each dynamic focused reflection matrix Rx^#, the following steps: - determination of a dual reflection matrix Rcx(z, # m) by forward projection directly or indirectly from the dynamic focused reflection matrix Rd # j to a correction basis (c), - calculation of the correction law 0(x, z) from the dual reflection matrix Rcx( z, # m), said correction law being a correction law, 0 = [^(e. x)] on the correction basis (c), - determination of the corrected focused reflection matrices r'Æ #J by back projection of the corrected dual reflection matrices R^x(z, # m) towards the focused basis (x), that is to say the set of points x at the depth z considered.

[0052]

[0053]

[0054] Preferably, the coefficients R^x(z_ # m) = C Z- #œ) ] of the corrected dual reflection matrix R^. can be determined by performing a term product to term between the dual reflection matrix Rcx(z, # m) and the phase conjugate of the correction law z), i.e.: R^. = R^o $*

[0055] where the symbol * denotes a phase conjugation operation, the symbol ° is the Hadamard product, such that:

[0056]

[0057] The step of calculating the correction law (fXx, z) may include the following steps:

[0058] - construction of a correlation matrix C(x,z) from reflection matrices dual Rcx(z,, # m) for each point (x,z) of a field of view,

[0059] - determination of the focusing law ^(XjZ) for each point of the field of vision by performing one of the following operations:

[0060] - eigenvalue decomposition of the correlation matrix C(x,z), the law of correction ®(x,z) being the first eigenvector Uj of the correlation matrix C(x,z) in the correction basis (c),

[0061] - singular value decomposition of the rearranged dual reflection matrix of in the following way:

[0062] z) = [r (c # mx, z)]

[0063] the correction law 0(x5 z) being equal to the first singular vector of the dual reflection matrix Rc#, i.e. z) = U1 - solving the following equation: 0(x, z) =exp(j argfC^x, z)xO(x, z)} )

[0064] iteratively by the following expression, which corresponds to an iterative phase-reversal calculation:

[0065] (x, z) - exp( j arg{Ccc x 0W(x, z)})

[0066] Where x is the matrix product, with ¢0 an arbitrary wavefront,

[0067] the correction law z) being obtained by:

[0068] 0(x, z) =limG„(x, z). tr^ao - solving the following equation: W(x, z) = exp(j arg[C##(x, z) x W(x, z)}) where x is the matrix product,

[0069] iteratively by the following expression:

[0070] W„+ ] ( x, z ) = exp ( / arg {C## ( x, z ) x W„ ( x, z )} ) where x is the matrix product, with Let Wo be an arbitrary wavefront, which allows us to obtain the following vector W(x, z):

[0071] W(x, z) = limW„(x, z) n-^ao

[0072] the correction law 0(x. z) being obtained by:

[0073] $(c,x,z) =exp{jxarg{Y,xzR(x,c,^ #a)W\x,c,z #,„)}}•

[0074] Advantageously, the correlation matrix C(x,z) can be determined in the correction basis c and in the frequency domain, by the following calculation of the elements of the correlation matrix C = C cc:

[0075] C(c,c',x,z) = YmR(x,c,z, #m) R"{x,c,z, #m)

[0076] where * is the conjugation operator.

[0077] c and c' being points of the correction basis c

[0078] This operation allows the reflected fields in the correction basis to be correlated for each virtual source in (x,z). This correlation is averaged over the different frames # m, that is to say the different realizations of the speckle, in order to overcome the random reflectivity of the medium and thus synthesize a coherent guide star from the different realizations # m of the speckle.

[0079] Alternatively, the correlation matrix C(x,z) can be determined in the number base # m, by the following calculation of the elements of the correction matrix C = C##;

[0080] C( # # i,X,z) =YtcR(xyC,Z, # m) # / )

[0081] * is the conjugation operator.

[0082] c and c' being points of the correction basis c

[0083] #m and # / denoting the same and same frames of the sequence of recorded reflection matrices.

[0084] Preferably, the step of calculating the correction law G(x, z) can be iterated at least twice with, at each iteration, the forward projection uses the corrected focused reflection matrix R'vv( # m ) obtained during the back projection of the previous iteration instead of the focused reflection matrix Rxx(z, # m).

[0085] The step of determining a dual reflection matrix Rcx(z, # m) can also be iterated at least twice with, at each iteration, the use of a different correction basis c.

[0086]

[0087] The step of determining a dual reflection matrix Rcx(z, # m) can otherwise be iterated at least twice with, at each iteration, the use of a forward projection either towards an input correction basis, or towards an output correction basis of the dynamic focused reflection matrix.

[0088]

[0089] According to an advantageous feature of the invention, the correction basis can be one of the following bases:

[0090] - a plane wave basis or spatial Fourier basis,

[0091] - a base of the u transducers,

[0092] - a base corresponding to the supposed location of aberrators in the medium,

[0093] - a basis corresponding to a plan determined by optimization.

[0094]

[0095]

[0096]

[0097]

[0098]

[0099]

[0100]

[0101]

[0102]

[0103]

[0104]

[0105] This allows projections to be made at the level of transducers, aberrators or other. The step of determining a dual reflection matrix Rcx(z, # m) can be performed by forward projection of the dynamic focused reflection matrix \ towards the correction basis (c) by considering a model propagator describing the propagation of waves from the focused basis (x) to the correction basis (c). According to one feature of the invention, the reflection matrix considered at the start of the process can be the broadband focused reflection matrix, which can be obtained by digital channel formation in the time domain or in the Fourier domain from frequency matrices: RlL« = 0, z,f, # m y This broadband reflection matrix corresponds to a focused reflection matrix windowed around the ballistic time. The step of determining a broadband dual reflection matrix Rcx(z, # m) can be performed by forward projection of the broadband focused reflection matrix obtained from the dynamic focused reflection matrix R^^ t — 0 # m) onto the correction basis (c) by considering a propagator at the single center frequency y — + f} / 2. The resulting aberration law then does not exhibit any frequency dependence, ^(c, f,x,z) = <p(c, x, z). Cette variante est par exemple à considérer si les aberrations sont assimilables à de simples décalages temporels d’amplitude inférieure à la résolution temporelle de la mesure ultrasonore, cette dernière variant comme l’inverse de la bande passante. Au-delà, une approche multi fréquence est à privilégier. If the aberrations are only axial (sound speed varying only with depth), the step of determining a dual reflection matrix R(x,f,m) is not necessary. The correction law we are looking for is indeed only frequency-dependent: ∂(x,f,x,z) - ∂(f,x,z). In this case, the matrix considered at the start of the process is the confocal signal of the dynamic focused reflection matrix R^j-#, this confocal signal being such that: S(x,z,f,m) = R(x,x,z,f,m) In this case, the correlation matrix for obtaining the frequency correction law, O(x, z) = [¢(f, x, z)], can be given by: C(f,f',x,z) =lL#S(x,z,f, #m) S~(x,Z,f', #m)

[0106] where the symbol * denotes a phase conjugation operation

[0107] f and f' denote frequencies of the bandwidth of the ultrasonic signal, the correction basis being directly named here (f), i.e. (c ) = ( f),

[0108] S(x, z, f, # m) being the confocal signal of the dynamic focused reflection matrix # J.

[0109] The step of calculating the correction law z), can be carried out:

[0110] - by defining a correlation matrix for obtaining the correction law, <î>( x, z) = [0(f,x,z)], by: [cm] cv. / '.xz) »„) *„)

[0112] - and by averaging the correlation matrix over the numbers #m but also on adjacent pixels belonging to the same isoplanetary patch C(f, f', Xp, Zp) = Z, f, # m) s\x, z, f, # m)

[0113] f and f' denote a frequency of the confocal signal in the bandwidth of the probe

[0114] xp is the transverse coordinate of the central point of the isoplanetary patch for which we seek to estimate the zero aberration law

[0115] zp is the axial coordinate of the central point of the isoplanetary patch for which we seek to estimate the aberration law 0

[0116]

[0117]

[0118] &p: being an isoplanetism patch centered on the point (Xp, zp). According to the invention, the step of determining corrected focused reflection matrices R^^ # can be iterated at least twice with, at each iteration, the forward projection uses the corrected focused reflection matrix R^( # m ) obtained during the back projection of the previous iteration instead of the focused reflection matrix RXXL, # m).

[0119] Preferably, step g) of constructing an image from the dynamic confocal signals may include a step of: - determining the intensity of each dynamic confocal signal to construct a confocal image of the medium, or - determining the intensity of each dynamic confocal signal to construct a power Doppler image by summing, for each point of spatial position (x, z), the intensities of several confocal signals, or - determination of a Fourier transform of the dynamic confocal signal Sc(x, Z, #) according to the time dimension #m: Xc(x, z,v) = IL SAx, z, # where the frequency v is the variable conjugate to the acquisition time #, to construct a directional Doppler image of the medium,

[0120] - determination of the local dynamics for each ultrasound image point of spatial position (x, z), by measuring the average frequency VD (x,z) of the Fourier transform Xc(x, z, v) of the complex confocal signal:

[0121] vD ( X, Z ) = -----—

[0122] to construct a map of the axial velocity of the diffusers at each point of the image.

[0123] According to an advantageous feature of the invention, the focused reflection matrix R xx(z,#m) can depend on a frequency f of the ultrasonic signals, the correction law being determined as a function of this frequency f and the coordinates of the correction plane.

[0124]

[0125] According to an advantageous feature of the invention, the focused reflection matrix R xx(z,#m) can be integrated over the entire bandwidth, the correction law being determined solely as a function of the coordinates of the correction plane (c). The focused reflection matrix, as well as the dynamics, is considered only at the ballistic time.

[0126]

[0127] According to another aspect of the invention, an ultrasonic system for constructing a confocal image of a dynamic medium is provided, the system comprising:

[0128] - an array of transducers adapted to generate a series of ultrasonic waves incident in an area of ​​interest of the medium, and to measure over time the ultrasonic waves backscattered by said area of ​​interest; and

[0129] - a computing unit connected to the transducer network and adapted to put into implements the process described above.

[0130]

[0131] A computer program product is also envisaged, comprising instructions which, when the program is executed by a computer, cause the computer to carry out the steps of the process described above.

[0132]

[0133] A computer-readable medium is also provided, comprising instructions which, when executed by a computer, lead the computer to carry out the steps of the process described above. Brief description of the drawings

[0134] Other advantages and features of the invention will become apparent from the detailed description of implementations and embodiments, which are by no means limiting, and the following attached drawings.

[0135] Fig. 1 is a schematic view illustrating the problem of aberrations in the acquisition of ultrasound images according to the prior art;

[0136] Fig. 2 is a schematic view illustrating an example of an ultrasonic construction system for implementing the method according to the present invention;

[0137] The [Fig.3] is a diagram of the method for constructing an ultrasonic image according to the present invention;

[0138] Fig. 4 illustrates several schematic views 2a to 2f illustrating emission / reception sequences used for ultrasonic imaging and characterization of a medium;

[0139] Fig. 5 is a schematic view showing the focusing principle in the method according to the invention;

[0140] Fig. 6 is a schematic view illustrating the operation of extracting the dynamic component of ultrasonic signals for a transcranial imaging experiment of a sleeping sheep brain;

[0141] Fig. 7 is a schematic view illustrating the operation of extracting aberration laws by iterative phase reversal in the dynamic speckle;

[0142] Fig. 8 represents several images illustrating the effect of applying different local transverse aberration laws in the sheep brain;

[0143] The [Fig.9] include several confocal images before and after correction of a part of the sheep's brain;

[0144] The [Fig. 10] include two "power Doppler" type visualizations showing the gain in contrast and resolution provided by an ultra-local correction of aberrations;

[0145] Figure 11 contains two images of the type of local dynamic frequency maps in the sheep brain, before and after correction of transverse aberrations. Detailed description of the figures

[0146] It is understood that the embodiments described below are in no way limiting. In particular, variants of the invention may be conceived comprising only a selection of the features described below, isolated from the other features described, if this selection of features is sufficient to confer a technical advantage or to differentiate the invention from the prior art. This selection includes at least one preferably functional feature without structural details, or with only a portion of the structural details if this portion alone is sufficient to confer a technical advantage or to differentiate the invention from the prior art.

[0147] The various embodiments and aspects described in this disclosure can be combined or simplified in numerous ways. In particular, the steps of the various processes can be repeated, reversed, and / or carried out in parallel, unless otherwise specified.

[0148] This disclosure relates to methods and systems for the ultrasonic characterization of a medium, and is particularly applicable to medical imaging of living or non-living tissues. The medium may be, for example, a heterogeneous medium that one seeks to characterize in order to, for example, identify and / or characterize heterogeneities. These construction techniques are notoriously non-invasive to the medium, which is advantageously preserved, particularly in its nature and integrity.

[0149]

[0150] Usual approach to ultrasound imaging

[0151] In the field of ultrasound imaging, the aim is often to construct an image of the reflectivity of a medium from echoes backscattered by heterogeneities in the medium. This is the principle of the ultrasound scanner used in medical imaging, which notably allows visualization of the internal anatomy of an individual or an animal. For the sake of simplification, to enable the construction of an ultrasound image, the medium is considered homogeneous, with a constant speed of sound c₀.

[0152] Conventional ultrasound methods generally use an array of piezoelectric transducers that can emit and / or receive ultrasonic signals independently or almost independently, each transducer being at a position u in the array supporting said array. The array of transducers, placed opposite a medium, allows insonification and the construction of a representative image of the medium in various ways. A conventional method consists of insonifying the medium using focused emissions by a technique called beamforming. This method consists of applying to the signals emitted by each transducer a set of appropriate delays x(u_in, x_in, z, co_n) based on a homogeneous velocity model c_0, in order to constructively interfere the wavelets produced by each transducer at the targeted focal point of spatial position (x_in, z).Due to the physical limitations of diffraction, ultrasounds are emitted through the aperture of the ultrasound probe, concentrated in an area often called a "focal spot", with a lateral width ôx.

[0153] In order to potentially construct an ultrasound image illustrating the characteristics of the medium under study, a digital focusing step is also performed at the receiver. The echoes captured by the transducers of the array are realigned by temporally shifting them. The delays r(u out, x out, z, c 0 ) The parameters are identical to those applied during transmission, with the variable uOllt denoting the position of each transducer. During transmission, all signals interfere at the position point (xin, z) at the ballistic time t = z / c0 if the velocity model c0 used corresponds to the reality of the medium being studied. During reception, the signals originating from this same point (xout = xin) interfere by summation at the echo time t = 2z / c0. This summation yields the final focusing result during reception. This confocal method with dual focusing at both transmission and reception allows for direct imaging of the medium's reflectivity with lateral resolution δx and good contrast. However, this method is time-consuming because it requires physical focusing at each point of the medium, or at least at a given depth, on each line of the constructed image representing the medium.

[0154]

[0155] The present invention aims to improve known ultrasonic sounding methods, in particular to correct aberrations.

[0156]

[0157]

[0158] Ultrasonic construction system

[0159] Figure 2 illustrates an example of an ultrasonic imaging system 1 for implementing the ultrasonic imaging method of a medium such as a heterogeneous medium M, according to the present invention. This system and the method allow the formation of an ultrasound image of at least a part (area of ​​interest or field of view) of the medium.

[0160] System 1 comprises:

[0161] - a sounding device 20 or probe 20,

[0162] - a computing unit 30 for calculating an image from the signals received from the probe 20,

[0163] - a control panel 40 connected to the computing unit 30, this control panel including, for example, buttons 41 and a touchpad 42,

[0164] - a display device 50 for viewing an image and various elements or measures.

[0165] The probe 20 is connected to the computing unit 30 via a cable 21 or via a wireless connection, and is capable of emitting ultrasonic waves W into the medium M and receiving ultrasonic waves W from the medium M, said ultrasonic waves being the result of reflections of the emitted ultrasonic waves on diffusing or diffusing particles inside the medium.

[0166] The probe 20 may include an array 10 comprising a plurality of transducers 11. The array 10 may be, for example, a linear, curved, two-dimensional, or matrix array. The transducers 11 are capable of converting an electrical signal into a vibration and vice versa. The transducers 11 may be, for example, piezoelectric ultrasonic transducers in the form of a rigid bar made in direct or indirect contact with an external surface of the medium M to be coupled to the medium and to vibrate and emit and receive ultrasonic waves W. The array 10 of transducers 11 of the probe 20 is then associated with the computing unit 30. The array 10 of transducers may include one hundred or more transducers 11.

[0167] The processing unit 30 may include a housing 31 comprising receiving devices for amplifying and / or filtering the signals received from the probe 20, and converters (analog-to-digital converters and digital-to-analog converters) for transforming the signals into representative signal data. The data may be stored in a memory of the processing unit 30 and / or directly processed to calculate intermediate data (lane formation data or other data). The processing unit 30 may implement any known method for constructing an image from the signal data received from the probe 20, such as lane formation.

[0168] The calculated image can be:

[0169] - a middle image (B-mode image) usually in greyscale for visualize organs in the environment, and / or

[0170] - an image showing a velocity or flow in the medium (color image) by a useful example for visualizing blood vessels in the medium, and / or

[0171] - an image showing a mechanical characteristic of the medium (elasticity) by useful example for identifying tumors within the medium.

[0172] The term "connection" or "link" between the sounding device 20, the computing unit 30, and the display device 50 means any type of wired connection, whether electrical or optical, or any type of wireless connection using any protocol such as WiFi™, Bluetooth™, or others. These connections or links may be one-way or two-way. The associated display device 50 may be of any type, such as a touchscreen or non-touchscreen, connected or not.

[0173] The display device 50 is a screen for viewing the image calculated by the processing unit 30. The display device 50 can also display other information such as image scales, configuration information for calculation or processing, or any measurement or assistance information. The screen 50 can be articulated on a support arm 51 for improved user positioning. A 50-inch screen is usually a large screen (at least 20 inches) for better viewing for the user.

[0174] The control panel 40 is, for example, a portion of a system enclosure, said portion comprising a panel enclosure having a substantially flat surface 40a inclined towards the user for one-handed operation. As shown in [Fig. 2], the control panel 40 may include a control screen 49 for displaying various configuration information.

[0175] The computing unit 30 is configured for implementing calculation and / or processing steps, in particular for implementing process steps according to this disclosure. By convention, as shown in [Fig. 5] depicting an array 10 of transducers 11 on a surface of a medium M, a spatial coordinate system of the medium M is defined by taking a first X-axis and a second Z-axis perpendicular to it. For simplification, the first X-axis corresponds to the transverse direction in which the transducers 11 are aligned in the example of a linear array, and the second Z-axis corresponds to the depth of the medium M relative to this array 10 of transducers 11.This definition can be adapted to the context and thus, for example, extended to a three-axis spatial frame in the case of a two-dimensional array, or to a polar frame in the case of a curved array, or to any other frame adapted to and / or dependent on the structure and shape of the ultrasonic transducer array. Therefore, in the remainder of this disclosure, we will use a Cartesian XZ frame, corresponding to a linear probe, for the sake of simplicity in the explanations, but a specialist in the field could easily generalize and apply the results to any type of frame.

[0176] In the remainder of this disclosure, reference is made to an array 10 of transducers 11 for transmission and reception, it being understood that, more generally, several transducer arrays may be used simultaneously. The transducers 11 may be both transmitters and receivers, or only transmitters for some and only receivers for others. Similarly, an array 10 may consist of one (1) to N transducers 11, of the same type or of different types.

[0177] The array 10 of transducers 11 serves, for example, both as a transmitter and as a receiver, or is made up of several sub-arrays of transducers, some being dedicated to the transmission, others to the reception of ultrasonic waves. By array of transducers is meant at least one transducer, an aligned or non-aligned sequence of transducers, or a two-dimensional distribution of transducers (for example, a matrix of transducers), or any spatial distribution of transducers.

[0178] Where reference is made in this disclosure to calculation or processing steps for the implementation of process steps, it is understood that each calculation or processing step may be implemented by software, hardware, firmware, microcode, or any suitable combination of these or related technologies. When software is used, each computational or processing step can be implemented by computer program instructions or code that can be interpreted or executed. These instructions can be stored or transmitted to a storage medium readable by a computer (or computing unit) and / or executed by a computer (or computing unit) to implement these computational or processing steps.

[0179]

[0180] Figure 3 shows a diagram of the main steps according to the invention. A step a) of acquiring a series of canonical reflection matrices Rui(t, μ) is distinguished. In b), a set of focused reflection matrices Rxx is determined. (z, #m) of the medium by a focusing process for several points of transverse spatial position (x, y) of a region and axial position z=c0t / 2 from the canonical reflection matrix R ui(t, #m) for a sound velocity model cO,

[0181]

[0182]

[0183]

[0184] In step c) we construct dynamic focused reflection matrices rIL # by removing the static component. Step d) allows the calculation of correction laws z) for each pixel of the image from the dual reflection matrices, the latter being obtained by projection at the output of the dynamic focused reflection matrices R^^ #m) into a correction basis c of the aberrations, In step e), we determine corrected focused reflection matrices R^^ # In step f), we determine a dynamic confocal signal S c (x,z- # m) for every point of spatial position (x, z).

[0185] Then construction in step g) of an image from the dynamic confocal signals.

[0186] These steps are described in more detail below.

[0187]

[0188] According to the invention, the ultrasonic construction method implemented by the computing unit 30 of system 1 comprises a series of acquisitions of reflection matrices at different times. Each realization of the reflection matrix will be called a frame hereafter, and the i-th frame will be denoted #m. Each reflection matrix can be acquired in the following manner:

[0189] - a step of generating a series of incident ultrasonic waves USin in a zone of said medium, by means of an array 10 of transducers 11, said series of incident ultrasonic waves being an emission basis i; and

[0190] - for each emitted wave i in, the field reflected by the medium is measured by each transducer and is denoted 7?(u out,im), where t is the echo time and the vector u out marks the position of each transducer. Each field is stored in a series of canonical reflection matrices R ui(t, #m)=[Æ(u „„t,i imt,# ,„)] defined between the transmitting basis i at the input and a receiving basis u at the output.

[0191] One possible method for measuring this canonical reflection matrix is ​​to successively emit an ultrasonic pulse from each transducer of the array, whose position is located by the coordinate uin, as shown schematically in [Fig. 4](a). This results in a diverging cylindrical (or spherical) incident wave. This wave is reflected by the scatterers of the medium, and the backscattered field is recorded by each transducer as a function of time, as shown in [Fig. 4](b). By repeating this operation with each transducer used successively as a source, the canonical reflection matrix Ruu(t, #m) expressed in the basis of the transducers is determined. This matrix is ​​composed of all the impulse responses A*(uout, uinj, #m) between each transducer. This matrix is ​​thus rich in information about the medium under study. However, the method assumes that the medium remains stationary throughout the measurement period.Furthermore, the recorded signals have a poor signal-to-noise ratio because the medium is insonified by a single transducer.

[0192] A second way of constructing this canonical reflection matrix consists of insonifying the medium with a plane wave series basis. This method overcomes the previous problems. [Fig. 4](c) illustrates the principle of this plane wave illumination. A delay law r' is applied to each signal at the emission stage to form a wavefront inclined at an angle θin with respect to the transducer array. At the reception stage, illustrated in [Fig. 4](d), the field backscattered by the medium, θ(uout, θin, t, μm), is measured by all the position sensors uout for each incident plane wave θin. All of these responses form a canonical reflection matrix Ruü(t, #,„)=[ 7?(u out, 0 in, t, #,„)]. This method gave rise to ultrafast imaging and elastography, and is described, for example, in the document:

[0193] “Coherent plane-wave compounding for very high frame rate ultrasound and transient elastography”, G. Montaldo et al. (IEEE Trans. Ultrason., Ferroelect. Freq. Control 56 489-506, 2009).

[0194] A third way to create this canonical reflection matrix is ​​to insonify the medium with a basis of diverging waves, as shown in [Fig.4](e) and [Fig.4](f), which allows the acoustic field to be illuminated more broadly than by using plane waves. This basis is identified by the position s in the virtual source associated with each diverging wave. This technique, used particularly in super-resolution imaging, is explained in the document:

[0195] “Ultrafast imaging of the heart using Circular Wave Synthetic Imaging with Phased Arrays”, Couade et al., IEEE International Ultrasonics Symposium (2009).

[0196] Each canonical reflection matrix R ui(t,# m) recorded can be a "real" matrix, that is, composed of real coefficients in the time domain, the electrical signals recorded by each of the transducers being real numbers. Alternatively, this matrix can be a "complex" matrix, that is, composed of complex values, for example in the case of demodulation for in-phase and quadrature beamforming (known in English as "beamforming IQ").

[0197]

[0198] Focusing of the reflection matrix

[0199] After or in parallel with the acquisition sequence, a channel-forming process is applied independently at the input and output of the measured reflection matrices. The result is a series of focused reflection matrices R xx (z,#m) which includes responses R(x out, x in,z,#m) of the midpoint between a virtual input transducer TVin of spatial position (x in,z) and a virtual output transducer TVout of spatial position (x out,z).

[0200] The responses of the focused reflection matrix R xx(z,#m) correspond to an acoustic pressure field calculated between all points in the middle of lateral positions x in and x out, located at the expected depth z=c0t / 2 and at an echo time t, and for a speed of sound assumed c0.

[0201] In the step of determining the focused reflection matrix R xx(z,# m), we apply with reference to [Fig.5]:

[0202] - an input focusing process from each reflection matrix canonical R ui(t,# m) which uses a time of flight on the way forward of the waves between the emission base i and the virtual input transducer TVin and which creates a focal spot called the input spot around the first point PI of spatial position r in=(x in,z), said input focal spot corresponding to the virtual input transducer TVin,

[0203] - an output focusing process from the canonical reflection matrix R m(t,# m) which uses a return time of the waves between the virtual output transducer TVout and the receiving base transducers u and which creates a focal spot called output around the second point P2 of spatial position r out =(x out,z), said focal spot output corresponding to the virtual output transducer TVout.

[0204] These input and output focusing processes form an input-output focusing process, referred to in the remainder of this disclosure as a focusing process or simply focusing.

[0205] In other words, in this ultrasonic construction method, the virtual input transducer TVin corresponds to an ultrasonic "virtual source" located at the spatial position r in the medium, and the virtual output transducer TVout corresponds to an ultrasonic "virtual sensor" located at the spatial position r out. This virtual source and sensor are spatially separated by the difference in their spatial positions Ax = x out - X in.

[0206] The maximum distance Ax max between x out and x in is dictated by the number of illuminations used to acquire the reflection matrix. In the plane wave basis, for example, the 80° angular sampling of the illumination sequence imposes Ax max -2 / (250) to prevent the focused reflection matrix from being polluted by grating lobes. Ideally, Ax is chosen according to the level of aberrations in order to limit the amount of data to be acquired and stored in memory. It is typically on the order of a few millimeters in ultrasound imaging.

[0207] The expected depth of virtual transducers is the parameter z used in the focusing law for a sound velocity model c0. Their actual depth is dictated by the axial position (in depth) of the isochronous volume, i.e., by the echo time t and by the sound velocity distribution c(r) in the medium. The lateral dimension of the virtual transducers is dictated by the focal spot produced by focusing at this actual depth.

[0208] Each focused reflection matrix R xx(z,# m) can be determined or calculated:

[0209] - either in the time domain, in which case it can be explicitly denoted with the time parameter t, that is, denoted R xx(z,t,#m), and

[0210] In practice, the data of the focused reflection matrix R xx(z,t,# m) are calculated between two predetermined time instants; and

[0211] In practice, it can only be considered at the expected ballistic time (i.e. at t=0 in the focused basis)

[0212] - either in the frequency domain, in which case it can be explicitly denoted with the angular frequency parameter w which corresponds to a frequency f by ai = 2æ / , that is to say denoted R xx(z, a',# m), and

[0213] In practice the data of the focused reflection matrix R xx(z, m) are calculated between two pulsations of a frequency bandwidth, for example between a lower pulsation and an upper pulsation, for a central pulsation 'A.

[0214] Thus, the subsequent calculations of the process can be carried out in the time domain or in the frequency domain.

[0215]

[0216] In the first time-domain calculation, each focused reflection matrix R xx(z,t,# m), in the midpoint between the virtual input transducer TVin and the virtual output transducer TVout, is obtained by focusing using a channel-forming calculation at the inputs and outputs. The coefficients of this focused reflection matrix R xx(z,t,# m) can be determined by: [02 17] # m) — (L- '■ ) ^ ( L* + G / j ( Xoup L j # m)

[0218] in which:

[0219] Nin is a first normalization coefficient,

[0220] Nout is a second normalization coefficient,

[0221] R ui(t, # m) is the reflection matrix, of which each coefficient R(u out, i in, t, # m) is the field recorded by the spatial position transducer u Ou following the emission of index i in in the emission basis (i) and at time t, to which the delay times rin and rout have been applied;

[0222] Ain ( ïifP xinz ) and Aout ( Zout, z) are apodization coefficients which are predefined, for example to keep a constant digital aperture for both transmission and reception;

[0223] The first normalization coefficient Nin can, for example, be defined by: N- ( X-1 Z ) = A'h ( î / n' X- , Z ) similarly, the second normalization coefficient Nout can be defined by: M Jr 71= P 4 Au v 1 1 l'ont\AotO' uou& A oufr )

[0224] Tin (XilP Z ) is the expected time of flight for each incident wave to reach the first focal point of spatial position (xin, z ) in a model sound velocity medium c0

[0225] Tout^out' xout- z) is the expected time of flight for a wave reflected from the second spatial position focal point (xout, z) to the position transducer

[0226] These delay times Tin and Tout are usually calculated by a person skilled in the art from an established model of the speed of sound. A relatively simplifying assumption is to assume a homogeneous medium with a constant speed of sound c0 in the medium. In this case, the times of flight are directly obtained from the distances between the probe transducers and the virtual transducers. Thus, these delay time calculations depend on the type of wave, the assumed speed of sound, and the geometry of the transducer array.

[0227] For example, in the particular case of a plane wave with an emission angle:

[0228] - the emission delay time ^in can be obtained by:

[0229]

[0230]

[0231]

[0232]

[0233]

[0234]

[0235]

[0236]

[0237]

[0238]

[0239]

[0240]

[0241]

[0242]

[0243]

[0244]

[0245]

[0246]

[0247] . Xi„smf)in+zincos0;a ^in ( ^iw % ) ~ - The delay time at reception 'out' can be obtained by: / \ out \ — co Thus, these examples of delay time calculations clearly show that they are a function of the type of wave and the speed of sound, assumed here to be constant in the medium. The number of elements in the Nin transmission base is, for example, greater than or equal to one (1), and advantageously greater than or equal to two (2). The number of elements in the Nout reception base is, for example, greater than or equal to two (2). Finally, each focused reflection matrix R xx(z,t,# m) expressed in the time domain can be transformed in the frequency domain into a focused reflection matrix R xx(z, # m) by a Fourier transform, that is: Rxx(z, ox #m') = i^dt Rxx(z, t, # m) This Fourier transform can be implemented by any type of discrete Fourier transform, normalized or not. In the second case of calculation in the frequency domain, each canonical reflection matrix R ui(t, # m) expressed in the time domain, since it is made up of the signals received by the transducers, can be transformed in the frequency domain into a canonical reflection matrix RMÎ( ux # m)) by a Fourier transform, that is to say by: Rai(W, #m) = fdt R^(t, #m)e-^ This Fourier transform can be implemented by any type of discrete Fourier transform, normalized or not. Thus, the focused reflection matrix R„(z, # m) of the medium can be obtained by focusing by the matrix calculation below, essentially equivalent to focusing by temporal path formation explained previously, that is to say by the following matrix product: Rjrx ( m) — ^ux ( ) Rj^i ( <a # x ( z-, )="R(x,X" dans lequel :

[0248] la matrice üj, m) est transformée de former chaque réflexion canonique r^.(t, m),

[0249] g^fz, w) passage réception adaptée pour le base (u) à focalisée (x) profondeur z et pulsation a’,

[0250] z, d’émission (i)

[0251] les symboles *et f désignent respectivement les opérations matricielles conjugaison transposition-conjugaison.

[0252] le symbole x désigne un produit matriciel.

[0253]

[0254] extraction composante dynamique

[0255] afin d’extraire signaux des diffuseurs en mouvement, une étape filtrage statique effectuée. la séparation composantes peut être réalisée diverses manières.

[0256] suivant premier mode réalisation, isolée soustrayant au signal moyenne glissante sur n="5)." frames (typiquement coefficients r sont ainsi obtenues

[0257] ( x, x, f, m>Z, f, # m) - Wx, X, Z, f.

[0258] ^R^C^ f # ) — [ / ^(x XZ f $ )] is the dynamic component of the Reflection matrix. Note that the same filter can be applied upstream to the measured matrix Rui(t, # m), which would be computationally more advantageous. However, the benefit of applying it to the focused reflection matrix is ​​the ability to adjust the number N of frames over which the moving average is performed.

[0259] According to a second embodiment, a more sophisticated high-pass (or band-pass) filter is applied along the dimension #m of the reflection matrix in order to filter the static component of the reflection matrix.

[0260] According to a third embodiment, the dynamic component is isolated by performing a singular value decomposition of the focused reflection matrix rearranged in a two-dimensional form as follows: R - [R({x, x, z, f}, #,„)]• The singular value decomposition of this matrix is ​​written: p _V v YYR# - jLpXp^pX P

[0261] with Xp - [ U p( [x, x, z, f}, f ) ] corresponding to the singular vectors of the reflection matrix in the associated space {x, x, z, f}.

[0262] Yp=[VP(# m ) ] corresponding to the singular vectors of the focused reflection matrix R# in the frame basis,

[0263] Xp corresponding to the real and positive singular values ​​of the focused reflection matrix R# arranged in descending order: > ^2 > ' ' ' >

[0264] The first P eigenspaces associated with the highest singular values ​​(■^^max) are associated with the static component. The last Q eigenspaces associated with the lowest singular values ​​(Xp < Zmjn) are associated with the noise. The eigenspaces associated with the intermediate singular values ​​(Anin < Xp < ^max) are associated with the dynamic signal of interest:

[0265] K# - p

[0266]

[0267] The thresholds / mjn and Xtiiiix can be determined as inflection points of the distribution of singular values.

[0268]

[0269] For the sake of simplicity of notation, the dynamic component of the reflection matrix will be denoted R and not Rd in the following.

[0270] Figure 6 illustrates the extraction of the dynamic component of signals from data acquired on a sheep brain. The different representations of Figure 6 allow comparison of the evolution of the reflectivity of a point in a speckle zone of the medium over time before filtering, in image A, and that of the same point after filtering, image B. The temporal evolution of the complex reflectivity associated with the pixel designated by the white cross [A and B] is represented using a point cloud in the complex plane [C] before and after filtering of the dynamic part of the signals, respectively at the periphery and at the center. Image D is an enlargement of the central part of the complex plane, highlighting the quasi-random nature of the complex reflectivity after filtering.

[0271]

[0272] Method for correcting aberrations

[0273] The process aims to correct aberrations, these aberrations being for example due to variations in structures in the medium which induce variations in the speed of sound and variations in reflectivity.

[0274] The method includes a correction process comprising steps of:

[0275] - determining a set of focusing laws from the responses of the medium obtained for the different realizations of the speckle measured at times # m.

[0276] - determination of corrected reflection matrices R'xx ( z, t, # m ) of the medium by applying the correction law to the measured reflection matrices R.xx ( Z, t, # m ) .

[0277] Thanks to these provisions, the method advantageously allows for local probing of the medium and correction of the focused reflection matrix with respect to aberrations, in particular by determining a correction law for each point in the middle and, optionally, for each frequency of the ultrasonic wave.

[0278] In addition, the process may include a step of:

[0279] - determination of an intensity E of an ultrasound image point corresponding to a point of spatial position r = (x, z) from the diagonals of the corrected reflection matrices R'M(z, A #m)

[0280] - Determination of the intensity of a point in a power Doppler image by summing the intensity of the ultrasound image at the same point over all or part of the frames #m

[0281] - determination of a colored Doppler image allowing access to the directionality of movement in tissues by performing a Fourier transform along the temporal dimension #m of the complex amplitude of the ultrasound image as described in the article:

[0282] E. Mace, G. Montaldo, B. -F. Osmanski, I. Cohen, M. Fink and M. Tanter, "Functional ultrasound imaging of the brain: theory and basic principles," in IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 60, no. 3, pp. 492-506, March 2013,

[0283]

[0284]

[0285] Determination of the Law of Focusing®

[0286] The step of determining a correction law ® then includes substeps carried out at each depth z and each frequency / , of:

[0287] - determination of a set of dual reflection matrices Rcx(z, / # m) by forward projection of the focused reflection matrix Rxx(z, / # m) to a correction basis (c),

[0288] - calculation of the correction law z) from the dual reflection matrix RMf # m), said correction law being determined on the basis of correction (c), x)], so that said correction law 0(x, z) is a spatio-frequency correction law,

[0289] - determination of a corrected dual reflection matrix R^. around the reference point and whose coefficients are written according to #m)= [Rc(xcz f # m ) ] ' determined by performing the term-by-term product between the dual reflection matrix Rcx(z, / # m) and the phase conjugate of the correction law (p[x, z], that is, by: r'x = Rcxo 0*

[0290] where

[0291] The symbol * denotes a phase conjugation operation

[0292] the symbol 0 is the Hadamard product, such that:

[0293] rcx(X, c, Z, f, #m)= Rcx(xc, z, f, # m) &(xc, z, f) •

[0294] The final step of the process includes determining a set of focused reflection matrices corrected by back projection of the matrices of corrected dual reflections R^ (z, / # m) towards the focused basis (x).

[0295] Thanks to these provisions, the method advantageously allows local probing of the medium and correction of the focused reflection matrix with respect to aberrations, in particular by determining a focusing law for each point of the medium, and for each frequency of the ultrasonic wave.

[0296] This correction is performed in a correction basis c adapted to the aberrations to be corrected. The correction basis c is an input correction basis or an output correction basis.

[0297] Basic examples of correction are:

[0298] - a plane wave basis or spatial Fourier basis,

[0299] - a base of u transducers,

[0300] - a basis corresponding to the supposed location of the aberrators in the medium, that is to say for example a plane between the transducer plane (transducer base u) and the focusing plane (focusing base x),

[0301] - a basis corresponding to a plan determined by optimization, for example by a correlation matrix whose first eigenvalue is maximal.

[0302]

[0303] Dual reflection matrix Rcx(z, / # m)

[0304] According to one embodiment of the process of this disclosure, the forward projection allows the determination of a dual reflection matrix Rc(z, / # m)- This forward projection can be performed by:

[0305] a matrix product between a change-of-basis matrix P and the focused reflection matrix Rxx(z, / that is: R^, / , # m ) = P(z, / ) x R. Jz, f , # J

[0306] where:P(z,f ) = [P(c, x, z, f) ] is the change-of-basis matrix at each frequency / between the focused basis (x) at depth z and the correction basis (c).

[0307]

[0308] The change-of-basis matrix P depends on the correction basis c used.

[0309]

[0310] In the case of a correction basis corresponding to a plane wave basis (c=k), the change-of-basis matrix P is the Fourier transform operator.

[0311] In the case of a linear transducer array 10 for generating a two-dimensional image, the coefficients of this transition matrix P can be written as follows: P(fc„i,z / ) = P{kx,x) =exp{-ikxx)

[0312] where kx, the transverse component of the wave vector k associated with each plane wave.

[0313] In the case of a matrix-type array 10 of transducers 11 for generating a three-dimensional image, the coefficients of this pass-through matrix P can be written as: P(k^p, z, f) = P(p, Z, f) = exp(-ik^p)

[0314] where

[0315] is the transverse component of the wave vector k associated with each plane wave, and

[0316] p = (x, y), the transverse position vector.

[0317]

[0318] In the case of a correction basis corresponding to a basis of the transducers (c=u), the coefficients of the transition matrix P correspond to the normal derivative of the Green's function relating each focal point of spatial position (x, z) and each transducer of spatial position (u, 0).

[0319] In the case of a linear transducer array for generating a two-dimensional image, the coefficients of the transition matrix P can be written as follows:

[0320] where V is the gradient projected along the depth direction z, and

[0321] G2D(u, r) is the 2D Green's function that relates each transducer u = 0 to each point r of the midpoint M, with: r, / ) = - { Ho ( k01U - r I )

[0322] where = 2zr / / c0 is the wave number,

[0323] Ho is the first-order Hankel function whose asymptotic expression is: u-ri / C(J)

[0324] In the case of a matrix-type transducer network for generating a three-dimensional image, the coefficients of this pass-through matrix P can be written as: P(li, X, z, f) = VzG3n(u,r, / )

[0325] where G3D(u, r) is the 2D Green's function that connects each transducer U_Q^ to each point r — (p, z) of the midpoint M, with: G. n ( u, r, / ) =----..... ■WJ ' 4-^Up|*+z2

[0326] The coefficients of the change-of-basis matrix P are therefore written in this case as follows: ex^-ik^-f^z2) P(u,x, z, f) = -¾-----

[0327]

[0328] Correlation matrix

[0329] According to a first embodiment, this calculation step comprises:

[0330] - the construction of a correlation matrix C(x,z) from the reflection matrix dual Rcx(z,f,, # m) for each point (x,z) in the field of view

[0331] - the analysis of this correlation matrix C(x,z) to determine the law of spatio-frequency focusing d>(x,z) for each point of the field of vision.

[0332]

[0333] According to a first variant, the correlation matrix C(x,z) is determined in the correction basis c and in the frequency domain, by the following calculation of the elements of the correlation matrix C = C cc: C({c, / }, = £# R(x,c,z,f, #m) Rv(x,c',zf\ #m)

[0334] Where * is the conjugation operator.

[0335] This operation allows the reflected fields in the correction basis to be correlated for each virtual source in (x,z). This correlation is averaged over the different frames # m (i.e., different realizations of the speckle) in order to eliminate the random reflectivity of the medium and thus synthesize a coherent guide star from the different realizations # m of the speckle.

[0336]

[0337] According to a second variant, the correlation matrix C is determined in the basis of frames # m, by the following calculation of the elements of the correction matrix C = C##; C( #m, # bx, z) — (x, c, z, f, # m) x R\x, c, z, f, # i)

[0338] * is the conjugation operator

[0339]

[0340] Focusing law estimation

[0341] According to a first variant, the analysis of the correlation matrix C(x,z) is carried out by an eigenvalue decomposition of the correlation matrix C(x,z), and the spatio-frequency correction law ^(x^) is the first eigenvector U] of the correlation matrix C(x,z) in the correction basis (c), that is C = C cc.

[0342] Since the correlation matrix is ​​Hermitian (ç _ q"), its eigenvalues ​​are real. and positive.

[0343] The correlation matrix C cc(x,z) can thus be written C _ Y tt Tri

[0344] or in terms of matrix coefficients: = Yp&püp(c, f )Up(c', f')

[0345] with Up corresponding to the eigenvectors of the correlation matrix C

[0346] corresponding to the real and positive eigenvalues ​​of the correlation matrix C cc(x,z) arranged in descending order; >cfv

[0347] We then have the spatio-frequency correction law z) which is equal to the first eigenvector, i.e. • ®(xz) — Uj; or its standardized version, O(x, z) = exp(jarg{U1} ), i.e. a spatio-frequency correction law of which The coefficients have a unit amplitude but a phase equal to that of U (the symbol arg{X] denotes the phase of the vector X; or to an inverse filter-type correction, z) = exp(arg{Uj}) / | The first option is preferable if there is a poor signal-to-noise ratio (matched filter). In general, however, the second option is preferred so that the correction does not act as an amplitude filter but only corrects phase distortions. Finally, the third option is relevant when the aberrant medium inhomogeneously attenuates certain components and / or frequencies of the field that we wish to enhance in order to obtain a more accurate estimator of the reflectivity.

[0348]

[0349] According to a second variant, the analysis of the correlation matrix C cc(x,z) is carried out by a singular value decomposition of the dual reflection matrix rearranged as follows: Rjx? z) = [fl( {c? / }, # x, Z)]

[0350] The eigenvalue decomposition of the correlation matrix Ccc(x,z) carried out in the first variant is indeed equivalent to the singular value decomposition (SVD) of each dual reflection matrix Rc#

[0351] Singular value decomposition applies to rectangular matrices, and when applied to the dual reflection matrix Rc# at each point (x,z), it is written as follows: ■» _ y ,C~tt W

[0352] or in terms of matrix coefficients: R( p(c, fWp( #m)

[0353] with Up = [ Up(c, f) ] corresponding to the singular vectors of the dual freflection matrix R ^x in the correction basis, or equivalently, to the eigenvectors of the matrix C Cc as defined in the first variant.

[0354] Vp=[Vp(# m ) ] corresponding to the singular vectors of the dual reflection matrix R(. in the frame basis,

[0355]

[0356] 2^ corresponding to the singular values ​​of the dual reflection matrix Rc# which are, by definition, equal to the square root of the eigenvalues ​​aP of the correlation matrix C as defined in the first variant: &p =

[0357] We then have the spatio-frequency correction law <J)(x, z) qui est égale au premier vecteur singulier de la matrice de distorsion duale Rc#, i.e. (x, z) = ; ou à sa normalized version, O(x, z) = exp( / arg{U|] ). i.e., a spatio-frequency correction law whose coefficients have unit amplitude but whose phase is equal to that of U, (the symbol arg{X] denotes the phase of the vector X); or to an inverse filter type correction, z ) exp(} ) / |UJ •

[0358] The advantage of the singular value decomposition of the dual reflection matrix Rc#, compared to an eigenvalue decomposition of the correlation matrix C cc, is the speed of computation of the numerical algorithms of the

[0359]

[0360]

[0361]

[0362]

[0363]

[0364]

[0365]

[0366]

[0367]

[0368]

[0369]

[0370]

[0371]

[0372] decomposition into singular values This search for the spatio-frequency correction law z) is also equivalent to solving the following equation: a <lXx, z, ) = C^x, z) x <$(x, z) où x est le produit matriciel et a est une constante iteratively using the following expression, which corresponds to an iterative time-reversal calculation: O„+1(x, z) = Ccc(x, z) with ¢0 an arbitrary wavefront, for example — [ [ ] ]1 Therefore, the spatio-frequency correction law ^x, z) is obtained by: Oix, z)= lirn<î>,,(x. zl' or its standardized version: z) = exp ( jarg {UmO^x, z) ] ) ' or its inverse filter version: .. X ... exp(jarg^^xx)}), p The iterative time-reversal algorithm converges to the same first eigenvector U of the matrix Ccc. In practice, there may be an advantage to using an iterative time-reversal algorithm rather than a SVD because it can converge after a few iterations, resulting in faster computation.

[0373]

[0374] According to a third variant, the analysis of the correlation matrix Ccc is carried out by solving the following equation: 0 ( x, z ) = exp ( j arg {Ccc(x, z) x <D ( x, z )} )

[0375] iteratively by the following expression, which corresponds to an iterative phase-reversal calculation: z ) = exp ( j arg {Ccc(x, z) X ( X, Z )} )

[0376] where x is the matrix product,

[0377] with: <î>0 an arbitrary wavefront, for example — [ ] j ]T.

[0378] Then, the spatio-frequency correction law z) is obtained by:

[0379] 0 (x, z) - limO„ (x. z),

[0380] Or its inverse filter version:

[0381] , Q(x. z) - ;----------tt

[0382] The advantage of an iterative phase reversal algorithm compared to the previous alternatives is that it is a more reliable estimator of the phase of the correction law <P(x, z) et donc d’accéder in fine à une meilleure compensation des distorsions de phase induites par l’aberrateur.

[0383]

[0384] According to a fourth variant, the analysis of the correlation matrix C## is carried out by solving the following equation: W(x, z) = exp(j arg{C##(x, z) xW(x,z)})

[0385] where x is the matrix product,

[0386] iteratively by the following expression: W„+i(x,z) =exp(j arg{C##(x,z) xW„(x,z)})

[0387] where x is the matrix product,

[0388] with Wo an arbitrary wavefront, for example Wo = [ 1 1 ]T

[0389] which allows us to obtain the following vector W(x, z): W(x, z) = limWH(x, z)

[0390] This vector W(x, z) = [ W( # m X, z)] defined in the frame basis contains the phase of each inconsistent guide star synthesized by focusing at point (x,z) for each frame # m.

[0391] The phase conjugate of this vector W(x, z) can then be used to rephase each incoherent virtual star so that they can be recombined coherently and thus obtain an estimator of the spatio-frequency correction law ^x. z) unbiased by the random reflectivity of the medium. Mathematically, this operation is written as follows: <tfic,f,x,z) ^exp{jxarg\^xsR(xc,z, f, # ^W^x, z, #J}}

[0392] The advantage of this approach compared to an SVD of the distortion matrix R c# (second variant) or the iterative phase reversal algorithm (third variant) is to converge towards a correction law not biased by the larger amplitude of the ultrasonic signal on certain frame of the reflection matrix (large amplitude caused by the passage of a bright scatterer such as a bubble or an experimental problem).

[0393] Figure 7 illustrates an example of aberration law extraction by iterative phase reversal in the dynamic speckle. [A] The aberrated focused reflection matrix ~Rpp(z, #) [B] is projected at output into a correction basis (c), shown here by the base of the transducers (u). p - (x, y) is the transverse position vector. Each frame of the dual matrix # ) constitutes a dynamic realization of disorder. [C] This matrix is ​​rearranged along dimensions (c) and (#) to construct an n / \ matrix at each point (j) and calculate [D] the correlation matrix ■**■ / *43 P • ? ZI \£ f m / \ m / associated p ( Y An analysis by RPI (iterative phase reversal) of the ccviri v matrix \ allows us to extract from these wavefronts an estimate J \ laws of aberration at each point of the field of vision. The phase conjugate of the laws of aberration (V) is used as the adaptive focusing law Vin' v in reception for each frame, allowing the recovery of a corrected output focal spot. Although illustrated here in basis (c) = (u), the reasoning can be extended to another correction basis, in particular that of plane waves (k). This approach can be described by constructing a distortion matrix or schematized as in [Fig. 7] using the dual reflection matrix. The two approaches are equivalent: in the first case, an aberration law is extracted, while in the second, a focusing law (aberration and geometric curvature) is extracted.

[0394] Figure 8 illustrates examples of local transverse aberration laws extracted from the dynamic speckle. [A] Several regions of the field of view are chosen for illustrative purposes, and their corresponding transverse aberration laws are shown. The pairing is done by color code. The transparent outline of each area in image A defines the area on which the law is estimated, while the opaque area in the center of image A delimits the area on which this law will be applied during correction. These areas also have an extension along the y-axis, which is not shown here. The Strehl ratio "S" associated with each of the aberration laws in images B is indicated.

[0395]

[0396] Correction of aberrations in reflection matrices

[0397] From the obtained focusing law, aberrations will be compensated in post-processing by recalculating dual and focused reflection matrices from the estimated focusing law.

[0398] Corrected dual reflection matrix

[0399] The corrected dual reflection matrix xf Z # J is determined by performing the term-by-term product between the dual reflection matrix Rc(z, f) and the phase conjugate of the focusing law®, that is, by: R^ = RCX° ®*

[0400] where the exponent * represents the phase conjugation operation

[0401] The symbol ° is the Hadamard product, that is, the term-by-term matrix product of matrix coefficients, such that

[0402] R'(x,c,f, z,, #tn) = R(x,c,f,z, # m)(p' (x,c, f,z)-

[0403]

[0404] Corrected focused reflection matrix

[0405]

[0406] A corrected focused reflection matrix R^J^ y) is then determined by back projection of the corrected dual reflection matrix (z,f) to the focused basis (x).

[0407] The back projection is performed by a matrix product between the transition matrix P defined above and the focused reflection matrix r'x (z,f), that is: Rxx (z,f, #m)=P'(z,f)* iCfc f, #ni

[0408] Where the exponent t denotes the matrix operation of trans-conjugation.

[0409]

[0410] Correction processing iterations

[0411] According to the embodiment of the process, the steps of the correction treatment, that is to say the steps of:

[0412] - of determining the correction law O(x, z) said step comprising possibly the steps of determining the dual reflection matrix Rcx(z, f, # m), S160 of calculating the correction law d)(x, z) and of determining the corrected dual reflection matrix (zf # m ) 'Ct

[0413] - of determination of the corrected focused reflection matrix R^J^ #

[0414] are iterated several times (twice or more than twice).

[0415] At each iteration, the forward projection uses the corrected focused reflection matrix Rxx(z- # m) °obtained during the back projection of the previous iteration instead of the focused reflection matrix Rxx(z, # m)-

[0416] Thus, at each iteration, the spatio-frequency correction law is improved to take into account one or more aberrations of the environment better and better.

[0417]

[0418] According to a first variant of this iterative process, at each iteration of the step of determining the dual reflection matrix Rc x(z, / , # Ht), a different correction basis c is used, for example to correct different aberrations located in different places in the medium.

[0419] For example, the medium can be discretized or modeled by a succession of layers along the depth direction z, and the correction bases c of the iterations correspond to planes of these successive layers. In other words, we apply to the course of iterations of corrections corresponding to a plurality of aberrations in the environment.

[0420]

[0421] According to a second variant of this iterative process, at each iteration of the step of determining the dual reflection matrix, a forward projection is used either towards an input correction basis or towards an output correction basis of the reflection matrix. In the latter case, the projection of the matrix Rxx(zJ, #, H) onto the correction basis is performed as follows:

[0422] R„ (z, / .#„)= ?(z, / ) x R Jz. f. # J

[0423] where the symbol refers to the matrix transposition operation.

[0424] In the succession of iterations, one can alternate between using an input correction basis and an output correction basis. Thus, the spatio-temporal correction law O is improved at each iteration, and the correction of aberrations is improved.

[0425]

[0426] Variants of the aberration correction process

[0427] Depending on the experimental conditions, the aberration process can be simplified or made more complex according to several variants described below.

[0428]

[0429] Spatial compensation of aberrations

[0430] If the medium upstream of the focal plane does not induce temporal dispersion of the echoes (absence of reverberations, negligible multiple scattering), there is not necessarily any advantage in considering the frequency components of the confocal signals independently. In this case, one can consider the broadband focused reflection matrix,

[0431] {f+ ^(z,t = O, #m) #m)

[0432] as the starting point for the aberration correction process. The projection of the reflection matrix onto the correction basis is then performed by considering the propagator at the center frequency. The rest of the process is identical to that described above. Only the frequency dependence of the various expressed quantities is no longer relevant.

[0433]

[0434] Temporal compensation for aberrations

[0435] If, on the contrary, the medium induces a primarily temporal dispersion of the ultrasonic echoes and few spatial aberrations, there is no point in considering the focused reflection matrix but only its confocal signal:

[0436] S(x,z,f, #m) - R(x,x,z,f,#m)

[0437] The correlation matrix to be considered in this case for obtaining the correction law, 0 ( X) z ) = [O ( f, x, z )]> is given by:

[0438]

[0439]

[0440] C( / , / ',X,z) = E #f S'(xZ, / , # m ) S"'(x,zf, ^ m ) Exploitation of the spatial memory effect

[0441] The number of frames may be insufficient for a correct estimation of the focusing law. In this case, instead of the dual reflection matrix, one can consider the associated distortion matrix, Dxc(z, f, #,„) = [ Dxc(x, c, z, f, # m) ], and average the correlation matrix not only over the different frames but also over adjacent pixels belonging to the same isoplanetary patch ^p: C ( { C, f}, { C \ f' ] , Xp, Zp ) = D' ( X c \ z, f, # m )

[0442] Withy, # m) = R(x,c,zJ, #m)Rpef(x, c, zf)

[0443] and Rref = [ Rref (x,cy zf) ], a reference matrix of a model medium in in which the speed of sound is c, the expected speed of sound for the medium, and in which a plane reflector is positioned at depth z. Averaging the correlation matrix over several speckle grains can accelerate the convergence of the correlation matrix to the associated covariance matrix. This results in a more accurate estimation towards the focusing law.

[0444]

[0445] Construction of images of the environment

[0446] Dynamic confocal image

[0447] According to one embodiment of the ultrasonic construction process of this disclosure, the process further comprises a step of:

[0448] - determination of the dynamic confocal signal S c (x,z, # m) of a position point spatial (x, z), from the diagonal coefficients of the corrected focused reflection matrix R^x(z, # m) integrated over the bandwidth of the ultrasonic signals, that is, by combining the confocal responses of the point at several frequencies. We thus have, for example, the following calculation: g çx z # ) — xzf # )

[0449] The intensity of this quantity,

[0450] h^z #m}^\Sc(x,z #M)|2'

[0451] determined in a plurality of points (x,z) allows to construct a corrected confocal image of the medium, which corresponds to a classic ultrasound image free from the problems of aberrations, reverberations and frequency dispersion of the speed of sound in the medium studied.

[0452] Figure 9 illustrates a comparison of dynamic confocal images before and after correction. Several examples of original confocal images (left) and their corrected equivalents (right) are shown for different frames, with or without microbubbles [AC] and without [D]. The benefits of the method are particularly noticeable in the microbubble image (bright dots), which appear better resolved and more contrasted after correction.

[0453]

[0454] Power Doppler image

[0455] According to one embodiment of the ultrasonic construction process of this disclosure, the process further comprises a step of:

[0456] - determination of a 1P (x,z) Power Doppler image of a position point spatial (x, z), from the sum over the different frames of the dynamic confocal intensity Sc (x, z, #m) measured at each point (x,z). For example, we thus have the following calculation: I (x, z) = Ic(x, z, #m) * X Ul /

[0457] The previous intensity determined at a plurality of points allows a corrected power Doppler image of the medium to be constructed, which quantifies the movement of the tissues during the acquisition sequence.

[0458] Figure 10 illustrates two power Doppler images. The contrast gain provided by the correction is visible. [A] is a PWD display of the 400 frames of the application block before and [B] after correction.

[0459]

[0460] Directional Doppler Image

[0461] According to one embodiment of the ultrasonic construction process of this disclosure, the process further comprises a step of:

[0462] - determination of a 1D (x,z) directional Doppler image of an image point spatial position ultrasound (x, z), from the Fourier transform of the complex confocal signal Sc ( x, z, # ) according to the temporal dimension # m of the frames:

[0463] Xc(x, z, v) Sc(x, z, #

[0464] Where the frequency v is the variable conjugate to the acquisition time #.

[0465] The initial power Doppler image can thus be separated into a positive component and a negative component of the displacement:

[0466] f, . , ,2 I+(x,z) =J0 av\Xc(x,z, v) | LUH-o / JJ IV / XI 2 l (X, Z) = ixdv I Xc (X, Z, V) |

[0468]

[0469] These two images can then be combined to form a coloured Doppler image ID (x,z) in which blue encodes the negative component and road the positive component of tissue movement at each point (x,z).

[0470]

[0471] Map of local dynamics

[0472] According to one embodiment of the ultrasonic construction process of this disclosure, the process further comprises a step of:

[0473] - determination of a local dynamics map of an ultrasound image point spatial position (x, z), from the average frequency VD (x,z) of the Fourier transform Xf(x, z, v) of the complex confocal signal: [°4741 Vn( X, Z ) — "735---------

[0475] This average frequency is related to the average velocity of the scatterers at each point in the medium. Figure 11 shows these local dynamic maps, allowing a comparison before [A] and after [B] correction of transverse aberrations. Aberration correction allows for a more precise estimation of the scatterer movement velocity in the medium.

[0476]

[0477] Static confocal image

[0478] Beyond the dynamic component of the acquired data, the focusing laws obtained can also be used to correct the raw reflection matrix, i.e. the reflection matrix considered before filtering the static component of the data.

[0479]

[0480] The preferred field of application is ultrasound imaging as a whole, and in particular its Doppler and ULM (Ultrasound Localization Microscopy) modes. The invention can be used to image the human brain through a skull as well as the vascular network of the liver for medical imaging. It can be used to image dynamic processes using non-destructive testing. The invention also has potential applications in various fields of wave physics, such as optical coherence tomography for light or seismology for seismic waves.

[0481]

[0482] Of course, the invention is not limited to the examples just described. Many modifications can be made to these examples without departing from the scope of the present invention as described.

Claims

Demands

1. Ultrasonic construction method of a confocal image of a dynamic medium, the method comprising the following steps: a) acquisition, by means of a network of transducers, of a series of canonical reflection matrices R ui(t, #m)=[R(u out,im )] at different times, each canonical reflection matrix is ​​defined between an ultrasonic wave emission basis i at the input and a reception basis u at the output; the coefficients of this canonical reflection matrix correspond to the signals received by the transducers and induced by the ultrasonic waves reflected in the medium; t denoting the echo time and #m denoting the mth canonical reflection matrix;b) determination of a focused reflection matrix R xx(z,#m) for each canonical reflection matrix by input-output focusing for any point of at least one region of the medium, the coefficients of this focused reflection matrix are obtained by calculating an acoustic pressure field between all points of the region with lateral positions x in and x out, located at an expected depth z for a speed of sound assumed to be c0;c) determination of a dynamic component for each coefficient of each focused reflection matrix so as to constitute dynamic Ritz focused reflection matrices, # „) d) determination of a correction law ^x, z) for each point x and depth z of the medium from the dynamic focused reflection matrices, e) determination of corrected focused reflection matrices Pæ" application of the correction law 0(x, z) at every point of the medium, f) determination of a dynamic confocal signal S c (x,z, # m) of any point of spatial position (x, z) from the diagonal coefficients of the corrected focused reflection matrix # m), g) construction of an image from the dynamic confocal signals.;

2. A method according to claim 1, characterized in that the step of acquiring a series of canonical reflection matrices R ui(t, #m) includes the emission of an ultrasonic pulse from each transducer of the network whose position is located by the coordinate u in; this pulse gives rise to a divergent cylindrical or spherical incident wave which is reflected by diffusers of the medium; these reflected echoes form a backscattered field which is recorded by each of the transducers as a function of time; the canonical reflection matrix Ruu(t,#m) expressed in the basis of the transducers being composed of a set of impulse responses R(u out,um ) between transducers.

3. A method according to claim 1, characterized in that the step of acquiring a series of canonical reflection matrices R ui(t, #m) comprises an insonification of the medium with a series of plane waves with a delay r' applied to each signal at emission for the formation of a wavefront inclined at an angle 0 in with respect to the array of transducers, a backscattered field by the medium, R(u out, 0 in, t, #m) is measured by all the position transducers u out for each incident plane wave 0 in, the set of responses forming a canonical reflection matrix Ru 0(t, # ,„)=[ R(u out, 0 in, t, # ,„)].

4. A method according to claim 1, characterized in that the step of acquiring a series of canonical reflection matrices R ui(t, #m) comprises an insonification of the medium with a series of divergent waves.

5. A method according to any one of the preceding claims, characterized in that the step of determining the focused reflection matrix R xx(z,# m) comprises: - an input focusing process from each canonical reflection matrix R ui(t,# m) which uses a forward time of flight of the waves between the ultrasonic wave emission base i and a virtual input transducer TVin and which creates a focal spot called the input spot around a first point PI with spatial position r in = (x in ,z), said input focal spot corresponding to the virtual input transducer TVin, and - an output focusing process from the canonical reflection matrix R ui(t,# m) which uses a return time of flight of the waves between a virtual output transducer TVout and the receiving base transducers u and which creates a focal spot called the output spot around a second point P2 with spatial position r out =(x out,z), said output focal spot corresponding to the virtual output transducer TVout.

6. Method according to claim 5, characterized in that xout is located at a distance of xin less than or equal to a maximum distance Axmax which is a function of a number of ultrasonic waves generated during the acquisition of the reflection matrix.

7. A method according to any one of the preceding claims, characterized in that the focused reflection matrix is ​​determined in the time domain or in the frequency domain.

8. A method according to any one of the preceding claims, characterized in that the dynamic component is determined by subtracting from each coefficient of each focused reflection matrix a moving average over N focused reflection matrices.

9. A method according to any one of claims 1 to 7, characterized in that the dynamic component is determined by applying a high-pass or band-pass filter along a dimension of the numbers #m of the reflection matrices.

10. A method according to any one of claims 1 to 7, characterized in that the dynamic component is determined by performing a singular value decomposition of the focused reflection matrices rearranged in a two-dimensional matrix, one of whose dimensions is that of the numbers #m of the reflection matrices.

11. A method according to any one of the preceding claims, characterized in that the step of determining a correction law comprises, for each dynamic focused reflection matrix > the following steps: - determination of a dual reflection matrix Rcx(z, # m) by direct or indirect forward projection of the dynamic focused reflection matrix R^j- # onto a correction basis (c), - calculation of the correction law O(x, l) from the dual reflection matrix Rcx(z, # m), said correction law being a correction law, z) = [$,c, x)] onto the correction basis (c), - determination of the corrected focused reflection matrices P31 projection back of the corrected dual reflection matrices R^x(z, # m) to the focused basis (x).

12. A method according to claim 11, characterized in that the coefficients R^x(^ # = [r'c(x,c, z, # m) ] of the corrected dual reflection matrix are determined by performing a term-by-term product between the dual reflection matrix Rcx(z, # m) and the phase conjugate of the correction law Q(x, z), i.e.: Rcx = Rcx° O* where the symbol * denotes a phase conjugation operation, the symbol O is the Hadamard product, such that: Rcx(x,c,z, # = Rcx(x, c, z, # m)&"(x,c,z)-

13. A method according to claim 11 or 12, characterized in that the step of calculating the correction law 0(x, z) comprises the following steps: - construction of a correlation matrix C(x,z) from the dual reflection matrices Rcx(z, z, z) for each point (x,z) of a field of view, - determination of the focusing law ^(x^) for each point of the field of view by performing one of the following operations: - eigenvalue decomposition of the correlation matrix C(x,z), the correction law ^(x,z) being the first eigenvector U of the correlation matrix C(x,z) in the correction basis (c), - singular value decomposition of the rearranged dual reflection matrix as follows: Rc#(x, z) = (C, z) [x, z] the correction law <j>(x, z) being equal to the first singular vector of the dual reflection matrix Rc#, i.e., z) = Uj - solving the following equation: 0(x,z) = exp(j argîC^x, z)xO(x, z)} ) iteratively by the following expression, which corresponds to an iterative phase-reversal calculation: 0K+1(x, z) = exp( j arg{Cccx 0„(x, z)} ) Where x is the matrix product, with an arbitrary wavefront, the correction law (p(x, z)) being obtained by: ¢( x, z) = x, z). - solving the following equation: W(x, z) -exp(j arg{C##(x,z) x W(x, z)} ) where x is the matrix product, iteratively by the following expression: W„+1(x, z) = exp(j argfC^Cx, z) x W„(x, z)} ) where x is the matrix product, with Wo an arbitrary wavefront, which allows us to obtain the following vector W ( x, z ) : W(x, z) = limW„(x, z) the correction law / Jetant obtained by: 0 ( c, x, z ) = exp {jx arg {( x, c, z, # m ) W\x, c, z, # m)}} .

14. A method according to claim 13, characterized in that the correlation matrix C(x,z) is determined in the correction basis c and in the frequency domain by the following calculation of the elements of the correlation matrix C = C Cc : C(c, c', X, z) R(x,c, z, #R"(x,c',z, #m) where * is the conjugation operator. c and c' being points of the correction basis c

15. A method according to claim 13, characterized in that the correlation matrix C(x,z) is determined in the number basis #m by the following calculation of the elements of the correction matrix C = C##; C(#m, #bX,z) = HCR(X, C, z, #m) XR^X, C, Z, # / ) * is the conjugation operator. c and c' are points of the correction basis c, #m and # / denoting the same and same frames of the sequence of recorded reflection matrices

16. A method according to any one of claims 11 to 15, characterized in that the step of calculating the correction law z) is iterated at least twice with, at each iteration, the forward projection uses the corrected focused reflection matrix R^ ( z, # m ) obtained during the back projection of the previous iteration instead of the focused reflection matrix Rxx(z, # m).

17. A method according to any one of claims 11 to 16, characterized in that the step of determining a matrix of dual reflection Rcx(z, # m) is iterated at least twice with, at each iteration, the use of a different correction basis c.

18. A method according to any one of claims 11 to 16, characterized in that the step of determining a dual reflection matrix Rcx(z, # m) is iterated at least twice with, at each iteration, the use of a forward projection either to an input correction basis or to an output correction basis of the dynamic focused reflection matrix.

19. A method according to any one of claims 11 to 18, characterized in that the correction basis is one of the following: - a plane wave basis or spatial Fourier basis, - a transducer u basis, - a basis corresponding to the assumed location of aberrators in the medium, - a basis corresponding to a plane determined by optimization.

20. A method according to any one of claims 11 to 19, characterized in that the step of determining a dual reflection matrix Rcx(z, #,«) is carried out by forward projection of the dynamic focused reflection matrix R^Jjj # to the correction basis (c) considering a model propagator describing the propagation of waves from the focused basis (x) to the correction basis (c).

21. 21. A method according to any one of the preceding claims, characterized in that the step of determining corrected focused reflection matrices is iterated at least twice, with at each iteration, the forward projection using the corrected focused reflection matrix R* (# ) obtained during the back projection of the previous iteration instead of the focused reflection matrix RXXU # m)-

22. 22. A method according to any one of the preceding claims, characterized in that step g) of constructing an image from the dynamic confocal signals comprises a step of: - determining an intensity of each dynamic confocal signal to construct a confocal image of the medium, or - determining an intensity of each dynamic confocal signal to construct a power Doppler image by summing, for each point of spatial position (x, z), the intensities of several confocal signals, or - determination of a Fourier transform of the dynamic confocal signal Sc ( x, z, # ) according to the temporal dimension # m : Xc(x, z, v) =E* Sc(x, z, # where the frequency v is the variable conjugate to the acquisition time # to construct a directional Doppler image of the medium, - determination of the local dynamics for each point of ultrasound image of spatial position (x, z), by measuring the average frequency VD (x,z) of the Fourier transform Xc(x, z, v) of the complex confocal signal: Z) — rW to construct a map of the axial velocity of the scatterers at each point of the image.

23. A method according to any one of the preceding claims, characterized in that the focused reflection matrix R xx(z,#m) depends on a frequency f of the ultrasonic signals, the correction law being determined as a function of this frequency f and the coordinates of the correction plane.

24. 24. A method according to any one of the preceding claims, characterized in that the focused reflection matrix R xx(z,#m) is integrated over the entire bandwidth, the correction law being determined solely as a function of the coordinates of the correction plane (c).

25. 25. A method according to claim 13 and any one of claims 14 to 23, characterized in that the correlation matrix for obtaining the correction law, ¢(x, z) = [$(f, x, z) is given by: C(J,f\x,z) =E# S(x,z,f, S\x,z,f, #m) f and f' being frequencies in the bandwidth of the ultrasonic signal, the correction basis being directly named here (f), i.e. (c) = (f) S (x, zf, #m) being the confocal signal of the dynamic focused reflection matrix f- # Y

26. 26. A method according to claim 11 or 12, characterized in that the step of calculating the correction law ([Xx, z]) can be carried out: - by defining a correlation matrix for obtaining the correction law, 0 ( x, z ) = [$ ( f, x, z ) ], by: C(f, f',x,z) #m) S*(x,z,f', #m) - and by averaging the correlation matrix over the numbers #m but also over adjacent pixels belonging to the same isoplanetism patch Qp: C(ff xp, zp ) = E#( x, z, f, # m ) ( x, z, f, # m ) f and f' denote a frequency of the confocal signal in the bandwidth of the probe xp is the transverse coordinate of the central point of the isoplanetism patch for which we seek to estimate the aberration law zp is the axial coordinate of the central point of the isoplanetism patch for which we seek to estimate the aberration law ® ^p-. being an isoplanetism patch centered on the point (Xp. Zp).

27. ​​Ultrasonic construction system for a confocal image of a dynamic medium, the system comprising: - an array (10) of transducers adapted to generate a series of ultrasonic waves incident in an area of ​​interest of the medium, and to measure as a function of time the ultrasonic waves backscattered by said area of ​​interest; and - a computing unit (30) connected to the array of transducers and adapted to implement the method according to any one of claims 1 to 26.

28. Product computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the steps of the process according to any one of claims 1 to 26.

29. Computer-readable medium comprising instructions which, when executed by a computer, cause the computer to carry out the steps of the process according to any one of claims 1 to 26.< / j>

Images