Method for determining a field of an ultrasonic wave through a plurality of walls

Abstract

The aim of the present invention is to provide a semi-analytical method for accurately simulating the transmission of an ultrasonic field through a cranial wall by means of algorithms having computation times that are compatible with a preoperative context. More generally, the invention consists in simulating an ultrasonic field through one or more walls in order to focus a beam at a point located behind the walls or else to direct a beam in a predetermined direction behind that wall.

Description

DESCRIPTION Title of the invention: Method for determining a field of an ultrasonic wave through several walls

[0001] The invention relates to the field of simulating the propagation of ultrasonic waves in heterogeneous and / or anisotropic solid media of complex geometry. The invention applies in particular to transcranial ultrasound therapy but also to ultrasound imaging.

[0002] More specifically, it relates to a method of simulating an ultrasonic wave emitted by an ultrasonic transducer through the cranial wall in order to control the transmitted ultrasonic beam, and in particular the focal zone, before generating it.

[0003] Cerebral or transcranial ultrasound therapy (TUS) involves the use of focused ultrasound to treat certain diseases. Curing brain diseases remains very challenging, mainly due to poor access of pharmacological agents to the brain due to the blood-brain and blood-tumor barriers (BBB and BHT). Focusing low-intensity pulsed ultrasound into the brain, combined with the injection of microbubbles, significantly increases drug delivery to brain tissue, with an established therapeutic effect in numerous animal models as well as in patients. This permeabilization of the BBB or BHT is local and transient provided that the beam intensity is well controlled because the implosion of microbubbles could lead to hemorrhages.Therefore, in order to ensure the safety of the protocol, it is necessary to precisely control the dosimetry of the acoustic field after crossing the bone wall.

[0004] High-intensity focused ultrasound (HIFU) therapy is an innovative approach for the treatment of intracranial tumors. In the near future, this therapy could be an attractive alternative to brain surgery and radiotherapy. HIFU therapy induces a temperature rise at the focal point of the ultrasound beam, and a defined tissue volume can be thermally destroyed (thermal ablation). For a given frequency and acoustic parameters, the ablation volume depends on the temperature rise and the exposure time. (thermal dose). This technique requires perfect control of the focal zone of the transmitted ultrasound beam, in order to target the area to be destroyed (tumor) while ensuring its safety on surrounding healthy tissue.

[0005] In the different applications considered above as well as for ultrasound imaging applications, it is necessary to precisely control the propagation of ultrasonic waves especially since it varies according to the geometry and mechanical properties of the cranial wall.

[0006] Indeed, incident ultrasound beams can be deflected, attenuated or aberrated, or even a combination of two or three of these cases at the same time. As a result, the skull is then considered a complex environment that strongly disrupts the focusing of ultrasound for therapy, or image reconstruction in the case of imaging. This problem, inherent to the skull, is mainly due to a geometry and local acoustic properties that are very different in the same individual, but also due to a variation of these same properties depending on the individual. Personalized modeling of the skull then takes on its full importance to simulate the propagation of ultrasound and improve therapy and imaging tools.

[0007] Usually, simulations of ultrasound fields generated towards transcranial walls are performed by numerical methods which are considered to be more accurate, although more expensive in execution time and memory space, than semi-analytical methods.

[0008] Reference [1] lists, in a bibliographic study, the main methods of simulating the propagation of ultrasonic fields through a cranial wall, which are grouped into two main categories.

[0009] The first category of methods concerns numerical methods, for example finite element, finite difference or boundary element methods, which are considered as the reference. Indeed, since they consist of an exact numerical resolution of the equations of propagation of ultrasonic waves, these methods make it possible to simulate all the complex physical phenomena involved. However, they remain in most cases expensive in terms of calculation time and memory space occupied. Furthermore, they have strong constraints on the calculation area. Indeed, this must be sufficiently spatially sampled, to meet numerical convergence criteria, and it must contain both the source term, i.e. the ultrasound probe, and the area of ​​interest. Finally, the implementation generally requires significant expertise on the part of the user, in particular through the use of scene meshing tools, including both the probe and the skull. An example of a comparative study of the main numerical simulation tools is given in reference [2].

[0010] The second category concerns so-called semi-analytical methods, in particular ray tracing methods, which can also be implemented for transcranial ultrasound simulations. These methods are faster, less expensive in terms of memory load, and do not involve constraints on the computing area. On the other hand, since they are based on asymptotic approximations, they do not allow to simulate satisfactorily all the complexity of the physical phenomena involved during propagation. Nevertheless, recent studies have shown that these semi-analytical methods allow to obtain results as precise as numerical methods.

[0011] Among the semi-analytical methods, the so-called brush method, based on ray theory, is a fast and reliable method but is very sensitive to the regularity of the surfaces describing the different propagation media. Indeed, discontinuities can lead to over-intensity zones (caustics) or to zones where no ray can propagate (Fresnel shadow). These limits of validity are particularly reached when taking into account organs or industrial components of complex 3D geometry, described by meshes.

