Method for determining an ultrasonic field by a plurality of walls

Abstract

The object of this invention is to provide a semi-analytical method that allows for accurate simulation of the transmission of an ultrasound field through the skull wall using an algorithm with computation time compatible with the preoperative background. More generally, the invention includes simulating an ultrasound field through one or more walls in order to focus the beam at a point located behind the wall, or indeed to guide the beam in a predetermined direction behind the wall.

Description

Technical Field

[0001] This invention relates to the field of simulating the propagation of ultrasound waves through heterogeneous and / or anisotropic solid media with complex geometries. The invention is particularly applicable to transcranial ultrasound therapy, but also to ultrasound imaging.

[0002] More precisely, the present invention relates to a method for simulating ultrasound waves emitted by an ultrasound transducer and passing through the skull wall, with the aim of characterizing the transmitted ultrasound beam prior to its generation, particularly its focal region. Background Technology

[0003] Transcranial ultrasound therapy involves the use of focused ultrasound to treat certain diseases. Curing brain diseases remains very difficult, primarily because pharmacological agents cannot easily penetrate the brain due to the blood-brain barrier and blood-tumor barrier (BBB and BTB). Focused low-intensity pulsed ultrasound in the brain, combined with microbubble injection, significantly increases drug distribution to brain tissue, and its therapeutic effects have been established in numerous animal models and patients. This permeation of the BBB or BTB is localized and transient, provided the beam intensity is well controlled, as microbubble implosion can lead to hemorrhage. Therefore, to ensure the safety of the protocol, precise dosing of the sound field after penetrating the bone wall is necessary.

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

[0005] In the various applications envisioned above, as well as in ultrasound imaging applications, it is necessary to precisely control the propagation of ultrasound waves, as they vary according to the geometry and mechanical properties of the skull wall.

[0006] Specifically, incident ultrasound beams may be deflected, attenuated, or aberrated, or even experience a combination of two or three of these simultaneously. Therefore, the skull is considered a complex medium that can strongly interfere with ultrasound focusing in therapeutic situations or with image reconstruction in imaging situations. This inherent problem with the skull stems primarily from the significant differences in geometry and local acoustic properties between given individuals, and also from the individual variations in the same properties. Therefore, personalized modeling of the skull is crucial for simulating ultrasound propagation and improving therapeutic and imaging tools.

[0007] Traditionally, ultrasound fields emitted through the transcranial wall have been simulated using numerical methods, which are considered more accurate than semi-analytical methods, although they are more expensive in terms of execution time and memory.

[0008] Reference [1] lists the main methods for simulating the propagation of ultrasound fields through the skull wall in the literature review, which can be divided into two main categories.

[0009] The first class of methods are numerical methods, such as the finite element method, the finite difference method, or the boundary element method, which are considered reference methods. Specifically, these methods allow the simulation of all complex physical effects at work because they consist of exact numerical solutions to the ultrasonic propagation equations. However, in most cases, they remain expensive in terms of computation time and memory usage. Furthermore, they are highly restrictive in terms of computational domain. Specifically, the latter must undergo sufficient spatial sampling to meet the numerical convergence criteria and must simultaneously include the source term (i.e., the ultrasonic probe) and the region of interest. Finally, implementation often requires substantial expertise on the user's side, particularly in the use of tools to mesh the scene, including the probe and the skull. Reference [2] provides an example of a comparative study of major numerical simulation tools.

[0010] The second category involves so-called semi-analytical methods, particularly ray tracing, which can also be implemented in transcranial ultrasound simulations. These methods are faster, cheaper in terms of memory load, and have no limitations on the computational domain. In contrast, because they are based on asymptotic approximations, they do not allow for a satisfactory full simulation of the complexity of the physical effects at work during propagation. Nevertheless, recent studies have shown that these semi-analytical methods can achieve results as accurate as numerical methods.

[0011] Among semi-analytical methods, the so-called pencil-tracing method based on ray theory is fast and reliable. However, its drawback is its high sensitivity to the regularity of surfaces describing various propagation media. Specifically, discontinuities can lead to overly intense regions (caustics) or areas where no rays can propagate (Fresnel shadows). These limitations are particularly pronounced when dealing with industrial parts described by meshes or components with complex 3D geometries.

[0012] Implementing the pen-beam tracing method requires defining the ultrasonic path between each source point and the point of interest, and obeying the laws of refraction and / or reflection at each interface according to the selected propagation mode. In configurations involving industrial parts described by meshes or components with complex 3D geometries, determining these paths is no trivial task.

