Method for determining the field strength of an ultrasonic wave through several walls

Abstract

The object of the present invention is to propose a semi-analytical method for accurately simulating the transmission of an ultrasound field through a cranial wall using algorithms with compatible computation 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 the wall. Figure 1

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 simulation of 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] It relates more specifically to a method of simulating an ultrasonic wave emitted by an ultrasonic transducer through the cranial wall with a view to controlling the transmitted ultrasonic beam, and in particular the focal zone, before generating it.

[0003] Cerebral or transcranial ultrasound therapy consists of the use of focused ultrasound to treat certain diseases. The cure of brain diseases remains very difficult, mainly due to the poor access of pharmacological agents to the brain due to the blood-brain and blood-tumor barriers (BBB and BHT). The focusing of low-intensity pulsed ultrasound in the brain, combined with the injection of microbubbles, significantly increases the delivery of drugs to brain tissues, with an established therapeutic effect in many animal models as well as in patients. This permeabilization of the BBB or BHT is local and transient provided that the intensity of the beam is well controlled because the implosion of the 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 interesting alternative to brain surgery and radiotherapy. HIFU therapy induces a rise in temperature 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 the surrounding healthy tissues.

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

[0006] Indeed, the 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 cranium is then considered as a complex medium that strongly disrupts the focusing of ultrasound for therapy, or the reconstruction of images 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 individuals. 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 ultrasonic fields generated towards transcranial walls are carried out 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 of the finite element, finite difference or boundary element method type, 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 sampled spatially, to respect numerical convergence criteria, and it must contain both the source term, i.e. the ultrasonic probe, and the area of ​​interest.Finally, the implementation generally requires significant expertise from the user, notably 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 calculation area. On the other hand, since they are based on asymptotic approximations, they do not allow the full complexity of the physical phenomena involved during propagation to be simulated satisfactorily. Nevertheless, studies Recent studies have shown that these semi-analytical methods make it possible to obtain results as precise as numerical methods.

[0011] The object of the present invention is to propose a semi-analytical method allowing the transmission of an ultrasound field through a cranial wall to be precisely simulated using algorithms with 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.

[0012] 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.

[0013] 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.

[0014] 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,

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

[0016] 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 medium crossed.

[0017] 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.

[0018] 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 originating from the central point of each elementary surface in the calculation of the impulse response of the ultrasonic field.

[0019] 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.

[0020] 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.

[0021] 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.

[0022] 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 for obtain an ultrasonic field transmitted by said element.

[0023] 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.

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

[0025] 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.

[0026] Other characteristics and advantages of the present invention will appear more clearly on reading the description which follows in relation to the following appended drawings.

[0027] [Fig.l] 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

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

[0029] [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,

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

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

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

[0033] [Fig.7] represents an example of an impulse response of an ultrasonic field,

[0034] [Fig.l] illustrates, on a flowchart, the main stages of a method of de termination of 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.

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

[0036] 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 ultrasonic probe 102 for use in generating an ultrasonic field.

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

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

[0039] 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.

[0040] 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].

[0041] 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 aberrant points in the point cloud obtained at extraction 104. This method makes it possible to obtain smoother surfaces than with a technique based on an interpolation, which makes it possible to improve the convergence of the ultrasonic path calculation algorithms described below. It also has the advantage of reducing noise in the point cloud.

[0042] The approximation method is for example a so-called multi-level 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 f (u. v) — (xy ^j) where (¾¼) corresponds to the parametric coordinates of the surface described by means of splines. In the case of the MBA method, we have (Uj, Vj) = (x, y).

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

[0044] The pre-processing phase also includes a modeling phase of the ultrasonic probe. The probe used may 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.

[0045] In step 106, each element of the emitting surface of the probe is spatially sampled to obtain a set of points, distributed uniformly over the entire surface of the element.

[0046] The number of points is an input parameter of the method. This sampling can be obtained for example from a Voronoi diagram. The sum of the elementary surfaces associated with each point is equal to the total surface of the element.

[0047] [Fig. 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.

[0048] 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 (u, v) of a direct orthonormal reference frame (m, v, δî) with n a vector normal to the elementary surface 401.

[0049] 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 (w, v, n) by the coordinates: £ I r1 !! IMI

[0050] Jâ s 7 — VT 2 T 2

[0051] SE^-C ? 2 113] | 2 d | v11