[0012] The implementation of the brush method requires defining the ultrasonic path between each source point and the point of interest, and respecting the laws of refraction and / or reflection at each interface depending on the chosen propagation mode. For configurations involving organs or industrial components of complex 3D geometry, described by meshes, the step of determining these paths is not trivial.

[0013] The state-of-the-art brush method consists of blindly searching for these paths starting from the point of interest and sending paths in all directions of space to keep only those reaching the emitting surface of the probe. This method is very expensive in terms of computation time and the complexity increases with the number of interfaces to manage. It does not ensure that all ultrasonic paths have been found exhaustively. The calculation of the surface area carried by each brush at the source point remains approximate and requires complex matrix calculations. Finally, it does not ensure regular sampling of the emitting surface, shadow or overlapping areas can be obtained. All these limitations lead to uncertainties in the ultrasonic field calculation.

[0014] The object of the present invention is to propose a semi-analytical method for precisely simulating the transmission of an ultrasound field through a cranial wall by means of algorithms having compatible calculation times in a preoperative context. More generally, the invention consists of simulating an ultrasound field through one or more walls to focus a beam at a point located behind the walls or to direct a beam in a predetermined direction behind this wall.

[0015] The invention is notably based on a different implementation of the brush method which makes it possible to significantly increase the reliability of the results and to very significantly reduce the calculation times for configurations involving complex 3D surfaces. It is based on spatial sampling of each ultrasonic element and on regularization and parameterization of the 3D surfaces allowing the implementation of an algorithm for calculating ultrasonic paths between each source point and each point of interest by a method for minimizing flight times.

[0016] The subject of the invention is a method, implemented by computer, for determining a field of an ultrasonic wave emitted by an ultrasonic transducer comprising at least one ultrasonic element towards an area of ​​interest, the transducer being separated from the area of ​​interest by several media separated by walls, the method comprising the steps of: - Receive a model of each wall in the form of a parameterized 3D surface obtained by approximation of multi-level B-splines, - Receive a decomposition of the surface of each ultrasonic element into a plurality of elementary surfaces each being described by a set of points comprising a central point and four points corresponding to the vertices of a square centered on the central point, - Determine an ultrasonic field transmitted by each ultrasonic element and for each point of interest in the area of ​​interest by means of the sub-steps of: • For each mode of propagation of the ultrasonic wave among several modes and for each source point of the set of points characterizing an elementary surface among the set of elementary surfaces describing each ultrasonic element, determine a path between the source point and the point of interest so as to minimize the flight time of the ultrasonic wave between these two points, determine the attenuation of the ultrasonic wave on this path, calculate the impulse response of the ultrasonic field for each elementary surface from the paths and attenuations determined for the set of associated points, • Sum all the calculated impulse responses, • Determine the ultrasonic field transmitted by the ultrasonic element from the result of the sum.

[0017] According to a particular aspect of the invention, the ultrasonic transducer comprises several ultrasonic elements, the method further comprising a sum of the ultrasonic fields of each element weighted by a predefined amplitude law, each ultrasonic field being delayed via a predefined delay or phase law.

[0018] In an alternative embodiment, the method according to the invention comprises a prior step of decomposing the surface of each ultrasonic element comprising at least: - Decompose the surface into a plurality of contiguous elementary surfaces, each elementary surface being centered on a central point, - For each of the central points, form a square centered on the central point in a plane tangent to the surface, - Characterize the elementary surface by the set of points including the central point and the four vertices of the square,

[0019] In an alternative embodiment, the method according to the invention comprises a prior step of determining a model of each wall comprising at least: - Acquire an image of the walls, - Extract wall meshes in the form of a point cloud, - Apply a multi-level B-spline approximation to the point cloud in order to obtain a smoothed 3D parameterized surface

[0020] According to a particular aspect of the invention, the media separated by walls are anisotropic and the step of determining a path between the source point and the point of interest is carried out by taking into account the anisotropy of each media crossed.

[0021] According to a particular aspect of the invention, the step of determining a path between the source point and the point of interest is carried out by means of an algorithm for minimizing the sum of the segments making up the path, each segment corresponding to a medium crossed, the sum being weighted by the slowness of each medium.

[0022] According to a particular aspect of the invention, the propagation modes of the ultrasonic wave comprise at least one propagation mode involving reflections of the ultrasonic wave on at least certain walls and taking into account the reflection coefficient of the path coming from the central point of each elementary surface in the calculation of the impulse response of the ultrasonic field.

[0023] According to a particular aspect of the invention, the decomposition of a surface of an element into elementary surfaces is carried out using a Voronoi diagram.

[0024] According to a particular aspect of the invention, the calculation of the impulse response of the ultrasonic field comprises at least the calculation of the amplitude of the field and the minimum and maximum limits of the flight times corresponding to the calculated paths.

[0025] According to a particular aspect of the invention, the area of ​​interest is a transcranial region and the walls are the walls of a skull.

[0026] According to a particular aspect of the invention, the step of determining an ultrasonic field comprises at least the step of convolving the result of the sum with a reference signal associated with the ultrasonic element or of extracting the amplitude and the phase of a predefined frequency component of each impulse response to obtain an ultrasonic field transmitted by said element.