[0013] The pen-beam tracing method taught in existing techniques involves blindly searching for these paths by tracing them in all directions of space from the point of interest, aiming to retain only the paths reaching the probe's emitting surface. This method is very expensive in terms of computation time, and its complexity increases with the number of interfaces to be managed. It cannot guarantee finding every ultrasound path in an exhaustive manner. The calculation of the area irradiated by each pen beam at the source point remains approximate and requires complex matrix calculations. Finally, it cannot guarantee regular sampling of the emitting region, thus potentially resulting in shadowed or overlapping areas. All these limitations lead to uncertainties in the calculation of the ultrasound field. Summary of the Invention

[0014] The object of this invention is to provide a semi-analytical method that allows for accurate simulation of the transmission of an ultrasound field through the skull wall using an algorithm with computation time compatible with the preoperative background. More generally, the invention includes simulating an ultrasound field through one or more walls in order to focus the beam at a point located behind the wall, or indeed to guide the beam in a predetermined direction behind the wall.

[0015] This invention is specifically based on different implementations of the pen-beam tracking method, which allows for significantly improved reliability of results and significantly reduced computation time involving the configuration of complex 3D surfaces. Based on spatial sampling of each ultrasonic element and regularization and parameterization of the 3D surface, it allows for the implementation of an algorithm that calculates the ultrasonic path between each source point and each point of interest using a time-of-flight minimization method.

[0016] One aspect of the present invention is a computer-implemented method for determining an ultrasonic field emitted toward a region of interest by an ultrasonic transducer comprising at least one ultrasonic element, the transducer being separated from the region of interest by a plurality of media separated by walls, the method comprising the following steps: - Receive a model for each wall, in the form of a parametric 3D surface obtained through multi-level B-spline approximation. - The surface of each ultrasonic element is decomposed into multiple basic regions, each described by a set of points including a center point and four points corresponding to the vertices of a square centered on that center point. - The ultrasonic field transmitted by each ultrasonic element and the ultrasonic field for each point of interest in the region of interest are determined by the following sub-steps: For each propagation mode of ultrasound in multiple modes, and for each source point in a set of points representing a fundamental region describing each ultrasound element, determine the path between the source point and the point of interest to minimize the ultrasound flight time between these two points, determine the attenuation of the ultrasound along the path, and calculate the ultrasound field impulse response for each fundamental region based on the path and attenuation determined for the associated set of points. Sum all calculated impulse responses. • Determine the ultrasonic field transmitted by the ultrasonic element based on the summation result.

[0017] According to a particular aspect of the invention, an ultrasonic transducer includes a plurality of ultrasonic elements, and the method further includes summing the ultrasonic fields of each element and weighting them by a predefined amplitude law, wherein each ultrasonic field is delayed by a predefined delay law or a predefined phase law.

[0018] In a variant embodiment, the method according to the invention includes a prior step of decomposing the surface of each ultrasonic element, comprising at least: - The surface is decomposed into multiple adjacent basic regions, each centered at a central point. - For each center point, a square centered at that center point is formed in a plane tangent to the surface. - The basic region is characterized by a set of points including the center point and the four vertices of the square.

[0019] In a variant embodiment, the method according to the invention includes a prior step of determining a model of each wall, comprising at least: - Obtain an image of the wall. - Extract the mesh of the walls in the form of a point cloud. - Apply multi-level B-spline approximation to the point cloud to obtain a smooth 3D parametric surface.

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

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

[0022] According to a particular aspect of the invention, the propagation mode of the ultrasound includes at least one propagation mode involving the reflection of the ultrasound from at least some walls, and the reflection coefficient of the path originating from the center point of each basic region is taken into account in the calculation of the ultrasound field impulse response.

[0023] According to a particular aspect of the invention, the surface of an element is decomposed into basic regions using a Voronoi diagram.

[0024] According to a particular aspect of the invention, the calculation of the ultrasonic field impulse response includes at least the calculation of the field amplitude and the calculation of minimum and maximum limits of the flight time corresponding to the calculation path.

[0025] According to a particular aspect of the invention, the region of interest is a transcranial region, and the wall is the wall of the skull.

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

[0027] Another subject of the invention is a computer program that includes instructions for performing the method according to the invention when the program is executed by a processor.

[0028] Another subject of the present invention is a processor-readable recording medium having a program containing instructions for performing the method according to the present invention when the program is executed by a processor.

[0029] Another subject of the invention is an ultrasonic field simulator comprising a processor and a memory configured to perform the method steps according to the invention, and a display device for displaying the calculated field in the region of interest. Attached Figure Description

[0030] Other features and advantages of the invention will become clearer from the following description with reference to the accompanying drawings.

