Method and system for non-invasively characterizing a heterogeneous medium using ultrasound

12257103 ยท 2025-03-25

Assignee

Inventors

Cpc classification

International classification

Abstract

Method for non-invasively characterizing a heterogeneous medium using ultrasound, comprising a step of generating a series of incident ultrasonic waves, a step of recording an experimental reflection matrix R.sub.ui(t) defined between the input emission basis (i) and an output reception basis (u), a step of determining a response REP(r,r) of the medium between an input virtual transducer (TV.sub.in) calculated based on an input focusing of the experimental reflection matrix that creates an input focal spot around a first point (P1), and an output virtual transducer (TV.sub.out) calculated based on an output focusing of the experimental reflection matrix that creates an output focal spot around a second point (P2), said response being expressed as a function of a central point (PC) of spatial position (r) in the medium located midway between the first and second points (P1, P2).

Claims

1. Method for non-invasive ultrasound characterization of a heterogeneous medium, the method comprising: generating a series of incident ultrasonic waves in a region of the heterogeneous medium, generated by a transducer array, said series of incident ultrasonic waves constituting an input emission basis (i); recording an experimental reflection matrix R.sub.ui(t) defined between the input emission basis (i) and an output reception basis (u); determining a response REP(r,r) of the medium between an input virtual transducer (TV.sub.in) of spatial position r.sub.in calculated based on an input focusing of the experimental reflection matrix corresponding to an input focal spot around a first point (P1) and an output virtual transducer (TV.sub.out) of spatial position r.sub.out calculated based on an output focusing of the experimental reflection matrix corresponding to an output focal spot around a second point (P2), said response REP(r,r) being expressed as a function of a central point (PC) of spatial position r in the medium and as a function of a distance coordinate r from said central point, said central point (PC) being located midway between the first and second points (P1, P2) such that r=(r.sub.out+r.sub.in)/2 and such that r=(r.sub.outr.sub.in)/2, and said central point being the origin of a measurement axis (AXm) passing through said first and second points (P1, P2), said measurement axis (AXm) forming an angle relative to a first axis (X) of the medium, and the first point (P1) being at the distance coordinate +r on the measurement axis (AXm), and the second point (P2) being at an opposite distance coordinate r on the measurement axis (AXm).

2. Method according to claim 1, wherein, in determining the response REP(r,r), the input focusing of the experimental reflection matrix uses an outward time-of-flight of the waves between the input emission basis and the input virtual transducer, and the output focusing uses a return time-of-flight of the waves between the output virtual transducer and the output reception basis.

3. Method according to claim 1, wherein the response REP(r,r) of the medium is calculated by the following formula: REP ( r , r ) = 1 N i n N o u t .Math. i i n .Math. u o u t R u i ( u out , i in , ( r in , r out , u out , i i n ) ) in which N.sub.in is the number of elements of the input emission basis, N.sub.out is the number of elements of the output reception basis, R.sub.ui(t) is the experimental reflection matrix, in which R.sub.ui (u.sub.out, i.sub.in, (r.sub.in, r.sub.out, u.sub.out, i.sub.in)) is the element of the experimental reflection matrix R.sub.ui(t) recorded by transducer u.sub.out consecutive to emission i.sub.in at time , is a time which is the sum of the outward time-of-flight .sub.in of the ultrasonic wave between the input emission basis (i) and the first point (P1) and the return time-of-flight Tout of the ultrasonic wave between the second point (P2) and the output reception basis (u), as explained by the following formula:
(r.sub.in,r.sub.out,u.sub.out,i.sub.in)=.sub.in(r.sub.in,i.sub.in)+.sub.out(r.sub.out,u.sub.out).

4. Method according to claim 1, further comprising: determining a response profile PR(r) which is a plurality of responses REP(r,r) calculated for a plurality of values of the distance coordinate r and for the central point (PC), and for the measurement axis (AX.sub.m), corresponding to the angle , r being the distance of the second point from the central point, wherein the value of said distance being such that r=u.sub..Math.u.sub., u.sub. being a unit vector in the direction of the measurement axis AX.sub.m defined by the angle .

5. Method according to claim 4, further comprising: determining the resolution w(r) of the central point (PC) based on a modulus of the response profile PR(r), in which the resolution w(r) is a width of a peak of said modulus of the response profile PR(r), centered around a distance being null (r=0).

6. Method according to claim 5, wherein the width of the peak is estimated at a height which is a portion of a maximum height of said peak, said portion being half the maximum height of the peak.

7. Method according to claim 5, further comprising: determining a focusing criterion F(r) of the central point (PC) based on a resolution w(r) and a theoretical resolution w.sub.0(r), said theoretical resolution w.sub.0(r) being determined based on the input emission basis (i) and the output reception basis (u).

8. Method according to claim 7, further comprising: calculating the focusing criterion F(r) at a plurality of points in the medium, and at least one calculation parameter for the input focusing and/or output focusing is optimized by minimizing or maximizing a mean of said focusing criterion F(r) for said plurality of points.

9. Method according to claim 8, wherein said at least one calculation parameter comprises the speed of sound in the medium.

10. Method according to claim 7, wherein the theoretical resolution is determined by a technique comprised in the list comprised in: a first analytical calculation at the central point (PC), for a pulse (.sub.1) and the input emission basis (i) and the output reception basis (u), in which the theoretical resolution is calculated by the angle from which the transducer array is viewed from the central point (PC), a second analytical calculation at the central point (PC), for a pulse range () and the input emission basis (i) and the output reception basis (u), in which the theoretical resolution is an integral calculation over said pulse range of the angle from which the transducer array is viewed from the central point (PC) weighted by the frequency spectrum of signals from the experimental reflection matrix R.sub.ui(t), and a third calculation of wave propagation simulation, firstly between the first point of the medium corresponding to an input virtual transducer (TV.sub.in) and the input emission basis (i), and secondly between the second point of the medium corresponding to an output virtual transducer (TV.sub.out) and the output reception basis (u), said simulation using the response REP(r,r) and a model of wave propagation in the medium.