[0027] The invention also relates to a computer program comprising instructions for executing the method according to the invention, when the program is executed by a processor.

[0028] The invention also relates to a recording medium readable by a processor on which is recorded a program comprising instructions for executing the method according to the invention, when the program is executed by a processor.

[0029] The invention also relates to an ultrasonic wave field simulator comprising a processor and a memory configured to execute the steps of the method according to the invention and a display device for displaying the calculated field in an area of ​​interest.

[0030] Other features and advantages of the present invention will become more apparent upon reading the following description in relation to the following appended drawings.

[0031] Fig. 1 represents a flowchart detailing the steps of implementing a method for determining an ultrasonic field emitted through at least one wall according to one embodiment of the invention.

[0032] Fig. 2 represents a flowchart illustrating a step of calculating ultrasonic paths and attenuation according to an embodiment of the invention,

[0033] Fig. 3 represents a flowchart illustrating a step of calculating an impulse response of an emitted ultrasonic wave, according to an embodiment of the invention,

[0034] Fig. 4 represents a diagram illustrating an example of decomposition of an ultrasonic probe into several elements,

[0035] Fig. 5 represents a diagram illustrating the principle of the path calculation method of Figure 2,

[0036] Fig. 6 shows a diagram illustrating the principle of calculating an ultrasonic field by the brush method,

[0037] Fig. 7 shows an example of an impulse response of an ultrasonic field,

[0038] Figure 1 illustrates, on a flowchart, the main steps of a method for determining an ultrasonic field emitted from a source point of a surface of an ultrasonic probe to a target point located behind one or more walls. In the remainder of the description, the invention is presented in the context of a cranial wall for therapy or transcranial imaging applications but the invention is not limited to this type of wall and can be applied to other types of object walls, in particular for non-destructive testing applications.

[0039] The method comprises a preprocessing phase 100 followed by a field calculation phase 110.

[0040] The preprocessing phase 100 takes as input a model of a skull 101 which is for example obtained via an X-ray medical imaging technique and a model of the geometry of an ultrasound probe 102 intended to be used to generate an ultrasound field.

[0041] In step 103, the acoustic properties of the skull are extracted from the obtained model. This step is carried out, for example, using the extraction techniques described in references [3] or [4].

[0042] In step 104, meshes of the skull walls are extracted, for example using a technique for decomposing a surface into meshes. The meshes are, for example, made up of 3D point clouds.

[0043] The skull description elements obtained in steps 103 and 104 are used to model the external surfaces constituting the walls of the skull in step 105.

[0044] In order to model the different walls of the skull by smooth surfaces, a B-spline approximation method is used in step 106, for example from the techniques described in reference [5].

[0045] An advantage of this method is that it consists of an approximation and not an interpolation, which makes it possible to be robust to possible outliers in the point cloud obtained at extraction 104. This method makes it possible to obtain smoother surfaces than with a technique based on interpolation, which makes it possible to improve the convergence of the ultrasonic path calculation algorithms described later. It also has the advantage of reducing noise in the point cloud.

[0046] The approximation method is for example a so-called multilevel B-spline or “MBA” method. At the end of step 106, each wall of the skull is modeled via a 3D surface which constitutes an interface between two media of index i and i+1 and which is described by a function ft u^Vi) = (x,y,z) where (u t , Vt) corresponds to the parametric coordinates of the surface described by means of splines. In the case of the MBA method, we have (uj.Vj) = (xy).

[0047] The interfaces defining the walls of the skull are at least two in number but may be more than two.

[0048] The pre-processing phase also includes a modeling phase of the ultrasonic probe. The probe used can be a single or multi-element ultrasonic transducer. The invention is suitable for any type of 2D or 3D surface geometry for each element of the transducer.

[0049] In step 106, each element of the emitting surface of the probe is spatially sampled to obtain a set of elementary surfaces, characterized by a central point and a surface, distributed uniformly over the entire surface of the element.

[0050] The number of elementary surfaces is an input parameter of the method. This sampling can be obtained, for example, from a Voronoi diagram. The sum of the elementary surfaces is equal to the total surface area of ​​the element.

[0051] Figure 4 illustrates an example of element 400 of a transducer whose shape is circular. The element 400 is decomposed into several elementary surfaces 401 with centers C.

[0052] In step 107, the coordinates of 4 vertices of a square with a center equal to point C, which is the center of the elementary surface, are determined for each elementary surface. The virtual square thus formed has a surface A with the same area as that of the element 401. In the case of a surface shaped in three dimensions, the virtual square is formed in a plane (ït, v) of a direct orthonormal reference frame (ït, v, n) with n a vector normal to the elementary surface 401.

[0053] This step is illustrated on the right of figure 4 which shows, for the example of element 401, the virtual square with center C and whose four vertices are defined, in a direct orthonormal reference frame (ït, v, n) by the coordinates:

[0058] These four points are thus defined for each elementary surface to serve as source points in the calculation of paths prior to field calculations. This modeling of surface elements makes it possible to avoid the complex and time-consuming matrix calculations present in the state-of-the-art brush method, to determine the coefficient of divergence and the duration of the time slot intended for calculating the impulse response associated with this surface element. This description of each surface element in the form of a square surface also makes it possible to address any probe geometry.

[0059] At the end of the pre-processing phase 100, we therefore obtain a model of the walls of the skull via three-dimensional parametric functions on the one hand and on the other hand the coordinates of five points (the center and the four vertices of the virtual square) characterizing each elementary surface resulting from the decomposition of an element of the ultrasound probe.

[0060] All of these parameters are provided as input to a second phase 110 of calculating ultrasonic fields. The particular example in Figure 1 is given for calculating the ultrasonic field corresponding to an ultrasonic element of a transducer. If the transducer comprises several ultrasonic elements, an additional combination step must be applied from the ultrasonic fields calculated for each of the elements. This step will be described in more detail later.

[0061] The method for calculating the object ultrasonic field of the second phase 110 comprises three iteration loops. The first loop aims to traverse all the points of the calculation zone, i.e. the points of the brain in the case of a transcranial application or the points of an area of ​​an object to be controlled in the case of a non-destructive testing application. For each target point of the calculation zone, an ultrasonic field emitted through the walls modeled in step 105 is calculated.

[0062] This calculation requires a first calculation of paths by minimizing flight time 111 and a second calculation of impulse response 112, these two steps being iterated on a second iteration loop going through several modes of propagation of the ultrasonic waves and a third iteration loop going through each point characterizing each elementary surface of each element of the probe. These points correspond to those determined in step 107.

[0063] Figure 2 details the calculation 111 of ultrasonic paths and Figure 3 details the calculation 112 of the associated impulse response.

[0064] The goal of 1 1 1 ultrasonic path calculation is to find the shortest path between a point A on the probe and a target point B located on the other side of the cranial wall modeled by several surfaces.

[0065] This principle is illustrated in Figure 5 for an example of a wall modeled by two surfaces fi and Ï2 separating three media characterized by ultrasonic wave propagation speeds respectively equal to Ci, C2 and C3. The aim of the calculation step 1 1 1 is to find the intersection points and l 2 of the ultrasonic path emitted from A to B which allows the shortest flight time to be obtained.

[0066] Let us note 11 = (x1,y1,z1) the interface point 1 and I2 = (x2,y2,z2) of the interface 2 through which the ray passes, we have (x^y^ z and (x 2 ,y 2 , z 2 ) = f 2 (u 2 ,v 2 ) because the interfaces are described by splines.

[0067] The 4 parameters to vary are therefore:

[0069] The objective of the 1 1 1 path calculation is to minimize the travel time given by the relation:

[0070] The path calculation step 11 1 receives as input a source point 210, that is to say a point on the surface of the emitting probe and a calculation point 21 1 , that is to say a point in the calculation zone. It also receives as input the interfaces 213 modeling the wall to be crossed as well as the speeds 212 associated with each medium separated by an interface.

[0071] Step 201 consists of a ray tracing between the source point A and the calculation point B so as to define the points of intersection of the ray with the different interfaces. The intersection points belong to the approximate surfaces obtained from splines. Step 201 can be carried out by means of a ray tracing algorithm as described for example in the Applicant's patent application FR31 14032.

[0072] Step 202 then consists of finding the coordinates of the intersection points which minimize the flight time given by the previous relation.

[0073] This optimization step can be carried out by any suitable optimization algorithm, for example a quasi-Newtonian minimization algorithm such as the L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shannon) algorithm, or by the method described in patent application FR31 14032.