[0031] Figure 1 A flowchart is shown, detailing the implementation steps of a method for determining an ultrasonic field emitted through at least one wall according to an embodiment of the present invention.

[0032] Figure 2 A flowchart is shown illustrating the steps for calculating the ultrasound path and attenuation according to an embodiment of the present invention.

[0033] Figure 3 A flowchart is shown illustrating the steps for calculating the impulse response of transmitted ultrasound according to an embodiment of the present invention.

[0034] Figure 4 A schematic diagram is shown illustrating an example of breaking down an ultrasound probe into multiple components.

[0035] Figure 5 A diagram is shown, illustrating... Figure 2 The principle of path calculation method,

[0036] Figure 6 The diagram illustrates the principle of calculating the ultrasonic field using the pen beam tracing method.

[0037] Figure 7 An example of ultrasonic field impulse response is shown. Detailed Implementation

[0038] Figure 1 The main steps of a method for determining the emitted ultrasound field from a source point on the surface of an ultrasound probe to a target point located behind one or more walls are illustrated in flowchart form. In the remainder of the description, the invention is described in the context of a skull wall, i.e., an application to transcranial imaging or therapy; however, the invention is not limited to this type of wall and can be applied to walls of other types of objects, particularly in non-destructive testing applications.

[0039] The method includes a preprocessing stage 100, followed by a field calculation stage 110.

[0040] The preprocessing stage 100 receives a skull model 101 as input, which is obtained, for example, by medical X-ray imaging technology, and an ultrasound probe geometry model 102 designed to generate an ultrasound field.

[0041] In step 103, the acoustic properties of the skull are extracted from the obtained model. This step is performed, for example, by the extraction technique described in either of references [3] and [4].

[0042] In step 104, a mesh is extracted from the walls of the skull, for example, using a decomposition technique for breaking down a surface into a mesh. The mesh is, for example, a 3D point cloud.

[0043] The elements describing the skull obtained in steps 103 and 104 are used to model the outer surface forming the skull wall in step 105.

[0044] In order to model the various walls of the skull by smoothing the surface, a B-spline approximation method is used in step 106, such as the method based on the technique described in reference [5].

[0045] One advantage of this method is that it consists of approximations rather than interpolation, which makes it robust to any point cloud anomalies obtained during extraction 104. This method enables the acquisition of smoother surfaces than interpolation-based techniques, thereby improving the convergence of the algorithm described below for calculating the ultrasonic path. It also has the advantage of reducing noise in the point cloud.

[0046] This approximation method is, for example, the so-called "multilevel B-spline" or "MBA" method. At the end of step 106, each wall of the skull is modeled as a 3D surface considered as the interface between two media indexed i and i+1, and this surface is defined by a function... Description, in which These correspond to the parametric coordinates of the surface described by the spline. In the case of the MBA method, .

[0047] There must be at least two interfaces defining the skull wall, but there can be more than two.

[0048] The preprocessing stage also includes a stage of modeling the ultrasonic probe. The probe used can be a single-element or multi-element ultrasonic transducer. This invention is applicable regardless of the 2D or 3D surface geometry of each transducer element.

[0049] In step 106, spatial sampling is performed on each element of the probe's emitting surface to obtain a set of basic regions characterized by a center point and area, and uniformly distributed across the entire surface of the elements.

[0050] The number of basic regions is an input parameter to this method. This sampling can be performed, for example, using a Voronoi diagram. The sum of the basic regions equals the total area of ​​the component surface.

[0051] Figure 4 An example of a circular transducer element 400 is illustrated. The element 400 is decomposed into multiple basic regions 401 centered at C.

[0052] In step 107, for each basic region, the coordinates of the four vertices of a square centered at point C, which is the center of the basic region, are determined. The resulting virtual square has the same area A as element 401. In the case of a three-dimensional surface, the virtual square is formed in a directly orthogonal coordinate system (…). , , ) plane ( , In ) It is a vector perpendicular to the basic region 401.

[0053] This step is in Figure 4 The right side of the diagram provides an explanation. For the example of component 401, the diagram shows a virtual square centered at C, with its four vertices in a directly orthogonal coordinate system (…). , , The coordinates in the coordinate system are defined by the following:

[0054]

[0055]

[0056]

[0057]

[0058] Therefore, these four points are defined for each basic region to serve as source points in the path calculation prior to the field calculation. This modeling of surface elements allows for the avoidance of computationally expensive and complex matrix calculations employed in existing pen-beam tracing methods to determine the divergence coefficients and durations of the time-dependent rectangular functions used when calculating the impulse response associated with the surface element. Furthermore, describing each surface element as a square region allows for handling any possible probe geometry.