11. Method according to claim 1, further comprising: determining a level of symmetry (r) which is the mean correlation coefficient between two reciprocal responses, said mean being calculated for distance coordinate values of modulus greater than a predetermined resolution w.sub.d(r) and/or being calculated for a range of angle values or for a predetermined angle value .sub.d, determining a first response REP1(r,r) of the medium between an input virtual transducer (TV.sub.in) of spatial position r.sub.in calculated based on an input focusing of the experimental reflection matrix corresponding to an input focal spot around a first point (P1) and an output virtual transducer (TV.sub.out) of spatial position r.sub.out calculated based on an output focusing of the experimental reflection matrix corresponding to an output focal spot around a second point (P2), and determining a second response REP2(r,r) of the medium between an input virtual transducer (TV.sub.in) of spatial position r.sub.in calculated based on an input focusing of the experimental reflection matrix corresponding to an input focal spot around the second point (P2) and an output virtual transducer (TV.sub.out) of spatial position r.sub.out calculated based on an output focusing of the experimental reflection matrix corresponding to an output focal spot around the first point (P1).

12. Method according to claim 11, wherein the level of symmetry is calculated by the following formula: ( r ) = .Math. R E P 1 ( r , r ) .Math. REP 2 ( r , - r ) .Math. ( | r |> w d ( r , ) , ) .Math. .Math. "\[LeftBracketingBar]" REP 1 ( r , r ) .Math. "\[RightBracketingBar]" 2 .Math. .Math. .Math. "\[LeftBracketingBar]" REP 2 ( r , - r ) .Math. "\[RightBracketingBar]" 2 .Math. or by the following formula: ( r ) = .Math. Re [ REP 1 ( r , r ) .Math. REP 2 ( r , - r ) ] .Math. ( | r | > w d ( r , ) , ) .Math. .Math. "\[LeftBracketingBar]" REP 1 ( r , r ) .Math. "\[RightBracketingBar]" 2 .Math. .Math. .Math. "\[LeftBracketingBar]" REP 2 ( r , - r ) .Math. "\[RightBracketingBar]" 2 .Math. in which Re is a real part mathematical operator, |.Math.| is a modulus mathematical operator, <.Math.> is a mean mathematical operator, and * is a complex conjugation operator.

13. Method according to claim 11, further comprising: determining a first multiple scattering indicator (r) which is calculated by: ( r ) = ( r ) 1 - ( r ) (r) being the level of symmetry defined at the central point (PC) of the medium.

14. Method according to claim 11, further comprising: determining an afocal intensity I.sub.off(r) which is the mean of a squared modulus of the responses REP(r,r), the mean being calculated for distance values of modulus greater than a predetermined resolution w.sub.d(r) and/or being calculated for a range of angle values or for the predetermined angle value .sub.d, for example calculated by the following formula:
I.sub.off(r)=custom characterI.sub.off(r,r)custom character.sub.(|r|>w.sub.d.sub.(r)) determining a multiple scattering intensity I.sub.M(r) which is the product of the level of symmetry (r) and the afocal intensity I.sub.off(r), wherein:
I.sub.M(r)=(r).Math.I.sub.off(r) determining a noise intensity I.sub.N(r) being the product of one minus the level of symmetry (r) and the afocal intensity I.sub.off(r), wherein:
I.sub.N(r)=(1(r)).Math.I.sub.off(r).

15. Method according to claim 14, wherein the first multiple scattering indicator (r) is calculated by a ratio between the multiple scattering intensity I.sub.M(r) and the noise intensity I.sub.N(r), wherein: ( r ) = I M ( r ) I N ( r ) .

16. Method according to claim 14, further comprising: determining a confocal intensity I.sub.on(r) which is the value of a squared modulus of the response REP(r, |r|=0) for a distance coordinate value of zero (|r|=0), wherein for a point of the medium for which the first point (P1), the second point (P2), and the central point (PC) are coincident, determining a single scattering intensity I.sub.S(r) calculated based on the following equation:
I.sub.on(r)=I.sub.S(r)+2I.sub.M(r)+I.sub.N(r).

17. Method according to claim 14, further comprising: determining a second multiple scattering indicator (r) which is calculated by: ( r ) = I M ( r ) I S ( r ) .

18. Method according to claim 1, further comprising: determining an image of a local characterization parameter of the medium, said local characterization parameter being determined based on the response REP(r,r).

19. System for non-invasive ultrasonic characterization of a heterogeneous medium, the system comprising: a first array of transducers suitable for generating a series of incident ultrasonic waves in a region of the heterogeneous medium, said series of incident ultrasonic waves constituting an input emission basis (i), and suitable for recording as a function of time the ultrasonic waves backscattered by said region; and a computing unit coupled to the first array of transducers and configured for: recording an experimental reflection matrix R.sub.ui(t) defined between the input emission basis (i) and an output reception basis (u); determining a response REP(r,r) of the medium between an input virtual transducer (TV.sub.in) of spatial position r.sub.in calculated based on an input focusing of the experimental reflection matrix corresponding to an input focal spot around a first point (P1) and an output virtual transducer (TV.sub.out) of spatial position r.sub.out calculated based on an output focusing of the experimental reflection matrix corresponding to an output focal spot around a second point (P2), said response REP(r,r) being expressed as a function of a central point (PC) of spatial position r in the medium and as a function of a distance coordinate r from said central point, said central point (PC) being located midway between the first and second points (P1, P2) such that r=(r.sub.out+r.sub.in)/2 and such that r=(r.sub.outr.sub.in)/2, and said central point (PC) being the origin of a measurement axis (AX.sub.m) passing through said first and second point s (P1, P2), said measurement axis (AX.sub.m) forming an angle relative to a first axis (X) of the medium, and the first point (P1) being at the distance coordinate +r on the measurement axis (AXm), and the second point (P2) being at an opposite distance coordinate r on the measurement axis (AXm).

Description

BRIEF DESCRIPTION OF FIGURES

(1) Other features and advantages of the technique presented above will be apparent from reading the detailed description below, presented in a non-limiting manner for illustrative purposes, made with reference to the figures in which:

(2) FIGS. 1A to 1C (already described) illustrate known emission/reception mechanisms for ultrasound imaging and quantification;

(3) FIG. 2 (already described) illustrates the impact of aberrations in ultrasound imaging according to the prior art;

(4) FIG. 3 illustrates an example of a system for ultrasonic characterization, for implementing the methods for ultrasonic characterization according to this description;

(5) FIG. 4 illustrates the definitions used in the method for ultrasonic characterization according to this description;

(6) FIG. 5 illustrates an example image representing the modulus of a response matrix REP(r,r) according to the method for ultrasonic characterization as in FIG. 4;

(7) FIG. 6 illustrates an example response profile PR(r) corresponding to the response matrix of FIG. 5;

(8) FIG. 7 illustrates ultrasound images of three heterogeneous media, part A of this figure corresponding to a test medium referred to as phantom, part B of this figure corresponding to a medium having a layer of meat placed on the test medium, and part C of this figure corresponding to a medium which is an in vivo liver;