[0074] The general objective of the optimization step 202 is to search for the path between the source point A and the calculation point B through N interfaces such that the time-of-flight function t is minimal. As each medium is homogeneous with slowness s f = - the flight time is the sum of the flight times in each medium crossed by the wave. To find this minimum, this amounts to finding the points of intersection of the ultrasonic path on each of the interfaces. We therefore use a minimization algorithm such as the L-BFGS algorithm to solve the following problem: min t(

[0075] fi(ui, vt) is the parametric surface that describes the interface between the media of index i and i+1.

[0076] For anisotropic media, the slowness s depends on the direction of propagation and the previous relationship is then written:

[0077]

[0078] Taking into account anisotropic media in the calculation of the path of the ultrasonic wave makes it possible to better characterize the object crossed by the wave, in particular the human skull which has properties specific to each individual.

[0079] For the anisotropic media considered, the cost function to be minimized and characterizing the flight time of the ultrasonic path mint( > must be parameterized to account for the dependence of the slowness on the propagation direction. This parameterization can be carried out from an analytical function S] describing the variation of the slowness with the propagation direction, expressed in 3D, which must be evaluated at each iteration. For optimization purposes, this parameterization can also be obtained from a tabulation of the slowness as a function of the propagation direction, expressed in 3D, established in a preliminary step.

[0080] After optimization, we obtain a description 203 of the ultrasonic path between the source point and the calculation point.

[0081] According to another aspect of the invention, the attenuation of the wave through each medium crossed is also calculated in step 204.

[0082] To calculate the total attenuation on a path, we add the attenuations on each of the portions of the path in each medium crossed. If each medium has an attenuation coefficient a t at frequency f o i , and an attenuation exponent fa, the total attenuation is obtained via the following relation:

[0084] The attenuation values ​​of each medium are input data of the method.

[0085] In one embodiment of the invention, the path calculations 111 are performed for different propagation modes. In particular, certain propagation modes involve reflections on the internal interfaces of the heterogeneous structure crossed (for example the skull).

[0086] In this case, the reflections are considered as additional path segments. We then apply the same method as that described previously by adding as many additional interfaces as there are reflections considered which correspond to multiple propagations of the wave in the same medium.

[0087] Depending on other propagation modes, the values ​​of the speeds and / or attenuation of the media crossed may differ from one mode to another.

[0088] The paths thus calculated in step 111 for each pair of points (source point, calculation point) and each propagation mode are then transmitted to a step 112 for calculating the impulse response of the emitted ultrasonic field.

[0089] More precisely, the paths are calculated for each center of an elementary surface of an element of the probe as well as for the four vertices of the virtual square formed from this center (as explained in step 107).

[0090] Figure 3 details the implementation of step 112 of impulse response calculation.

[0091] This calculation step receives as input, for each elementary surface of an element of the probe, the path coming from the central point 301 as well as the paths 302 coming from the four vertices of the virtual square.

[0092] The calculation of the impulse response 112 is partly based on the so-called brush method described in reference [6].

[0093] First, the principle of the brush method is described.

[0094] To reconstruct the emitted wavefront and therefore calculate the Rayleigh integral associated with the ultrasonic field, it is possible to use the brush method.

[0095] Physically, a brush represents a flux of acoustic energy originating from the point energy source O and propagating towards the surface S. Since O is the only energy source considered in the model, the law of conservation of energy tells us that the total energy flux crossing the surface of the brush is zero. Thus, the surface energy crossing the surface S is inversely proportional to its area.

[0096] The central ray of the brush follows the Fermat trajectory, obtained from the path calculation step 111 . The paraxial ray follows a different direction from that of the central ray, but which is considered close enough for the paraxial approximation to be valid (the paraxial approximation being an approximation of order 1 , the error increases as one moves away from the axis). It is characterized by four differential quantities expressed in a frame included in the plane orthogonal to the central axis of the brush, allowing to define a paraxial radius of the brush by a vector with 4 coordinates: paraxial ray and the surface S of the brush. dSx and dSy give the projection of the slowness vector of the paraxial ray onto the surface S.

[0098] In the context of paraxial approximation, the evolution of a beam vector in a homogeneous medium can be characterized by a 4x4 matrix that depends on the propagation medium. Similarly, 4x4 matrices characteristic of the interfaces between two materials can be defined. These matrices are called propagation matrices and are usually denoted using a block notation: L = ) where the matrices A, B, C, and D are 2x2 matrices.

[0099] Considering a beam passing through a sequence of homogeneous media and interfaces characterized by the propagation matrices (D) for 1 < i < n, the resulting beam can be written in the form: = (n =1 L f ) ijj

[0100] Since a heterogeneous medium can be modeled by a sequence of homogeneous media and the interfaces separating them, these matrices are sufficient to describe the evolution of a beam when passing through a heterogeneous medium.

[0101] The use of brushes for ultrasonic field calculation allows to simulate the spatial propagation of an ultrasonic wave. In order to obtain the final amplitude of the brush, its divergence must be taken into account. Indeed, the law of conservation of energy inside the brush means that the energy flux decreases with the opening of the brush, which is, in the framework of the approximation of a brush with a small opening, distributed uniformly on the base surface S.

[0102] The ratio between the final surface intensity Is and the intensity per unit of initial solid angle IQ is therefore written as the ratio between the initial solid angle dQ and the final base surface

[0103] In addition to the amplitude factor related to the divergence of the brush during its propagation, each interface crossed reduces the energy of the brush since different modes are generated at each interface. The energy distribution between the different propagation modes (transmitted, reflected and evanescent for each type of polarization) is obtained from the continuity conditions at the interfaces and results in the calculation of a Fresnel coefficient.

[0104] To obtain the value of the amplitude transmission coefficient, directly applicable to the value of the displacement and therefore of the amplitude of the ultrasonic field, the ratio between the brush sections before and after the interaction must be considered.

[0105] For each brush emitted from the considered field point and whose final surface intersects the surface of the ultrasonic transducer, it is possible to extract an elementary impulse response. The summation of these elementary responses constitutes the global impulse response which is then used to determine the amplitude value of the field over time. The elementary impulse response of a brush is assimilated to a time step whose amplitude, time spread and central time are extracted from the final characteristics of the brush. The amplitude value corresponds to an energy. Consequently, it is distributed over the time spread of the brush, which amounts to adding an inverse factor of the time spread to the amplitude value.