[0059] At the end of the preprocessing stage 100, a model of the skull wall was thus obtained, on the one hand by a three-dimensional parametric function, and on the other hand by characterizing the coordinates of five points (the center and four vertices of the virtual square) of each basic region generated by the decomposition of the ultrasonic transducer elements.

[0060] All these parameters are passed as input to the second stage 110 for calculating the ultrasonic field. Figure 1 The specific example given is for calculating the ultrasonic field corresponding to one ultrasonic element of the transducer. If the transducer comprises multiple ultrasonic elements, an additional combination step must be applied to the ultrasonic field calculated for each element. This step will be described in more detail below.

[0061] The method for calculating the ultrasound field, which is the subject of Phase 2 (110), comprises three iterative loops. The first loop aims to scan all points in the computational region, i.e., all points in the brain in the case of a transcranial application, or all points in the region of the object to be examined in the case of a non-destructive testing application. For each target point in the computational region, the ultrasound field emitted through the wall modeled in step 105 is calculated.

[0062] The calculation first requires calculating path 111 by minimizing the time of flight, and then calculating the impulse response 112. These two steps are iterated in a second iteration loop that scans multiple propagation modes of the ultrasound waves, and in a third iteration loop that scans each point in each basic region characterizing each element of the probe. These points correspond to the points determined in step 107.

[0063] Figure 2 The calculation of the ultrasound path is explained in detail 111. Figure 3 The calculation of the associated impulse response is explained in detail 112.

[0064] The purpose of ultrasound path calculation 111 is to find the shortest path between point A on the probe and target point B on the other side of the skull wall, which is modeled by multiple surfaces.

[0065] Figure 5The principle is illustrated using an example of a wall modeled by two surfaces, f1 and f2, which separate three media characterized by ultrasonic propagation speeds of c1, c2, and c3, respectively. The purpose of calculation step 111 is to find the intersection points I1 and I2 of the ultrasonic path transmitted from A to B, thereby obtaining the shortest flight time.

[0066] If I1 = (x1, y1, z1) represents interface point 1 through which the ray passes, and I2 = (x2, y2, z2) represents interface point 2, then and This is because the interface is described by splines.

[0067] Therefore, the four parameters that need to be changed are:

[0068]

[0069] The objective of path calculation 111 is to minimize the flight time given by the following relation: .

[0070] The path calculation step 111 receives a source point 210 (i.e., a point on the surface of the transmitting probe) and a calculation point 211 (i.e., a point in the calculation region) as inputs. It also receives the interface 213 of the wall to be modeled and the velocity 212 associated with each medium separated by the interface as inputs.

[0071] Step 201 includes the source point With calculation point Ray tracing is performed between the ray points to define the intersections of the ray with each interface. The intersections belong to an approximate surface obtained using splines. Step 201 can be performed using a ray tracing algorithm, such as the one described in the applicant's patent application FR3114032.

[0072] Step 202 then involves finding the intersection coordinates that minimize the flight time given by the aforementioned relation.

[0073] This optimization step can be performed using any suitable optimization algorithm, such as a quasi-Newton minimization algorithm like the L-BFGS (Limited Memory Broyden-Fletcher-Goldfarb-Shannon) algorithm, or using the method described in patent application FR3114032.

[0074] The overall goal of optimization step 202 is to find the source point. With calculation point Passing between The path of each interface makes the time-of-flight function Minimum. Due to the slowness of each medium. The wave is homogeneous, and its flight time is the sum of the flight times in each medium it passes through. To find this minimum, we need to find the intersection of the ultrasonic path with each interface. Therefore, a minimization algorithm such as the L-BFGS algorithm is used to solve the following problem: in:

[0075] These are points on a parameterized surface describing the interface between media with indices i and i+1.

[0076] For anisotropic media, the slowness s i Depending on the direction of propagation, the above relationship can be written as:

[0078] Considering anisotropic media in ultrasound path calculations allows for a better characterization of the objects through which the wave passes, especially the human skull, which has unique characteristics for each individual.

[0079] For the anisotropic medium under consideration, the ultrasonic path time of flight is characterized. The cost function to be minimized must be parameterized to account for the dependence of slowness on the propagation direction. This parameterization can be achieved using an analytical function s that describes the variation of slowness with the propagation direction. i (In 3D representation) This parameterization is performed and must be evaluated in each iteration. For optimization purposes, this parameterization can also be obtained using a tabulation established in the initial steps to show how the slowness varies with the direction of propagation (in 3D representation).

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

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