(9) FIG. 8 illustrates an image of resolution w(r) (on the left) and an image of theoretical resolution w.sub.0(r) (on the right), established for the test medium (of part A) of FIG. 7;

(10) FIG. 9 illustrates images of the focusing criterion F(r), established for the three heterogeneous media of FIG. 7 (parts A, B, and C respectively);

(11) FIG. 10 illustrates images of the first multiple scattering indicator (r) for focusing, established for the three heterogeneous media of FIG. 7 (parts A, B, and C respectively);

(12) FIG. 11 illustrates images of the second multiple scattering indicator (r) established for the three heterogeneous media of FIG. 7 (parts A, B, and C respectively);

(13) FIG. 12 illustrates three calculations of optimal speeds of sound, carried out based on three models of medium B of FIG. 7 and based on the defined focusing criterion F(r).

(14) In the various embodiments which will be described with reference to the figures, similar or identical elements bear the same references.

DETAILED DESCRIPTION

(15) In the following detailed description, only certain embodiments are described in detail to ensure clarity of the description, but these examples are not intended to limit the general scope of the principles that emerge from this description.

(16) The various embodiments and aspects described in this description can be combined or simplified in multiple ways. In particular, the steps of the various methods can be repeated, inverted, and/or executed in parallel, unless otherwise specified.

(17) Certain elements depicted in the specification and claims are denoted in bold, meaning these elements are vectors. A person of ordinary skill in the art would recognize the context in which the following elements are used, and would understand that these elements represent vectors, even if said elements are not depicted in bold in every instance in the publication of the present application, or in any patent that may issue therefrom. These elements include, but are not necessarily limited to: R.sub.ui (e.g., as used in the context of experimental reflection matrix R.sub.ui(t)) i (e.g., as used in the context of input emission basis (i), as well as in i.sub.in); u (e.g., as used in the context of output reception basis (u), as well as in u.sub.out); r (e.g., as used in the context of spatial positions r, r.sub.in, and r.sub.out); r (e.g., as used in the context of distance coordinates +r and r); and u.sub. (e.g., as used in the context of r.Math.u.sub.).

(18) FIG. 3 illustrates an example of a system 40 for ultrasonic characterization, for implementing methods for ultrasonic characterization of a heterogeneous medium 20, according to this description. The system 40 comprises at least a first array 10 of transducers 11, for example a linear or two-dimensional array; the transducers are for example piezoelectric ultrasonic transducers which may be in the conventional form of a rigid bar placed in contact with the medium 20. The array of transducers is, for example, part of a probing device 41, hereinafter also referred to by the more common term probe; the array of transducers is connected to a computing unit 42, which itself may be connected to a display device 43; the computing unit transmits and records electrical signals to/from each of the transducers 11. The ultrasonic transducers then convert these electrical signals into ultrasonic waves and vice versa. Connection or link between the probing device 41, the computing unit 42, and the display device 43, is understood to mean any type of wired connection of the electrical or optical type, or any type of wireless connection using any protocol such as WiFi, Bluetooth, or others. These connections or links are one-way or two-way.

(19) The computing unit 42 is configured to implement calculation or processing steps, in particular to implement the steps of methods according to this description. By convention, a spatial reference system of the medium 20 is defined by taking a first axis X and a second axis Z perpendicular to the first. For simplicity, the first axis X corresponds to the direction in which the transducers 11 are aligned for a linear array, and the second axis Z corresponds to the depth of the medium 20 relative to this array 10 of transducers 11. This definition can be extended to a three-axis spatial reference system in the case of a two-dimensional array.

(20) In FIG. 3, as in the rest of the description, reference is made to an array of transducers for emission and reception, it being understood that, in a more general case, several arrays of transducers could be used simultaneously. They can be both transmitter and receiver, or only a transmitter for some and only a receiver for others.

(21) When reference is made in this description to calculation or processing steps, in particular for implementing steps of the methods, it is understood that each calculation or processing step can be implemented by software, hardware, firmware, microcode, or any suitable combination of these technologies. When software is used, each calculation or processing step may be implemented by computer program instructions or software code. These instructions may be stored in or transmitted to a storage medium readable by a computer (or computing unit) and/or executed by a computer (or computing unit) in order to implement these calculation or processing steps.

(22) Defining the Analysis of a Point of the Medium

(23) This description describes methods and systems for non-invasive ultrasonic characterization of a heterogeneous sample. These methods and systems are based on definitions shown in FIG. 4: a central point PC of spatial position r is defined in the spatial reference system of the medium, and located midway between a first point P1 and a second point P2. A measurement axis AX.sub.m is defined, passing through the first point P1 and the second point P2, and forming an angle with the first axis X of the array of transducers 11. The central point PC is located on the origin of the measurement axis AX.sub.m (distance coordinate of zero on the measurement axis). The first point P1 is at a distance coordinate r and the second point P2 is at a distance coordinate +r from the central point PC, origin of the measurement axis.

(24) The spatial position r and the distance coordinate r are denoted in bold, meaning that these elements are vectors of the position and offset relative to a position, vectors in the spatial reference system of the medium (X,Z). The distance coordinate vector r thus takes into account the direction of the measurement axis AX.sub.m, and its angle relative to the first axis X. Other definitions of the positions of points relative to other points are possible and accessible to any specialist in the field of ultrasound. In particular, the first and second points can be identified by a distance |r| and the angle , or by another position identifier.

(25) These two points P1 and P2 can be at a short distance from each other, i.e. a few millimeters from each other, and for example 20 millimeters or less.

(26) As shown in FIG. 4, the method for non-invasive ultrasound characterization comprises: a step of generating a series of incident ultrasonic waves US.sub.in in a region of said heterogeneous medium, by means of an array 10 of transducers 11, said series of incident ultrasonic waves being an emission basis i; and a step of recording an experimental reflection matrix R.sub.ui(t) defined between the input emission basis i and an output reception basis u; a step of determining a response REP(r,r) of the medium between an input virtual transducer TV.sub.in of spatial position r.sub.in at the first point P1 and an output virtual transducer TV.sub.out of spatial position r.sub.out at the second point P2, said response being expressed as a function of the central point PC of spatial position r.

(27) As the central point PC is midway between the two points P1, P2, we have the following relations:
r=(r.sub.out+r.sub.in)/2 and r=(r.sub.outr.sub.in)/2

(28) The input emission basis i being for example a basis of waves each generated by one of the transducers 11 of the array 10 or a basis of plane waves of angular inclination with respect to the axis X, as described above in the description of FIGS. 1A to 1C.