[0052] „ _ -, \[Â ow - c - 2 - 2

[0053] 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 the probe elements makes it possible to address any probe geometry.

[0054] At the end of the pre-processing phase 100, a model of the walls of the skull is therefore obtained 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.

[0055] All of these parameters are provided as input to a second phase 110 of calculating ultrasonic fields. The particular example of [Fig.l] 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.

[0056] The method for calculating the ultrasonic field which is the subject of the second phase 110 comprises three iteration loops. The first loop aims to go through all the points in the calculation zone, i.e. the points of the brain in the case of a trans-

[0057]

[0058]

[0059]

[0060]

[0061]

[0062]

[0063]

[0064]

[0065]

[0066]

[0067] cranial 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 area, an ultrasonic field emitted through the walls modeled in step 105 is calculated. 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. [Fig.2] details the calculation 111 of ultrasonic paths and [Fig.3] details the calculation 112 of the associated impulse response. The goal of the 111 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. This principle is illustrated in [Fig.5] for an example of a wall modeled by two surfaces fi and f2 separating three media characterized by ultrasonic wave propagation speeds respectively equal to cb c2 and c3. The aim of the calculation step 111 is to find the intersection points L and I2 of the ultrasonic path emitted from A to B which make it possible to obtain the shortest flight time. Let II = (xl,yl,zl) be the interface point 1 and 12 = (x2,y2,z2) of interface 2 through which the ray passes, we have Zj) = and (x^y^) = f2(u2, V2) because the interfaces are described by splines. The 4 parameters to vary are therefore: Wj U2 V2 The objective of the journey calculation 111 is to minimize the journey time given by the relation: The path calculation step 111 receives as input a source point 210, that is to say a point on the surface of the emitting probe and a calculation point 211, 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. Step 201 consists of a ray trace 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 using a ray tracing algorithm as described for example in the Applicant's patent application FR3114032.

[0068] Step 202 then consists of searching for the coordinates of the intersection points which minimize the flight time given by the previous relationship.

[0069] This optimization step can be carried out by any appropriate 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 FR3114032.

[0070] 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 function1 is minimal. As each medium is homogeneous in slowness s = the time-of-flight is the sum of the times-of-flight 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 (B ? L with : (U-, V.) , = sd v i) 11 + 1 'ne{LN} / 111 JD 1 f7 1 1 ~ f ^i+1) fê-f U^, ) | |

[0073] / .( uj, Vj) is the parametric surface which describes the interface between the media of index i and i+1.

[0074] For anisotropic media, the slowness s; depends on the direction of propagation and the previous relationship is then written: that the L-BFGS algorithm to solve the following problem:

[0071]

[0075] hi,., s „ (s J / B-.flu-, / \ii l| + v n) 11

[0076]

[0077] 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.

[0078] After optimization, a description 203 of the ultrasonic path between the source point and the calculation point is obtained.

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

[0080] To calculate the total attenuation on a path, we sum the attenuations on each of the portions of the path in each medium crossed. If each medium has an attenuation coefficient ai at the frequency fOi, and an attenuation exponent, the total attenuation is obtained via the following relationship:

[0081] ., ,

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

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

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

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

[0086] 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.

[0087] 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).

[0088] [Fig.3] details the implementation of step 112 of impulse response calculation.

[0089] 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.

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

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

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

[0093] Physically, a brush represents a flow of acoustic energy coming from the point energy source O and propagating towards the surface S. Since O is the only energy source considered in the modeling, 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.

[0094] The central ray of the brush follows the Fermât 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 a paraxial ray of the brush to be defined by a vector with 4 coordinates:

[0095] dx dy dSx dx and dy characterize the position of the point of intersection between the ray \dSyl

[0096]

[0097]

[0098]

[0099] paraxial and the surface S of the brush. dSx and dSy give the projection of the slowness vector of the paraxial ray on the surface S. In the framework of the paraxial approximation, the evolution of a brush vector in a homogeneous medium can be characterized by a 4x4 matrix dependent on the propagation medium. We can similarly define 4x4 matrices characteristic of interfaces between two materials. We call these matrices propagation matrices and we usually denote them with a block notation '■ r _ ! AB \ where the matrices \cd) A, B, C and D are matrices of size 2x2. Considering a brush crossing a sequence of homogeneous media and interfaces characterized by the propagation matrices (L;) for l <i<n, le pinceau résultant peut s’écrire sous la forme : Since a heterogeneous medium can be modeled by a sequence of homogeneous media and interfaces separating them, these matrices are sufficient to describe the evolution of a brush when crossing a heterogeneous medium. 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.

[0100] The ratio between the final surface intensity Is and the intensity per initial solid angle unit IQ is therefore written as the ratio between the initial solid angle dQ and the final base surface of the brush dS: h _ dQ Kl dS

[0101] 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 distribution of energy 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.

[0102] 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.

[0103] For each brush emitted from the field point considered and whose final surface intersects the surface of the ultrasonic translator, 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.

[0104] The amplitude A is therefore given by the relation 4 _ AnxDFxTAxdS " At

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

[0106] 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.

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

[0108] For a multi-element translator using delay or amplitude laws, the central time and amplitude of each brush are modified in coherence with the element of the intersected translator.