[0082] To calculate the total attenuation along the path, the attenuation of each segment of the path through each medium is summed. If each medium... In frequency The lower has an attenuation coefficient and decay index The total attenuation is then obtained through the following relationship:

[0083]

[0084] The attenuation value for each medium is the input data for this method.

[0085] In one embodiment of the invention, path calculation 111 is performed for various propagation modes. In particular, some propagation modes involve reflections from the internal interfaces of heterogeneous structures (e.g., skulls) traversed.

[0086] In this case, reflection is considered as an additional path segment. The same method is then applied, but with an additional interface added as many times as the number of reflections, corresponding to the multiple propagations of the wave through a given medium.

[0087] According to other propagation embodiments, the velocity and / or attenuation value in the medium traversed may vary from mode to mode.

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

[0089] More precisely, calculate the path for each center of the basic region of the probe element, and the four vertices of the virtual square formed around that center (as explained with respect to step 107).

[0090] Figure 3 The implementation of step 112 for calculating the impulse response is explained in detail.

[0091] This calculation step receives the path originating from the center point 301 and the path 302 originating from the four vertices of the virtual square as inputs for each basic region of the probe.

[0092] The calculation of the impulse response 112 is based on the so-called pen beam tracking method described in reference [6] (referred to as "méthode des pinceaux" in the reference).

[0093] First, the principle of the pen beam tracing method is described.

[0094] To reconstruct the emitted wavefront and thus calculate the Rayleigh integral associated with the ultrasonic field, the pen-beam tracing method can be used.

[0095] Physically, the beam represents the acoustic energy flow originating from a 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 flow through the surface of the beam is zero. Therefore, the energy per unit area through surface S is inversely proportional to its area.

[0096] The central ray of the pen bundle follows a Fermat minimum time path, obtained in path calculation step 111. The paraxial rays follow a different direction than the central ray, but are considered close enough for the paraxial approximation to be valid (since the paraxial approximation is a first-order approximation, the error increases with distance from the axis). The paraxial rays are characterized by four differential components represented in a coordinate system in a plane perpendicular to the central axis of the pen bundle, thus allowing the paraxial rays of the pen bundle to be defined by a vector containing four coordinates:

[0097] dx and dy characterize the intersection point of the paraxial ray and the pencil beam surface S, while dSx and dSy give the projection of the slowness vector of the paraxial ray onto the surface S.

[0098] Under the paraxial approximation, the variation of the beam vector in a homogeneous medium can be characterized by a 4×4 matrix that depends on the propagation medium. A characteristic 4×4 matrix of the interface between two materials can also be defined. These matrices are called propagation matrices, and they are typically written in block notation as follows: , where matrices A, B, C and D are 2×2 matrices.

[0099] Considering that the pen beam passes through a series of uniform media and interfaces characterized by the propagation matrix (Li) (for 1≤i≤n), the resulting pen beam can be written in the following form: .

[0100] Since heterogeneous media can be modeled by a series of homogeneous media and the interfaces separating them, these matrices are sufficient to describe the changes in the pen beam as it passes through the heterogeneous media.

[0101] Using a pen beam in ultrasonic field calculations allows for the simulation of the spatial propagation of ultrasonic waves. To obtain the final amplitude of the pen beam, its divergence must be considered. Specifically, the law of conservation of energy within the pen beam means that the energy flow decreases as the pen beam aperture increases, and in the context of a small-aperture pen beam approximation, the energy flow is uniformly distributed on the substrate surface S.

[0102] Therefore, the ratio between the final intensity IS per unit area and the intensity IΩ per unit initial solid angle can be written as the ratio between the initial solid angle dΩ and the final base area dS of the pen beam: .

[0103] In addition to the amplitude factor related to the divergence of the beam during its propagation, each interface traversed reduces the beam's energy because various modes are generated at each interface. The energy distribution among the various propagation modes (transmission, reflection, and dissipation modes for each polarization type) is obtained from the continuity condition at the interface and quantified by calculating the Fresnel coefficients.

[0104] To obtain the value of the amplitude transmission coefficient (which is directly applicable to the displacement value and therefore to the amplitude of the ultrasonic field), the ratio between the beam sections before and after the interaction needs to be considered.

[0105] For each pulse emitted from the considered field point and whose final surface intersects the ultrasonic transducer surface, a fundamental impulse response can be extracted. The summation of these fundamental responses yields the total impulse response, which is then used to determine the amplitude value of the field as a function of time. The fundamental impulse response of the pulse can be compared to a time-dependent step function, whose amplitude, time span, and center time are extracted from the final characteristics of the pulse. The amplitude value corresponds to energy. Therefore, it is distributed over the time span of the pulse, which is equivalent to adding a factor to the amplitude value as the reciprocal of the time span.