(29) The reception basis u is usually the basis of the transducers 11.

(30) Thus, the step of generating ultrasonic waves is understood to be between the emission basis i and the reception basis u. This ultrasonic generation step is therefore defined for any type of ultrasonic wave of the focused or unfocused type, such as plane waves.

(31) In the recording step, the experimental reflection matrix R.sub.ui(t) is therefore defined between the input emission basis i and an output reception basis u. This matrix contains all the temporal responses of the medium, measured at time t by each transducer 11 of spatial coordinate u.sub.out for each emission i.sub.in. It is understood that the elements named with the index in refer to emission (i.e. the input) and the elements named with the index out refer to reception (i.e. the output).

(32) In the step of determining a response REP(r,r), we apply: an input focusing process based on the experimental reflection matrix R.sub.ui(t) that creates an input focal spot around the first point P1, said input focal spot corresponding to the input virtual transducer TV.sub.in, and an output focusing process based on the experimental reflection matrix R.sub.ui(t) that creates an output focal spot around the second point P2, said output focal spot corresponding to the output virtual transducer TV.sub.out.

(33) The input focusing process uses an outward time-of-flight of the waves between the emission basis (i) and the input virtual transducer TV.sub.in. The output focusing process uses a return time-of-flight of the waves between the output virtual transducer TV.sub.out and the transducers of the reception basis (u). These input and output focusing processes in fact form an input-output focusing process, hereinafter referred to as the focusing process.

(34) The first point P1 being relative to the input virtual transducer TV.sub.in, it is therefore located at a coordinate r on the measurement axis AX.sub.m in relation to the central point PC, and the second point P2 being relative to the output virtual transducer TV.sub.out, it is therefore located at a coordinate +r on the measurement axis AX.sub.m in relation to the central point PC.

(35) The input focusing (at emission) is configured to concentrate the acoustic wave around point P1 over a spatial extent corresponding to the input focal spot. The scatterers located inside this region in the medium then generate a wave which is backscattered towards the probe. This region, characterized by the focal spot at emission and the reflectivity of the corresponding medium is called the input virtual transducer TV.sub.in and can be interpreted as a virtual source.

(36) The output focusing (at reception) is configured to select the echoes generated by scatterers located around point P2 over a spatial extent corresponding to the output focal spot. This region, characterized by the output focal spot (reception) and the reflectivity of the corresponding medium is called the output virtual transducer TV.sub.out and can be interpreted as a virtual sensor.

(37) The response REP(r,r) can therefore be interpreted as an estimate of the pressure field coming from position r.sub.out for a focusing at position r.sub.in.

(38) In other words, in this method of non-invasive ultrasonic characterization, the input virtual transducer TV.sub.in corresponds to an ultrasonic virtual source at spatial position r.sub.in in the medium and the output virtual transducer TV.sub.out corresponds to an ultrasonic virtual sensor at spatial position r.sub.out. Thus, the method is able to probe the medium around point P1 and/or point P2, spatially, which makes it possible to obtain new information about the propagation of the waves.

(39) For example, a calculation of the response REP(r,r) of the medium between the input virtual transducer TV.sub.in and the output virtual transducer TV.sub.out by a focusing process, which is for example an improved beamforming method, which can be expressed by the following simplified formula:

(40) REP ( r , r ) = 1 N i n N o u t .Math. i i n .Math. u o u t R u i ( u out , i i n , ( r i n , r o u t , u o u t , i i n ) ) ( Eq . 1 )
in which N.sub.in is the number of elements of the emission basis i, N.sub.out is the number of elements of the output reception basis u, R.sub.ui(t) is the experimental reflection matrix, in which R.sub.ui(u.sub.out, i.sub.in, (r.sub.in, r.sub.out, u.sub.out, i.sub.in)) is the element of the experimental reflection matrix R.sub.ui(t) recorded by transducer u.sub.out consecutive to emission i.sub.in at time .

(41) The time is the sum of the outward time-of-flight .sub.in of the ultrasonic wave between the transducers of the emission basis i and the first point P1 and the return time-of-flight .sub.out of the ultrasonic wave between the second point P2 and the transducers of the reception basis u, as explained by the following formula:
(r.sub.in,r.sub.out,u.sub.out,i.sub.in)=.sub.in(r.sub.in,i.sub.in)+.sub.out(r.sub.out,u.sub.out)(Eq. 2)

(42) The times-of-flight .sub.in and .sub.out are calculated based on a speed of sound model. The hypothesis consists of considering a homogeneous medium with a constant speed of sound c.sub.0. In this case, the times-of-flight are obtained directly, based on the distances between the probe and the virtual transducers.

(43) The number of elements of the emission basis N.sub.in is for example greater than or equal to two. The number of elements of the reception basis N.sub.out is for example greater than or equal to two.

(44) This improved beamforming formula is therefore a double summation of the temporal responses recorded in the experimental reflection matrix R.sub.ui: a first summation linked to the input focusing (emission) according to emission basis i at point P1 of position spatial r.sub.in, and a second summation related to the output focusing (reception) according to reception basis u at point P2 of spatial position r.sub.out. This calculation is therefore carried out for the spatial coordinates of the two points P1 and P2 (of spatial positions r.sub.in, r.sub.out). The result of this improved beamforming formula is therefore a time signal which is a pressure field for these two spatial coordinates (r.sub.in, r.sub.out).

(45) Note that the particular case of calculating a response REP(r, r=0) corresponds to the situation in which the points P1 and P2 are coincident at the same spatial position r.sub.in=r.sub.out=r. This configuration corresponds exactly to the case of the usual confocal ultrasound imaging in which each pixel of the image results from a process of input focusing (at emission) and output focusing (at reception) at a same point in the medium, and for the points of images of the medium. The time-of-flight then corresponds to the outward and return time-of-flight required for a wave to propagate from the probe to the single point of spatial position r, then from this point r to each transducer of the probe.

(46) By means of these arrangements, the method makes it possible to probe the medium very locally in any direction corresponding to the measurement axis AX.sub.m, in order to extract, via the input and output focusing, more local information about the medium at the central point PC of spatial position r, between the first point and second point of the heterogeneous medium 20. This local information is entirely contained within the values of the calculated response, the response REP(r,r) of the medium which can be used to characterize each point of the medium, for example in terms of resolution or in terms of multiple scattering. This local information is entirely contained within the values of the calculated temporal response which can be exploited to characterize each point of the medium.