[0106] The amplitude A is therefore given by the relation A = A ° XDFxFAXds

[0107] With A0 the initial amplitude, DF the divergence factor, T the amplitude transmission coefficient, dS the surface of the brush at arrival, At the time spread.

[0108] The time spread At is calculated as the maximum path difference evaluated on the brush, as a function of the angle of incidence of the brush on the translator and its spatial extent.

[0109] The central time is equal to the flight time of the central ray of the brush.

[0110] For a multi-element translator using delay or amplitude laws, the center time and amplitude of each brush are modified consistently with the intersected translator element.

[0111] Since Fresnel coefficients can be complex, the magnitude of the displacement is itself potentially complex. In addition, the displacement being a vector field, the calculated response is a signal with 3 complex dimensions, or 6 real dimensions.

[0112] Once the elementary impulse responses are calculated for a given field point, we sum them into a single impulse response. This corresponds to the field radiated by the translator towards the field point when each element of the translator emits a time-domain Dirac pulse.

[0113] Therefore, given an input signal S o to which the ultrasonic translator is subjected, the ultrasonic field at a field point is obtained by convolution of the complex impulse response signals with the input signal. For each x, y and z component of the displacement, the real part is convolved with the reference signal while the complex part is convolved with the Hilbert transform of the reference signal. The calculation of each scalar component si(t) of the field is therefore written, denoting Rli(t) the impulse response: si(t) = Rli(t) * SO(t)

[0114] The paraxial brush model is based on propagation matrices to provide a linear approximation of paraxial rays around the central ray of a brush.

[0115] The invention proposes to use an approximation of the paraxial rays included in the brush using interpolation. To do this, we simulate the propagation of the four rays r 0 , h, r 2 and r 3 delimiting the brush in addition to propagating the central ray rc, as shown in Figure 6. As explained previously, the four rays r 0 , r 1 ; r 2 and r 3 correspond to the four paths from the four vertices of the virtual square centered on the center of an elementary surface obtained by decomposition of an element of the ultrasonic transducer.

[0116] This approach makes it possible to avoid the calculation of propagation matrices along the path of the brush while making accessible the calculation of the quantities necessary for the integration of the elementary impulse response associated with the brush.

[0117] This approach allows the use of controlled, regular and homogeneous surface sampling of the emitting surface of each element. The position and surface covered by each brush is exact and controlled, which is not the case with the ray tracing starting from the calculation point and the estimation of the brush parameters by matrix calculation. The real emitting surface of each element is thus taken into account by the brush method.

[0118] The time window of the slot modeling the brush is determined so as to frame all the arrival times of the rays composing the brush. Thus, the limits t m in and t max of the time slot representing the impulse response of the brush are determined respectively as the minimum instant and the maximum instant among the set of calculated flight times.

[0119] Figure 7 represents the resulting time slot.

[0120] Using the four-ray brush model, the brush area is the area enclosed by the points from the four rays as shown in Figure 6. The initial solid angle is defined by the area of ​​the ray brush at unit distance.