[0106] Therefore, the magnitude A is given by the following relationship:

[0107] Where A0 is the initial amplitude, DF is the divergence factor, and T... A dS is the amplitude transmission coefficient, dS is the area of ​​the pen beam at arrival, and Δt is the time span.

[0108] The time span Δt is calculated as the maximum step difference evaluated on the pen beam, which is a function of the incident angle of the pen beam on the transducer and its spatial range.

[0109] The center time is equal to the flight time of the central ray of the pen beam.

[0110] For multi-element transducers using delay or amplitude laws, the center time and amplitude of each beam are modified according to the intersecting transducer elements.

[0111] Since the Fresnel coefficients can be complex, the magnitude of the displacement itself can also be complex. Furthermore, since the displacement is a vector field, the calculated response is a signal with three complex dimensions, or six real dimensions.

[0112] Once the fundamental impulse responses at a given field point are calculated, they are summed into a single impulse response. The latter corresponds to the field radiated by the transducer to the field point when each element of the transducer emits a time Dirac pulse.

[0113] Therefore, given the input signal S0 experienced by the ultrasonic transducer, the ultrasonic field at the field point is obtained by convolving the complex impulse response signal with the input signal. For each component of displacement x, y, and z, the real part is convolved with the reference signal, while the complex part is convolved with the Hilbert transform of the reference signal. Therefore, the calculation of each scalar component si(t) of the field (the impulse response denoted as Rii(t)) is written as: si(t) = Rii(t) S0(t).

[0114] The paraxial ray model is based on the propagation matrix, which allows for a linear approximation of the paraxial ray around the central ray of the ray.

[0115] This invention proposes using interpolation to approximate the paraxial rays contained within the pen bundle. To this end, in addition to the propagating central ray rc, the propagation of four rays r0, r1, r2, and r3 defining the pen bundle is simulated, such as... Figure 6 As shown above, the four rays r0, r1, r2, and r3 correspond to four paths originating from the four vertices of a virtual square centered on the basic region obtained by decomposing the ultrasonic transducer elements.

[0116] This method allows for avoiding the computation of the propagation matrix along the pen beam path, while making it feasible to compute the amount of the basic impulse response associated with the pen beam through integrals.

[0117] This method enables the sampling of the emission surface region of each element in a controlled, regular, and uniform manner. The location and area covered by each beam are precise and controlled, which is not the case when ray tracing starts from a calculation point and beam parameters are estimated through matrix calculations. Therefore, the beam tracing method takes into account the actual emission area of ​​each element.

[0118] Determine the time window of the rectangular function for modeling the pen beam so as to include the arrival times of all component rays of the pen beam. Therefore, the bound t of the time-dependent rectangular function representing the pen beam impulse response is... min and t max It is defined as the lowest and highest times in the calculated set of flight times.

[0119] Figure 7 The obtained rectangular function changes over time.

[0120] Using the four-ray pen beam model, the pen beam area is the area bounded by the points defined by the four rays, such as... Figure 6 As shown. The initial solid angle is defined by the area of ​​the ray beam at a unit distance.

[0121] The divergence coefficient can be calculated from these quantities and the initial slowness S0 of the axial ray using the following relationship: .

[0122] In this way, the cumulative amplitude transmission coefficient T of the central ray is taken into account. A And the initial amplitude A0 can be obtained in the same way as before through the relation. Distribute it over the time span of the pen bundle to calculate the displacement amplitude at the field point under discussion.

[0123] This simplified pen beam model allows for the reproduction of the fundamental pen beam impulse response without calculating the paraxial matrix, thus avoiding costly matrix computations (especially in anisotropic media), and is replaced by ultrasonic ray tracing operations. For an isolated pen beam, this method requires a total of 5 ultrasonic rays to generate the fundamental impulse response.

[0124] Therefore, step 112, which calculates the impulse response, is performed based on the above-described pen beam tracing method.

[0125] In step 303, polarization calculations are performed based on the central path 301.

[0126] In step 304, the cumulative transmission or reflection coefficient TA of the center path is calculated (depending on the propagation mode). In cases where the selected propagation mode involves both crossing various interfaces and reflecting from various interfaces, the cumulative coefficient is calculated based on the transmission and reflection coefficients of each interface in question.

[0127] In step 305, the flight times corresponding to the five paths 301 and 302 received as input are calculated. From these flight times, the time position and duration of the time-dependent rectangular function describing the impulse response are derived.

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

[0129] In step 307, the solid angle is determined using the relational formula. Calculate the divergence coefficient DF.