(47) Indeed, it is customary to deduce, from this temporal response after beamforming, an estimate of the reflectivity of the medium by considering the absolute value of the confocal signals characterized by the superposition of points P1 and P2, i.e. the superposition of the focal spots of the input and output virtual transducers (i.e. r.sub.in=r.sub.out and |r|=0). This reflectivity then corresponds to the value of a pixel of a conventional ultrasound-type image.

(48) The responses REP(r,r) can thus be determined for any set of real distance values |r|, for example between two limits such as r.sub.max and +r.sub.max, these limits being determined so that the input and output virtual transducers remain inside the medium 20. (|r|, ) are then the polar coordinates of the distance coordinates vector r.

(49) In this previous convention of spatial description around the central point PC of spatial position r, the response REP(r, r) corresponds to inverting the spatial positions of the input and output virtual transducers.

(50) The set of responses REP(r,r) can then be recorded in a matrix of the same name. This matrix of responses is a focused reflection matrix, which records a pressure field calculated at any point in the medium with the defined hypotheses.

(51) One therefore obtains a response matrix REP(r,r) (4-dimensional in the case of a linear probe, including two for r and r) which records focused time signals.

(52) FIG. 5 shows an image corresponding to a sub-matrix of the response matrix REP, said sub-matrix corresponding to a set of several central points PC of spatial position r in which the Z axis coordinate is fixed, and the angle is zero. Thus, in this image the abscissa corresponds to a variation along the X axis of the position of the central point PC and the ordinate corresponds to the distance coordinate r relative to this central point. The values of the points of this image (response) not on the abscissa axis of this image have a low (but not zero) value. The values of the points of this image (response) on the abscissa axis have a value corresponding to the intensity of the ultrasonic image point at the central point PC.

(53) Such an image can be extracted from the response REP(r,r) for variations in distance coordinates r on a single measurement axis AXm or several measurement axes, i.e. for one or more angle values .

(54) A polar image representing the variation of the modulus of the response as a function of the separation distance |r| and of the angle can also be constructed, which gives a representation of the variation of the response around a central point PC, and therefore of the focal spot at this point.

(55) Calculating a Response Profile PR(r)

(56) Having obtained the responses of the medium determined according to the above method, a step of determining a response profile PR(r) can be carried out, the response profile being a plurality of responses REP(r,r) calculated for a plurality of values of the distance coordinate r. This response profile PR(r) is considered for a same central point PC of spatial position r and along a same measurement axis AX.sub.m, corresponding to a same direction of angle . The response profile PR(r) is therefore determined for a plurality of distances r, the distance r being the abscissa of the second point P2 with respect to the central point PC, i.e. the value such that r=r.Math.u.sub., u.sub. being a unit vector in the direction of the measurement axis AX.sub.m defined by the angle . In other words, the response profile PR(r) is a vertical slice of the image of FIG. 5, and this response profile is the curve of said image slice.

(57) The responses REP(r,r) can be complex values, particularly when using a focusing formulation in complex values, as is known in in-phase/quadrature beamforming (known as IQ beamforming). Consequently, the response profile PR(r) can also be represented by any modulus of these complex responses.

(58) FIG. 6 shows a schematic example of the squared modulus of a response profile PR(r) that can be determined for a central point PC of the medium and for a direction of angle : here for =0 because the image of FIG. 5 is also established for this angle value.

(59) However, the angle can take any value between zero (0) and pi (), and therefore a response profile curve PR.sub.c can be plotted or determined for multiple angle values .

(60) The set of response profiles PR(r) or PR(r, ) (if several angles are used, but in the following we will only keep the spatial position for the sake of simplifying the description) or PR(r, r, ) (if also using the spatial position of the central point) can be recorded in a matrix of the same name.

(61) The response profile PR(r) presents: an over-intensity (maximum of the curve) for a low value of r which corresponds to a single scattering in which the width is related to the focal spot, and an incoherent background present for all values of r, which is due to echoes arising from the multiple scattering phenomenon and noise.

(62) The sub-matrix represented in FIG. 5 is therefore the set of response profiles PR(r) for a set of central points PC of spatial position r in which the Z axis coordinate is fixed, and the angle is zero (or constant).

(63) This response profile PR(r) is a basic representation making it possible to determine new parameters for the local characterization of the medium and/or for the performance of the ultrasound imaging process (i.e. beamforming). We will illustrate the results of these characterization parameters by images of the heterogeneous medium, constructed with said characterization parameters.

(64) These characterization parameters for three different cases of heterogeneous media are described in more detail below. FIG. 7 shows ultrasound images of these three heterogeneous media: in part A of this FIG. 7 (on the left), one can see an image of a test medium (medium A) referred to as phantom, which includes two cylinders of differing rigidity, in part B of this FIG. 7 (in the center), one can see an image of a medium having a layer of meat placed on the above test medium (medium B), and in part C of FIG. 7 (on the right), one can see an image of a medium which is an in vivo liver (medium C).
Calculating the Resolution of Point w(r)

(65) A step of determining the resolution w(r) of the central point PC in the direction of the measurement axis AX.sub.m of angle can then be performed, based on a modulus of the response profile. This resolution is therefore a local estimate of the resolution of the ultrasound image.

(66) Note that the modulus of the response profile as shown in FIG. 6 comprises a peak or maximum around the zero distance coordinate r (|r|=0). This peak or over-intensity of the response profile is linked to the single scattering echoes and appears when the focal spots of the two virtual transducers TV.sub.in and TV.sub.out overlap. The spatial extent of this peak is therefore strongly correlated with the spatial dimensions of the input and output focal spots along the direction of angle and therefore with the local resolution of the ultrasound image.

(67) The resolution w(r) can then be determined for example by the width of this peak. The width of this peak is for example determined at a height which is a portion of the maximum height of this peak. For example, the portion of the height will be one-half or one-third () or two-thirds () or any other ratio of the maximum height. The maximum height of the peak is in fact the intensity of the ultrasound image at the central point PC if one considers only the squared modulus of the response profile, i.e. |PR(r=0)|.sup.2, as is the case in the example illustrated in FIG. 6.

(68) It is understood that the resolution depends on the central point PC considered, but also on the angle .

(69) Therefore, the proposed method makes it possible to obtain at each point: the axial resolution for an angle of /2, i.e. by the value w(r) taken for a response profile corresponding to an angle =/2, and the lateral resolution for a zero angle , i.e. by the value w(r) taken for a response profile corresponding to an angle =0.

(70) The method makes it possible to define, at any point in the medium, the extent of the focal spot and therefore the resolution of the ultrasound method in each of the directions of angle .

(71) The image in the left part of FIG. 8 shows an example calculation of the resolution w(r, ) at each point of the test medium (medium A). Note that the resolution degrades with depth and when moving towards the edge of the image.