[0109] Since the Fresnel coefficients can be complex, the amplitude of the displacement is itself potentially complex. Moreover, since the displacement is a vector field, the calculated response is a signal with 3 complex dimensions, or 6 real dimensions.

[0110] Once the elementary impulse responses have been 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.

[0111] Therefore, given an input signal So 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 component x, y and z 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)

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

[0113] An approximation of the rays included in the brush can also be obtained using interpolation. To do this, we simulate the propagation of the four rays r0, rb r2 and r3 delimiting the brush in addition to propagating the central ray rc, as shown in [Fig.6]. As explained previously, the four rays r0, rb r2 and r3 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.

[0114] 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.

[0115] 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 tmin and tmax 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 flight times calculated.

[0116] [Fig.7] represents the time slot obtained.

[0117] Using the four-ray brush model, the final brush surface is the area bounded by the end points of the four rays as shown in [Fig.6]. The initial solid angle is defined by the surface of the ray brush at a unit distance.

[0118] The divergence coefficient can be calculated from these quantities and the initial slowness So of the axial ray via the following relation: _ 5o V «12

[0119] In this way, considering the cumulative amplitude transmission coefficient TA of the central ray and the initial amplitude Ao, we can calculate in the same way as pre the amplitude of the displacement at the field point considered by distributing it on the temporal spread of the brush via the relation _ AaxDFxTAxdS. HAS /

[0120] This simplified brush model allows reproducing an elementary brush impulse response without having to calculate paraxial matrices, removing the need for expensive matrix calculations (particularly 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.

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

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

[0123] In step 304, a cumulative transmission or reflection coefficient TA is calculated. (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.

[0124] In step 305, the flight times corresponding to the 5 paths 301, 302 received as input are calculated.

[0125] In step 306, the solid angle formed by the four paths r0, rb r2 and r3 corresponding to the vertices of the virtual square, at the calculation point, is calculated.

[0126] In step 307, the divergence coefficient DF is calculated from the solid angle via the relation jjp _ iS'o -r— /

[0127] In step 308, the values ​​of tmin and tmax are calculated in the manner described previously.

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

[0129] 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.

[0130] We then return to [Fig. 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.

[0131] In step 114, all of the calculated impulse responses are summed.

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

[0133] Finally, in step 116, the pressure field is obtained at any point in the calculation zone for an ultrasonic element.

[0134] 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 [Fig.l]) 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.

[0135] The invention has several advantages over the methods of the prior art. It is suitable for taking into account ultrasonic probes having surfaces of any geometry, which has a particular advantage in the context of transcranial applications for which the emitting surface of the probe may 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 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.

[0136] Furthermore, the invention also makes it possible to take into account anisotropic media in the calculation of the paths and to take into account the attenuation of each media 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 the ultrasonic paths.

[0137] 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.

[0138] Such a simulation tool allows assistance to an operator prior to the generation of the ultrasonic field by means of a transducer so as to allow better calibration of the field before its generation or, in the case of a multi-element probe, to calculate the phase or delay laws, making it possible to correct the aberrations undergone by the field when crossing the cranial wall. References

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

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

[0141] [3] J.-F. Aubry, M. Tanter, M. Pernot, J.-L. Thomas, and M. Fink, "Experimental dé monstration of noninvasive transskull adaptive focusing based on prior computed to-mography scans,” J. Acoust. Soc. Am., vol. 113, no. 1, pp. 84-93, 2003, doi: 10.1121 / 1.1529663.

[0142] [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.

[0143] [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.

[0144] [6] H. Chouh, "Simulations interactives de champ ultrasonore pour des confi gurations complexes de contrôle non destructif,” Université de Lyon, 2016.

Claims

Claims

1. A 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: - Receiving a model (105,213) of each wall in the form of a parameterized 3D surface obtained by multi-level B-spline approximation, - 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 of 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. A method of determining a field of an ultrasonic wave according to claim 1 wherein 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: - Decomposing (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, forming a square centered on the central point in a plane tangent to the surface, - Characterizing (107) the elementary surface by the set of points comprising 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: - Acquiring (101) an image of the walls, - Extracting (104) meshes of the walls in the form of a point cloud, - Applying (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. A method of determining a field of an ultrasonic wave according to claim 5 wherein the step (111) of determining a path between the source point and the point of interest is performed by means of an algorithm minimization (202) of 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 preceding claim wherein the area of interest is a transcranial region and the walls are the walls of a skull.

11. A method of determining a field of an ultrasonic wave according to any one of the preceding claims wherein 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. A computer program comprising instructions for executing the method according to any one of claims 1 to 11, when the program is executed by a processor.

13. A processor-readable recording medium having recorded thereon a program comprising instructions for executing the method according to any one of claims 1 to 11, when the program is executed by a processor.

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

Images