[0130] In step 308, t is calculated in the manner described above. min and t max The value of .

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

[0132] The impulse response of point pairs (source point, computation point) and propagation modes is given by the bounds, amplitude, and polarization of the time-dependent rectangular function.

[0133] More precisely, the ultrasonic field impulse response of the fundamental region and propagation mode under discussion is calculated by distributing the beam energy on the obtained time-dependent rectangular function using the attenuation coefficient, reflection coefficient, and / or transmission coefficient of the central ray, as well as the divergence coefficient.

[0134] Now refer to it again Figure 1The calculation continues in step 113 to determine the total magnitude of the signal, which is obtained by determining the minimum and maximum times among all time-dependent rectangular functions 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, all calculated impulse responses are summed.

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

[0137] Finally, in step 116, the pressure field at each point in the calculation region for an ultrasonic element is obtained.

[0138] When the transducer is multi-element, the algorithm 110 used to calculate the ultrasonic field must be iterated for each element of the transducer. Additional steps must also be applied ( Figure 1 (Not shown in the image) in order to combine the various fields emitted by each element into a total field associated with the transducer. For this purpose, predefined amplitude and delay laws are considered in this combination, applied according to the target. More precisely, this additional step involves calculating the sum of the ultrasonic fields of each element, weighting the sum according to the amplitude law, and delaying each field according to the delay law. In other words, the amplitude and delay laws respectively give the amplitude and delay to be applied to each field corresponding to each element.

[0139] Compared to existing methods, this invention offers numerous advantages. It is capable of handling ultrasound probes with surfaces of any geometry, a particular advantage in the context of transcranial applications where the probe's emitting surface may have a spherical profile, or more generally be shaped in 3D to match the shape of the skull. Similarly, in general, each element has a specific projected surface with a complex geometry, such as a polygonal type. Specifically, each element of any geometry is decomposed into basic regions and combined to form four points of interest associated with the element's center point, enabling adaptation to any geometry.

[0140] Furthermore, this invention allows for the consideration of anisotropic media in path calculations and the attenuation of each medium in ultrasonic field calculations. Finally, better modeling of the surface defining the skull also allows for improved accuracy and convergence time in ultrasonic path calculations.

[0141] The present invention can be implemented as a computer program including instructions for its execution, so as to allow simulation of the ultrasound field generated for each point in a target region. The computer program can be recorded on a processor-readable recording medium. The simulator may include a processor and memory for performing the method steps for generating an ultrasound field according to the present invention, and a display device for displaying the field generated for a region of interest (e.g., a transcranial region).

[0142] This simulation tool allows for assistance to the operator before the ultrasonic field is generated by the transducer, enabling better calibration of the field before its generation, or, in the case of a multi-element probe, calculating the phase or delay law that allows correction of aberrations experienced by the field as it passes through the skull 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 compressive 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. [4] L. Marsac et al., "Ex vivo optimization of a heterogeneous speed of sound model of the human skull for non-invasive transcranial focusedultrasound at 1 MHz", Int. J. Hyperth., vol. 33, no. 6, pp. 635–645, 2017, doi: 10.1080 / 02656736.2017.1295322. [5]BG Lee, JJ Lee, and J. Yoo, "An efficient scattered dataapproximation 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. [6] H. Chouh, “Simulations interactives de champ ultrasonore pour desconfigurations complexes de contrôle non destructif”, University of Lyon, 2016.

Claims

1. A computer-implemented method for determining an ultrasonic field emitted toward a region of interest by an ultrasonic transducer comprising at least one ultrasonic element, the transducer being separated from the region of interest by a plurality of media separated by walls, the method comprising the steps of: - Receive the model (105, 213) for each wall, which is in the form of a parametric 3D surface obtained by multi-level B-spline approximation. - The surface (400) of each ultrasonic element is decomposed (107) into multiple basic regions (401), each basic region being described by a set of points including a center point and four points corresponding to the vertices of a square centered on that center point. - The ultrasonic field transmitted by each ultrasonic element and the ultrasonic field for each point of interest in the region of interest are determined by the following sub-steps: For each propagation mode of ultrasound in multiple modes, and for each source point in a set of points representing a basic region describing each ultrasound element, determine (111) the path between the source point and the point of interest to minimize the ultrasound flight time between the two points, determine (204) the attenuation of the ultrasound along the path, and calculate (112) the ultrasound field impulse response for each basic region based on the path and attenuation determined for the associated set of points. Sum all calculated impulse responses (114). • The ultrasonic field transmitted by the ultrasonic element is determined based on the summation result (115).

2. The method for determining an ultrasonic field as claimed in claim 1, wherein the ultrasonic transducer comprises a plurality of ultrasonic elements, the method further comprising summing the ultrasonic fields of each element and weighting them by a predefined amplitude law, wherein each ultrasonic field is delayed by a predefined delay law or a predefined phase law.

3. The method for determining an ultrasonic field as described in any of the preceding claims, comprising the prior step of decomposing the surface of each ultrasonic element, including at least: - The surface is decomposed (106) into multiple adjacent basic regions, each with a center point as its center. - For each center point, a square centered at that center point is formed in a plane tangent to the surface. - The basic region (107) is characterized by a set of points including the center point and the four vertices of the square.

4. The method for determining an ultrasonic field as described in any of the preceding claims, comprising the prior step of determining a model of each wall, including at least: - Obtain the image of the (101) wall, - Extract the mesh of the walls (104) in the form of a point cloud. - Apply (105) multi-level B-spline approximation to the point cloud to obtain a smooth 3D parameterized surface.

5. The method for determining an ultrasonic field as described in any of the preceding claims, wherein the medium separated by the wall is anisotropic, and the step (111) of determining the path between the source point and the point of interest is performed taking into account the anisotropy of each medium traversed.

6. The method for determining an ultrasonic field as claimed in claim 5, wherein the step (111) of determining the path between the source point and the point of interest is performed by an algorithm (202) for minimizing the sum of the component segments of the path, each segment corresponding to a medium traversed, the sum being weighted by the slowness of each medium.

7. The method for determining an ultrasonic field as claimed in any of the preceding claims, wherein the propagation mode of the ultrasonic wave includes at least one propagation mode involving the reflection of the ultrasonic wave from at least some walls, and the reflection coefficient of the path originating from the center point of each basic region is taken into account in the calculation (112) of the ultrasonic field impulse response.

8. The method for determining an ultrasonic field as described in any of the preceding claims, wherein the surface of the element is decomposed (107) into basic regions using a Voronoi diagram.

9. The method for determining an ultrasonic field as described in any of the preceding claims, wherein the calculation of the ultrasonic field impulse response (112) includes at least the calculation of the field amplitude and the calculation of minimum and maximum limits of the flight time corresponding to the calculation path.

10. The method for determining an ultrasound field as described in any of the preceding claims, wherein the region of interest is a transcranial region and the wall is the wall of the skull.

11. The method for determining an ultrasonic field as described in any of the preceding claims, wherein the step of determining the ultrasonic field comprises at least the following steps: The summation result is convolved with a reference signal associated with the ultrasonic element (115), or the amplitude and phase of a predefined frequency component of each impulse response are extracted (115) to obtain the ultrasonic field transmitted by the element.

12. The method for determining an ultrasonic field as claimed in any of the preceding claims, wherein the step of determining the ultrasonic field transmitted by each ultrasonic element and the ultrasonic field for each point of interest in the region of interest further comprises the sub-step of determining (303) the polarization of the ultrasonic wave at the point of interest based on the determined path and propagation mode.

13. The method for determining an ultrasonic field as claimed in any of the preceding claims, wherein the step of determining the ultrasonic field transmitted by each ultrasonic element and the ultrasonic field for each point of interest in the region of interest further comprises the sub-step of determining (304) cumulative reflection and / or transmission coefficients for the path associated with the center point of the basic region.

14. The method for determining an ultrasonic field as claimed in any of the preceding claims, wherein calculating (112) the ultrasonic field impulse response comprises determining (305) the flight times associated with five paths associated with a base region, and determining the time position and duration of a time-dependent rectangular function defining the impulse response based on the flight times.

15. The method for determining an ultrasonic field as claimed in claim 14, wherein calculating (112) the ultrasonic field impulse response further includes determining (306) a solid angle formed by four paths corresponding to the four vertices of a square centered at the center point, and determining (307) a divergence coefficient for a beam defined by the four paths and the path corresponding to the center point.

16. The method for determining an ultrasonic field as claimed in claim 15, wherein the ultrasonic field impulse response (112) is calculated by distributing the beam energy over a time-dependent rectangular function, the energy being calculated based on the attenuation coefficient, the cumulative reflection and / or transmission coefficient of the path associated with the center point, and the divergence coefficient.

17. A computer program comprising instructions for performing the method as claimed in any one of claims 1 to 16 when executed by a processor.

18. A processor-readable recording medium having a program recorded thereon including instructions for performing the method as described in any one of claims 1 to 16 when the program is executed by a processor.

19. An ultrasonic field simulator comprising a processor and a memory configured to perform the method steps as described in any one of claims 1 to 16, and a display device for displaying a calculated field in a region of interest.

Images