(72) Calculating the Theoretical Resolution w.sub.0(r)

(73) According to a first variant, the theoretical resolution w.sub.0(r) is determined by a first analytical calculation at the central point (PC) for a pulse (.sub.1), the emission basis (i), and the reception basis (u): It is calculated by the angle from which the transducer array is viewed from the central point (PC). It depends on the maximum half angle of aperture used during emission to insonify the central point of spatial position r, or during reception to collect the echoes coming from this central point.

(74) According to a second variant, the theoretical resolution w.sub.0(r) is determined by a second analytical calculation at the central point (PC) for a pulse range () and the emission basis (i) and the reception basis (u). It is obtained by an integral calculation over said pulse range and over the angle from which the transducer array is viewed from the central point (PC) weighted by the frequency spectrum of recorded signals. The latter can be obtained by averaging the modulus of the Fourier transform of the elements of the experimental reflection matrix R.sub.ui(t).

(75) According to a third variant, the theoretical resolution w.sub.0(r) is determined by a third calculation of wave propagation simulation, firstly between the first point of the medium corresponding to an input virtual transducer (TV.sub.in) and the emission basis (i), and secondly between the second point of the medium corresponding to an output virtual transducer (TV.sub.out) and the reception basis (u), said simulation using the response REP(r,r) and a model of wave propagation in the medium. This third calculation reflects the double focusing step carried out in order to calculate the response profiles PR(r, , r). This third simulation calculation consists of generating a theoretical reflection matrix associated with a random medium in which the speed of sound corresponds exactly to the speed of sound model assumed for calculating the responses REP(r,r). This simulation then uses the same emission basis and the same reception basis as those used for the physical experiment. The set of operations carried out to determine the resolution w(r) are then repeated to calculate the theoretical resolution w.sub.0(r) based on the theoretical reflection matrix generated. All the diffraction phenomena are entirely taken into account and an estimate of the theoretical resolution of the medium without aberration is thus obtained. Note that the statistical properties of the medium such as the average reflectivity of a region, the spectrum of the backscattered echoes, can be deduced from the responses REP(r,r) in order to use a simulation which best models the experiment performed.

(76) The image in the right part of FIG. 8 shows an example calculation of the theoretical resolution w.sub.0(r) at each point of the test medium (medium A). Note that the resolution degrades with depth and also when moving towards the edge of the image. Note that this image is very similar to the image in the left part of this same figure. The calculation of the resolution carried out by the preceding method is therefore in agreement with the calculation of the theoretical resolution.

(77) Calculating a Focusing Criterion F(r)

(78) A step of determining a focusing criterion F(r) of the central point PC can then be performed, based on the resolution w(r) and a theoretical resolution w.sub.0(r). The theoretical resolution is for example determined based on an input emission basis i, an output reception basis u, and a modeling of the propagation of ultrasonic waves in the medium.

(79) Usually, the focusing criterion F(r) is a ratio of said resolution and theoretical resolution, or the reverse (a simple rule). In other words, we can obtain:
F(r)=w(r)/w.sub.0(r)(Eq. 3)
or
F(r)=w.sub.0(r)/w(r)(Eq. 4)

(80) FIG. 9 illustrates images of the focusing criterion F(r), established for the three heterogeneous media of FIG. 7 (media A, B, and C, which correspond between FIGS. 7 and 9).

(81) A value of one (1) for this focusing criterion corresponds to an identical resolution and theoretical resolution (light in the figure). A value close to zero (0) for this focusing criterion corresponds to divergent resolution values (dark in the figure), i.e. a degraded focusing.

(82) The image of medium A illustrates great homogeneity with a mean of this focusing criterion that is close to 0.97. This means that the ultrasound image is well formed and that the focusing assumptions are correct. The images of media B and C show notable degradations corresponding to heterogeneities located upstream of the ultrasonic wave propagation: layer of meat at the surface for medium B and adipose or muscle tissue for the liver of medium C. The image of medium C highlights very degraded regions (dark regions at the bottom left of this image) which means that the image produced in FIG. 7 (part C) also has the same problems at the same locations. Here, the image of the focusing criterion can indicate to the operator of the ultrasound system that the left part of their ultrasound image is of poor quality, especially in spatial resolution. This can help the user to interpret the ultrasound image in order to establish a diagnosis preferentially on correctly imaged regions or to modify the manner in which he or she generates said ultrasound image. For example, this can encourage the operator to modify the positioning of the probe (i.e. the array of transducers) so as to obtain a focusing criterion of good quality in the area(s) of interest for medical diagnosis.

(83) Advantageously, the imaging system which implements the present technique will be able to superimpose an ultrasound image on an image of the focusing criterion.

(84) Calculating a Level of Symmetry (r)

(85) A step of determining a level of symmetry (r) of the central point PC can then be performed, this level of symmetry being a mean correlation coefficient between two reciprocal responses: i.e. by exchanging positions r.sub.in and r.sub.out of the input virtual transducer TV.sub.in and output virtual transducer TV.sub.out. The criterion of acoustic reciprocity of the medium at the central point PC of spatial position r is thus tested. Note that the signals resulting from multiple scattering phenomena are reciprocal while those resulting from electronic noise are not correlated.

(86) For this method, the following are determined: a first response REP1(r,r) of the medium between an input virtual transducer TV.sub.in (spatial position r.sub.in) calculated based on an input focusing of the experimental reflection matrix that creates an input focal spot around the first point P1 and an output virtual transducer TV.sub.out (spatial position r.sub.out) calculated based on an output focusing of the experimental reflection matrix that creates an output focal spot around the second point P2, and a second response EP2(r, r) of the medium between an input virtual transducer TV.sub.in (spatial position r.sub.in) calculated based on an input focusing of the experimental reflection matrix that creates an input focal spot around the second point P2 and an output virtual transducer TV.sub.out (spatial position r.sub.out) calculated based on an output focusing of the experimental reflection matrix that creates an output focal spot around the first point P1.

(87) In addition, the mean correlation coefficient is calculated for distance coordinate values r of modulus greater than a predetermined resolution w.sub.d(r) (as shown in FIG. 6) and/or being calculated for a range of angle values or for a predetermined angle value .sub.d.