[0121] The divergence coefficient can be calculated from these quantities and the initial slowness So of the axial ray via the following relation: DF = s ° ( — > dS

[0122] In this way, considering the cumulative amplitude transmission coefficient T A of the central ray and the initial amplitude A o, we can calculate in the same way as previously the amplitude of the displacement at the field point considered by distributing it over the time spread of the brush via the relation . _ A o vDFvT A vdS / i — . At

[0123] This simplified brush model allows reproducing an elementary brush impulse response without having to calculate paraxial matrices, removing the need for expensive matrix calculations (especially in anisotropic media) and replacing them with ultrasonic ray tracing operations. For an isolated brush, this approach requires a total of 5 ultrasonic rays to produce an elementary impulse response.

[0124] The impulse response calculation step 112 is therefore carried out on the basis of the brush method described above.

[0125] In step 303, a polarization calculation is performed from the central path 301.

[0126] In step 304, a cumulative transmission or reflection coefficient T is calculated. A (depending on the propagation mode) of the central path. In the case where the propagation mode chosen involves both transmissions and reflections on the different interfaces, the cumulative coefficient is calculated from the transmission and reflection coefficients on each of the interfaces considered.

[0127] In step 305, the flight times corresponding to the 5 journeys are calculated. 301,302 received at input. From these flight times, we deduce the temporal position and the duration of the time slot describing the impulse response.

[0128] In step 306, we calculate the solid angle formed by the four paths r 0 , n, r 2 and r 3 corresponding to the vertices of the virtual square, at the calculation point.

[0129] In step 307, the divergence coefficient DF is calculated from the solid angle via the relation

[0130] In step 308, the values ​​of t are calculated m in and t ma x in the manner described previously.

[0131] Finally, in step 309, the amplitude A is calculated via the formula given above.

[0132] The impulse response for a pair (source point, calculation point) and a propagation mode is given by the limits of the time slot, its amplitude and its polarization.

[0133] More precisely, we calculate the impulse response of the ultrasonic field for the elementary surface and the propagation mode considered by distributing, over the time slot obtained, the energy of the brush calculated from the attenuation coefficient, the reflection and / or transmission coefficients of the central ray and the divergence coefficient.

[0134] We then return to Figure 1. The calculation continues at step 113 by determining the overall size of the signal which is obtained by recording the minimum and maximum time instants among all the time slots calculated for the two iterations (propagation mode, source point). This step is necessary to allocate sufficient memory space for the output signal.

[0135] In step 114, we sum all the calculated impulse responses.

[0136] In step 115, a convolution of the impulse response with the reference signal or a monochromatic calculation is performed. In the most general case, the reference signal is the one used to generate the ultrasonic signal by the probe. In the case of the monochromatic calculation, one implementation consists of calculating only the amplitude and the phase of the spectrum measured at the desired frequency.

[0137] Finally, in step 116 we obtain the pressure field at any point in the calculation area for an ultrasonic element.

[0138] In the case where the transducer is multi-element, the algorithm 110 for calculating the ultrasonic field must be iterated for each element of the transducer. An additional step (not shown in Figure 1) must also be applied in order to combine the different fields emitted by each element into a global field associated with the transducer. For this, an amplitude law and a delay law predefined according to the intended application are taken into account in this combination. More precisely, this additional step consists of calculating the sum of the ultrasonic fields of each element, weighted from the amplitude law, each field being delayed from the delay law. In other words, the amplitude law and the delay law give the amplitudes (respectively delays) to be applied to each field corresponding to each element.

[0139] The invention has several advantages over the methods of the prior art. It is suitable for taking into account ultrasound probes having surfaces of any geometry, which presents a particular advantage in the context of transcranial applications for which the emitting surface of the probe can have a spherical profile or more generally shaped in 3D to adapt to the shape of the skull. Similarly, in the general case, each element has a specific projected surface and of complex geometry, for example of the polygonal type. Indeed, the decomposition into elementary surfaces of each element of any geometry coupled with the formation of four points of interest associated with the central point of the element makes it possible to adapt to any geometric shape.

[0140] Furthermore, the invention also makes it possible to take into account anisotropic media in the calculation of paths and to take into account the attenuation of each medium in the calculation of the ultrasonic field. Finally, better modeling of the surfaces defining the skull also allows an improvement in the precision and convergence time in the calculation of ultrasonic paths.

[0141] The invention may be implemented as a computer program comprising instructions for its execution so as to enable a simulation of the ultrasonic fields generated for each point of a targeted area. The computer program may be recorded on a recording medium readable by a processor. The simulator may comprise a processor and a memory for executing the method of generating a field of an ultrasonic wave according to the invention as well as a display device for displaying the field generated for an area of ​​interest, for example a transcranial area.

[0142] Such a simulation tool provides assistance to an operator prior to generating the ultrasonic field using a transducer so as to enable better calibration of the field before its generation or, in the case of a multi-element probe, to calculate the phase or delay laws, enabling correction of the aberrations undergone by the field when crossing the cranial wall.

[0143] References

[0144] [1] Célestine Angla, Benoit Larrat, Jean-Luc Gennisson, Sylvain Chatillon, Transcranial ultrasound simulations: A review, First published: 01 September 2022

[0145] [2] JF Aubry, Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models, The Journal of the Acoustical Society of America 152, 1003 (2022);

[0146] [3] J.-F. Aubry, M. Tanter, M. Pernot, J.-L. Thomas, and M. Fink, “Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans,” J. Acoust. Soc. Am., vol. 113, no. 1 , pp. 84-93, 2003, doi: 10.1121 / 1 .1529663.

[0140] [4] L. Marsac et al., “Ex vivo optimisation of a heterogeneous speed of sound model of the human skull for non-invasive transcranial focused ultrasound at 1 MHz,” Int. J. Hyperth., vol. 33, no. 6, pp. 635-645, 2017, doi: 10.1080 / 02656736.2017.1295322.

[0141] [5] B. G. Lee, J. J. Lee, and J. Yoo, “An efficient scattered data approximation using multilevel B-splines based on quasi-interpolants,” Proc. Int. Conf. 3-D Digit. Imaging Model. 3DIM, pp. 110-117, 2005, doi: 10.1109 / 3DIM.2005.18.

[0142] [6] H. Chouh, “Simulations interactives de champ ultrasonore pour des configurations complexes de contrôle non destructif,” Université de Lyon, 2016.

Claims

CLAIMS 1. Computer-implemented method for determining a field of an ultrasonic wave emitted by an ultrasonic transducer comprising at least one ultrasonic element towards an area of ​​interest, the transducer being separated from the area of ​​interest by several media separated by walls, the method comprising the steps of: - Receive a model (105,213) of each wall in the form of a parameterized 3D surface obtained by approximation of multilevel B-splines, - Receiving a decomposition (107) of the surface (400) of each ultrasonic element into a plurality of elementary surfaces (401) each being described by a set of points comprising a central point and four points corresponding to the vertices of a square centered on the central point, - Determine an ultrasonic field transmitted by each ultrasonic element and for each point of interest in the area of ​​interest by means of the sub-steps of: • For each mode of propagation of the ultrasonic wave among several modes and for each source point of the set of points characterizing an elementary surface among the set of elementary surfaces describing each ultrasonic element, determine (111) a path between the source point and the point of interest so as to minimize the flight time of the ultrasonic wave between these two points, determine (204) the attenuation of the ultrasonic wave on this path, calculate (112) the impulse response of the ultrasonic field for each elementary surface from the paths and attenuations determined for the set of associated points, • Sum (114) the set of calculated impulse responses, • Determine (115) the ultrasonic field transmitted by the ultrasonic element from the result of the sum.

2. Method for determining a field of an ultrasonic wave according to claim 1 in which the ultrasonic transducer comprises several ultrasonic elements, the method further comprising a sum of the ultrasonic fields of each element weighted by a predefined amplitude law, each ultrasonic field being delayed via a predefined delay or phase law.

3. Method for determining a field of an ultrasonic wave according to any one of the preceding claims comprising a prior step of decomposing the surface of each ultrasonic element comprising at least: - Decompose (106) the surface into a plurality of contiguous elementary surfaces, each elementary surface being centered on a central point, - For each of the central points, form a square centered on the central point in a plane tangent to the surface, - Characterize (107) the elementary surface by the set of points including the central point and the four vertices of the square, 4. Method for determining a field of an ultrasonic wave according to any one of the preceding claims comprising a prior step of determining a model of each wall comprising at least: - Acquire (101) an image of the walls, - Extract (104) wall meshes in the form of a point cloud, - Apply (105) a multi-level B-spline approximation to the point cloud so as to obtain a smoothed 3D parameterized surface 5. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the media separated by walls are anisotropic and the step (111) of determining a path between the source point and the point of interest is carried out by taking into account the anisotropy of each medium crossed.

6. Method for determining a field of an ultrasonic wave according to claim 5 in which the step (111) of determining a path between the source point and the point of interest is carried out by means of an algorithm for minimizing (202) the sum of the segments making up the path, each segment corresponding to a medium crossed, the sum being weighted by the slowness of each medium.

7. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the propagation modes of the ultrasonic wave comprise at least one propagation mode involving reflections of the ultrasonic wave on at least certain walls and taking into account the reflection coefficient of the path originating from the central point of each elementary surface in the calculation (112) of the impulse response of the ultrasonic field.

8. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the decomposition (107) of a surface of an element into elementary surfaces is carried out using a Voronoi diagram.

9. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the calculation (112) of the impulse response of the ultrasonic field comprises at least the calculation of the amplitude of the field and the minimum and maximum limits of the flight times corresponding to the calculated paths.

10. A method of determining a field of an ultrasonic wave according to any one of the preceding claims wherein the area of ​​interest is a transcranial region and the walls are the walls of a skull.

11. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the step of determining an ultrasonic field comprises at least the step of convolving (115) the result of the sum with a reference signal associated with the ultrasonic element or of extracting (115) the amplitude and the phase of a predefined frequency component of each impulse response to obtain an ultrasonic field transmitted by said element.

12. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the step of determining an ultrasonic field transmitted by each ultrasonic element and for each point of interest of the area of ​​interest further comprises a sub-step of determining (303) the polarization of the ultrasonic wave at the point of interest as a function of the determined path and the propagation mode.

13. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the step of determining an ultrasonic field transmitted by each ultrasonic element and for each point of interest of the area of ​​interest further comprises a sub-step of determining (304) a cumulative transmission and / or reflection coefficient for the path associated with the central point of an elementary surface.

14. Method for determining a field of an ultrasonic wave according to any one of the preceding claims in which the calculation (112) of the impulse response of the ultrasonic field comprises the determination (305) of the flight times associated with the five paths associated with an elementary surface and the determination of a temporal position and a duration of a time slot defining the impulse response from said flight times.

15. Method for determining a field of an ultrasonic wave according to claim 14 in which the calculation (112) of the impulse response of the ultrasonic field further comprises the determination (306) of the solid angle formed by the four paths corresponding to the four vertices of the square centered on the central point and the determination (307) of the divergence coefficient for a brush defined by these four paths and the path corresponding to the central point.

16. Method for determining a field of an ultrasonic wave according to claim 15 in which the calculation (112) of the impulse response of the ultrasonic field is carried out by distributing over the time slot, the energy of the brush calculated from the attenuation coefficient, the cumulative transmission and / or reflection coefficient for the path associated with the central point and the divergence coefficient.

17. A computer program comprising instructions for executing the method according to any one of claims 1 to 16, when the program is executed by a processor.

18. A processor-readable recording medium on which is recorded a program comprising instructions for executing the method according to any one of claims 1 to 16, when the program is executed by a processor.

19. An ultrasonic wave field simulator comprising a processor and a memory configured to perform the steps of the method according to any one of claims 1 to 16 and a display device for displaying the calculated field in an area of ​​interest.

Images