(88) The predetermined resolution w.sub.d(r) is advantageously a value greater than half the resolution w(r). Preferably, the predetermined resolution w.sub.d(r) is a value greater than one time, two times, or three times the resolution w(r), in order to better exclude the values of the peak, as can be seen in FIG. 6. Thus, the mean correlation coefficient according to this variable indicates that the correlation coefficient is averaged for values of the distance coordinate r far from the peak, and therefore for virtual transducers far from the central point PC, which in effect makes it possible to test the acoustic reciprocity of this central point PC, i.e. the symmetry of the response matrix REP(r,r) relative to the abscissa axis (see an example of such a matrix in FIG. 5). Although the values beyond the predetermined resolution are small, comparison by a mean correlation of the symmetry allows reliably estimating a level of symmetry (r) of the central point PC.

(89) Advantageously, the mean correlation coefficient is calculated for a range of angle values or for a predetermined angle value .sub.d. Thus, the mean correlation coefficient according to this variable indicates that the correlation coefficient is averaged for one or more angle values , which makes it possible to test the angular symmetry or based on an angular sector the level of symmetry around the central point PC, i.e. the reciprocity of this central point PC.

(90) Advantageously, the mean correlation coefficient is calculated for distance coordinate values r of modulus greater than a predetermined resolution w.sub.d(r) and for a range of angle values . Thus, the correlation is averaged over distance coordinate values for which the single scattering contribution is zero, the latter only appearing for a modulus of the distance coordinate that is below the predetermined resolution.

(91) This also makes it possible to obtain a more robust and stable estimate of the level of symmetry (r).

(92) For example, the level of symmetry (r) can be calculated by the following correlation formula:

(93) ( r ) = .Math. R E P 1 ( r , r ) .Math. REP 2 ( r , - r ) .Math. ( | r | > w d ( r , ) , ) .Math. .Math. "\[LeftBracketingBar]" REP 1 ( r , r ) .Math. "\[RightBracketingBar]" 2 .Math. .Math. .Math. "\[LeftBracketingBar]" REP 2 ( r , - r ) .Math. "\[RightBracketingBar]" 2 .Math. ( Eq . 5 )
or for example by the following correlation formula:

(94) ( r ) = .Math. Re [ REP 1 ( r , r ) .Math. REP 2 ( r , - r ) ] .Math. ( | r | > w d ( r , ) , ) .Math. .Math. "\[LeftBracketingBar]" REP 1 ( r , r ) .Math. "\[RightBracketingBar]" 2 .Math. .Math. .Math. "\[LeftBracketingBar]" REP 2 ( r , - r ) .Math. "\[RightBracketingBar]" 2 .Math. ( Eq . 6 )
in which Re[.Math.] is the real part mathematical operator, |.Math.| is the modulus mathematical operator, <.Math.> is the mean mathematical operator, this operator being able to be applied according to one or more variables (in the above, for example, distance coordinate values greater than a predetermined resolution and angle values ), and * is the complex conjugation operator.

(95) Several correlation formulas can be used, and the specialist will modify this definition according to his or her requirements and the characteristics of the medium observed.

(96) Generally speaking, the level of symmetry (r) is close to zero (0) if the propagation of the ultrasonic waves does not behave reciprocally around the central point PC, and the level of symmetry (r) is close to one (1) if the propagation of ultrasonic waves behaves symmetrically or reciprocally around the central point PC.

(97) Thus, this level of symmetry tests the validity of the acoustic reciprocity for the portion of the signals corresponding to multiple scattering. This makes it possible to differentiate the multiple scattering from noise, as noise does not comply with the property of acoustic reciprocity.

(98) Calculating a First Multiple Scattering Indicator (r)

(99) A step of determining a first multiple scattering indicator (r) of the central point PC can then be performed, this multiple scattering indicator being for example calculated by the following formula:

(100) 0 ( r ) = ( r ) 1 - ( r ) ( Eq . 7 )
in which (r) is the level of symmetry defined at the central point PC of the medium.

(101) This first multiple scattering indicator (r) is zero if the level of symmetry is zero, and it tends towards infinity if the level of symmetry is close to one (1).

(102) FIG. 10 illustrates images of the first multiple scattering indicator (r) established for the three heterogeneous media of FIG. 7 (media A, B, and C, which correspond between FIGS. 7 and 10). These figures are graduated in decibels (logarithmic scale).

(103) A large value for this first multiple scattering indicator corresponds to a high multiple scattering proportion (light in the figure). A small value for this focusing criterion corresponds to little or insignificant multiple scattering (dark in the figure). Thus, for medium A, this first multiple scattering indicator shows a localized region of ultrasonic wave scattering behind one of the heterogeneity cylinders. For media B and C, this first multiple scattering indicator highlights the scattering in most of the volume of these media, either downstream from the meat portion of medium B, or downstream from the fatty tissues of the liver.

(104) Calculating an Afocal Intensity I.sub.off(r)

(105) It is then possible to carry out a step of determining an afocal intensity I.sub.off(r) of the central point PC, this afocal intensity being a mean of a squared modulus of the responses REP(r,r), the mean being calculated for distance values of modulus greater than a predetermined resolution w.sub.d(r) (i.e. |r|>w.sub.d(r)) and/or calculated for a range of angle values or for a predetermined angle value .sub.d, for example calculated by the following formula:
I.sub.off(r)=custom character(I(r,r)custom character.sub.(|r|>w.sub.d.sub.(r,),)(Eq. 8)
with
I(r,r)=|REP(r,r)|.sup.2(Eq. 9)
in which I(r,r) is the intensity of the point, i.e. the squared modulus of the response at this point, |.Math.| is the modulus mathematical operator, and <.Math.> is the mean mathematical operator, this operator being able to be applied according to one or more variables (for example in the above, distance coordinate values greater than a predetermined resolution and angle values ).

(106) The afocal intensity I.sub.off(r) is used to characterize the incoherent energy in the ultrasound image.

(107) The afocal intensity I.sub.off(r) results from the contribution of echoes originating from multiple scattering phenomena (which are reciprocal) or coming from electronic noise (which are random). It is therefore possible to subdivide the afocal intensity I.sub.off(r) into multiple scattering intensity and noise intensity by using the level of symmetry (r).

(108) Calculating a Multiple Scattering Intensity I.sub.M(r)

(109) A step of determining a multiple scattering intensity I.sub.M(r) can then be performed, this multiple scattering intensity being the product of the level of symmetry (r) and the afocal intensity I.sub.off(r), i.e.:
I.sub.M(r)=(r).Math.I.sub.off(r)(Eq. 10)
Calculating a Noise Intensity I.sub.N(r)

(110) A step of determining a noise intensity I.sub.N(r) can then be performed, this noise intensity being the product of one minus the level of symmetry (r) and the afocal intensity I.sub.off(r), i.e.:
I.sub.N(r)=(1(r)).Math.I.sub.off(r)(Eq. 11)
such that we have the following relation:
I.sub.off(r)=I.sub.M(r)+I.sub.N(r)(Eq. 12)

(111) Then, the first multiple scattering indicator (r) can also be calculated by a ratio between the multiple scattering intensity I.sub.M(r) and the noise intensity I.sub.N(r), i.e. according to the following formula:

(112) ( r ) = I M ( r ) I N ( r ) ( Eq . 13 )

(113) The proportion of multiple scattering in the ultrasound image being low due to the antenna gain obtained by the double focusing process, this indicator makes it possible to identify the regions in which the multiple scattering proportion in the ultrasound image becomes predominant over the electronic noise.

(114) For example, for the images of FIG. 10, obtained with an angle of zero (=0), the following observations can be made.

(115) In image A of FIG. 10, one can see that the increase in scatterer density in the disc (located at spatial position [X,Z]=[5,40] mm) highlights an increase in the first multiple scattering indicator (r). Note that this appears for a depth greater than that of the disc. This phenomenon arises from the fact that the echoes resulting from multiple scattering phenomena have a time-of-flight which is no longer proportional to the depth of the last scatterer having interacted with the multiply-scattered wave. Consequently, this criterion makes it possible to probe the multiple scattering which takes place upstream of the point considered.

(116) In image B of FIG. 10, at a depth Z of 38 mm, the image of the first multiple scattering indicator (r) shows a sharp increase in this first multiple scattering indicator. These signals arise from a double reflection of the echoes between the surface of the Phantom and the surface of the probe, which is parallel.

(117) In image C of FIG. 10, there is a lateral pattern that repeats axially. This phenomenon arises from echoes which have been successively reflected between two parallel walls of tissue.

(118) Calculating a Confocal Intensity I.sub.on(r)

(119) A step of determining a confocal intensity I.sub.on(r) can then be performed, this confocal intensity then being the value of a squared modulus of the response REP(r, r=0) for a distance coordinate value r of zero, i.e. |r|=0, meaning for a point in the medium for which the first point P1, the second point P2, and the central point PC are coincident.

(120) Calculating a Single Scattering Intensity I.sub.S(r)

(121) A step of determining a single scattering intensity I.sub.S(r) can then be performed, this single scattering intensity being calculated based on the following equation:
I.sub.on(r)=I.sub.S(r)+2I.sub.M(r)+I.sub.N(r)(Eq. 14)

(122) Note that the factor of 2 in this equation comes from the phenomenon of coherent backscattering.

(123) Calculating a Second Multiple Scattering Indicator (r)

(124) A step of determining a second multiple scattering indicator (r) can then be performed, this second multiple scattering indicator being calculated by the following formula:

(125) ( r ) = I M ( r ) I S ( r ) ( Eq . 15 )

(126) This second multiple scattering indicator (r) makes it possible to compare the proportion of multiple scattering to the proportion of single scattering. FIG. 11 illustrates images of the second multiple scattering indicator (r) established for the three heterogeneous media of FIG. 7 (media A, B, and C, which correspond between FIGS. 7 and 10). These figures are graduated in decibels (logarithmic scale).

(127) A large value for this first multiple scattering indicator corresponds to a large proportion of multiple scattering compared to single scattering (light in the figure). A low value for this second multiple scattering indicator corresponds to little multiple scattering compared to single scattering (dark in the figure). Thus, this second multiple scattering indicator shows that medium C is highly scattering.

(128) Although the images of FIG. 11 highlight the same phenomena observed in the corresponding images of FIG. 10 but with a different contrast, this second multiple scattering indicator makes it possible to determine a new criterion sensitive to the phenomenon of multiple scattering, independently of the attenuation of the medium. In effect, we are studying the echoes which have the same propagation time inside the medium: they therefore have been attenuated in the same manner.

(129) Optimizing the Ultrasound Image Calculation

(130) One possible use of the above calculations is to optimize the calculation of the ultrasound image.

(131) The ultrasound image is for example calculated by successive input and output focusing for all points in the medium. The set of points on the abscissa axis of FIG. 5 corresponds to a row of this ultrasound image calculated at the depth z used to form the image of this FIG. 5.

(132) However, this ultrasound image calculation strongly depends on the assumption of a homogeneous medium in which the speed of sound (the speed of propagation of ultrasonic waves) is well-known and constant. If this assumption is incorrect, the focusing delay rules do not correspond to the medium considered and the focusing is imperfect. The resolution of the ultrasound image is then degraded. Other focusing calculation parameters may possibly influence focusing.

(133) An image optimization step can then be carried out in which the focusing criterion F(r) is calculated at a plurality of points in the medium (for example a predetermined region of the image), and at least one calculation parameter for the input focusing and/or output focusing is optimized by minimizing or maximizing a mean of said focusing criterion F(r) for said plurality of points.

(134) For example, said at least one calculation parameter comprises at least the speed of sound in the medium. Optionally, this calculation parameter is the speed of sound.

(135) In the case of optimizing the focusing criterion F(r) based on a possible range of the speed of sound c, a curve is obtained similar to one of those presented in FIG. 12, and for which the focusing criterion F(r) exhibits a maximum for an optimum speed of sound c.sub.opt. This optimum speed of sound c.sub.opt corresponds in fact to a speed of sound integrated over the entire thickness of the medium traversed (between the probe and the focal plane of the predetermined region).

(136) In other words, the focus quality and therefore the focusing criterion is maximal when the speed of sound model used to carry out the focusing coincides with the true speed of sound of the medium.

(137) To improve the estimation of the speed of sound in the medium, it is necessary for example to introduce a multi-layered model of the medium, each layer being assumed to be homogeneous with a speed of sound c.sub.i of this layer. The thickness of each layer must be estimated, either with prior knowledge of the medium, or based on a first ultrasound image. By optimization of the focusing criterion, one therefore estimates the speed of sound in the outermost layer with a homogeneous model, then by optimization the next layer one estimates the speed of sound in the next layer with a two-layer model, and so on.

(138) FIG. 12 illustrates the optimizations carried out for a same predetermined region of medium B, using: a) a single-layer model in curve B1 with optimum speed of sound c1.sub.opt, b) a two-layer model in curve B2 with optimum speed of sound c2.sub.opt, c) a three-layer model in curve B1 with optimum speed of sound c3.sub.opt,

(139) The optimization method presented above and the multilayer modeling of the medium therefore make it possible to refine the speed of sound estimates, and therefore make it possible to determine the depthwise evolution of the speed of sound in the medium.