FR 3155142 A1 20250516 FR 2312399 A 20231113 L’objet de la présente invention est de proposer une méthode semi-analytique permettant de simuler précisément la transmission d’un champ ultrasonore à travers une paroi crânienne au moyen d’algorithmes présentant des temps de calcul compatibles dans un contexte préopératoire. Plus généralement l’invention consiste à simuler un champ ultrasonore à travers une ou plusieurs parois pour focaliser un faisceau en un point situé derrière les parois ou encore pour diriger un faisceau dans une direction prédéterminée derrière cette paroi. Figure 1 Méthode de détermination d’un champ d’une onde ultrasonore à travers plusieurs parois L’invention concerne le domaine de la simulation de la propagation d’ondes ultrasonores dans des milieux solides hétérogènes et/ou anisotropes de géométrie complexe. L’invention s’applique notamment à la thérapie ultrasonore transcrânienne mais aussi à l’imagerie par ultrasons. Elle porte plus précisément sur une méthode de simulation d’une onde ultrasonore émise par un transducteur à ultrasons à travers la paroi crânienne en vue de maîtriser le faisceau ultrasonore transmis, et notamment la zone focale, avant de le générer. La thérapie ultrasonore cérébrale ou transcrânienne consiste en l'utilisation d’ultrasons focalisés pour traiter certaines maladies. La guérison des maladies cérébrales reste très difficile, principalement en raison du mauvais accès des agents pharmacologiques au cerveau en raison des barrières hémato-encéphalique et hémato-tumorale (BHE et BHT). La focalisation des ultrasons pulsés de faible intensité dans le cerveau, combinée à l’injection de microbulles augmente considérablement la distribution de médicaments aux tissus cérébraux, avec un effet thérapeutique établi dans de nombreux modèles animaux ainsi que chez les patients. Cette perméabilisation de la BHE ou de la BHT est locale et transitoire à condition que l'intensité du faisceau soit bien maîtrisée car l'implosion des microbulles pourrait entraîner des hémorragies. Par conséquent, afin d'assurer la sécurité du protocole, il est nécessaire de maitriser précisément la dosimétrie du champ acoustique après le franchissement de la paroi osseuse La thérapie par ultrasons focalisés de haute intensité (HIFU) est une approche innovante pour le traitement des tumeurs intracrâniennes. Dans un avenir proche, cette thérapie pourrait être une alternative intéressante à la chirurgie cérébrale et à la radiothérapie. La thérapie HIFU induit une élévation de la température au point focal du faisceau ultrasonore et un volume tissulaire défini peut être détruit thermiquement (ablation thermique). Pour une fréquence et des paramètres acoustiques donnés, le volume d'ablation dépend de la montée en température et du temps d'exposition (dose thermique). Cette technique nécessite une maîtrise parfaite de la zone focale du faisceau ultrasonore transmis, afin de cibler la zone à détruire (tumeur) tout en assurant son innocuité sur les tissus sains environnants. Dans les différences applications envisagées ci-dessus ainsi que pour des applications d’imagerie par ultrasons, il est nécessaire de maitriser précisément la propagation des ondes ultrasonores d’autant qu’elle varie selon la géométrie et les propriétés mécaniques de la paroi crânienne. En effet, les faisceaux ultrasonores incidents peuvent être déviés, atténués ou aberrés, voire une combinaison des deux ou trois de ces cas à la fois. De ce fait la boite crânienne est alors considérée comme un milieu complexe perturbant fortement la focalisation des ultrasons pour la thérapie, ou la reconstruction des images dans le cas de l’imagerie. Ce problème, inhérent au crâne, est principalement dû à une géométrie et à des propriétés acoustiques locales très différentes chez un même individu, mais aussi dû à une variation de ces mêmes propriétés en fonction des individus. La modélisation personnalisée du crâne prend alors toute son importance pour simuler la propagation des ultrasons et améliorer les outils de thérapie et d’imagerie. Habituellement, les simulations de champs ultrasonores générés vers des parois transcrâniennes sont réalisées par des méthodes numériques qui sont considérées comme étant plus précises, bien que plus couteuses en temps d’exécution et en espace mémoire, que les méthodes semi-analytiques. La référence [1] recense, dans une étude bibliographique, les principales méthodes de simulation de la propagation de champs ultrasonores au travers d’une paroi crânienne, qui sont regroupées en deux grandes catégories. La première catégorie de méthodes concerne les méthodes numériques, par exemple du type éléments finis, différences finies ou méthodes aux éléments de frontières, qui sont considérées comme la référence. En effet, puisqu’elles consistent en une résolution numérique exacte des équations de propagation des ondes ultrasonores, ces méthodes permettent de simuler l’ensemble des phénomènes physiques complexes mis en jeu. Cependant, elles demeurent dans la plupart des cas coûteuses en temps de calcul et en place mémoire occupée. Par ailleurs, elles comportent des contraintes fortes sur la zone de calcul. En effet, celle-ci doit être suffisamment échantillonnée spatialement, pour respecter des critères de convergence numérique, et elle doit contenir à la fois le terme source, c’est-à-dire la sonde ultrasonore, et la zone d’intérêt. Enfin, la mise en œuvre nécessite généralement une expertise importante de la part de l’utilisateur, notamment au travers de l’utilisation d’outils de maillage de la scène, comprenant à la fois la sonde et le crâne. Un exemple d’étude comparative des principaux outils de simulation numérique est donné dans la référence [2]. La seconde catégorie concerne les méthodes dites semi-analytiques, notamment des méthodes du type lancer de rayon, qui peuvent également être mises en œuvre pour des simulations ultrasonores transcrâniennes. Ces méthodes sont plus rapides, moins coûteuses en charge mémoire, et ne comportent pas de contraintes quant à la zone de calcul. En revanche, puisqu’elles sont fondées sur des approximations asymptotiques, elles ne permettent pas de simuler de manière satisfaisante toute la complexité des phénomènes physiques mis en jeu lors de la propagation. Néanmoins, des études récentes ont montré que ces méthodes semi-analytiques permettaient d’obtenir des résultats aussi précis que les méthodes numériques. L’objet de la présente invention est de proposer une méthode semi-analytique permettant de simuler précisément la transmission d’un champ ultrasonore à travers une paroi crânienne au moyen d’algorithmes présentant des temps de calcul compatibles dans un contexte préopératoire. Plus généralement l’invention consiste à simuler un champ ultrasonore à travers une ou plusieurs parois pour focaliser un faisceau en un point situé derrière les parois ou encore pour diriger un faisceau dans une direction prédéterminée derrière cette paroi. L’invention a pour objet une méthode, mise en œuvre par ordinateur, de détermination d’un champ d’une onde ultrasonore émis par un transducteur ultrasonore comprenant au moins un élément ultrasonore vers une zone d’intérêt, le transducteur étant séparé de la zone d’intérêt par plusieurs milieux séparés par des parois, la méthode comprenant les étapes de : Recevoir un modèle de chaque paroi sous la forme d’une surface 3D paramétrée obtenue par approximation de B-splines multi-niveaux, Recevoir une décomposition de la surface de chaque élément ultrasonore en une pluralité de surfaces élémentaires étant chacune décrite par un ensemble de points comprenant un point central et quatre points correspondants aux sommets d’un carré centré sur le point central, Déterminer un champ ultrasonore transmis par chaque élément ultrasonore et pour chaque point d’intérêt de la zone d’intérêt au moyen des sous-étapes de : Pour chaque mode de propagation de l’onde ultrasonore parmi plusieurs modes et pour chaque point source de l’ensemble de points caractérisant une surface élémentaire parmi l’ensemble des surfaces élémentaires décrivant chaque élément ultrasonore, déterminer un trajet entre le point source et le point d’intérêt de sorte à minimiser le temps de vol de l’onde ultrasonore entre ces deux points, déterminer l’atténuation de l’onde ultrasonore sur ce trajet, calculer la réponse impulsionnelle du champ ultrasonore pour chaque surface élémentaire à partir des trajets et atténuations déterminés pour l’ensemble de points associés, Sommer l’ensemble des réponses impulsionnelles calculées, Déterminer le champ ultrasonore transmis par l’élément ultrasonore à partir du résultat de la somme. Selon un aspect particulier de l’invention, le transducteur ultrasonore comprend plusieurs éléments ultrasonores, la méthode comprenant en outre une somme des champs ultrasonores de chaque élément pondérée par une loi d’amplitudes prédéfinie, chaque champ ultrasonore étant retardé via une loi de retards ou de phases prédéfinie. Dans une variante de réalisation, la méthode selon l’invention comprend une étape préalable de décomposition de la surface de chaque élément ultrasonore comprenant au moins: Décomposer la surface en une pluralité de surfaces élémentaires contiguës, chaque surface élémentaire étant centrée sur un point central, Pour chacun des points centraux, former un carré centré sur le point central dans un plan tangent à la surface, Caractériser la surface élémentaire par l’ensemble de points comprenant le point central et les quatre sommets du carré, Dans une variante de réalisation, la méthode selon l’invention comprend une étape préalable de détermination d’un modèle de chaque paroi comprenant au moins : Acquérir une image des parois, Extraire des maillages des parois sous la forme d’un nuage de points, Appliquer une approximation par B-splines multi-niveau au nuage de points de sorte à obtenir une surface paramétrée 3D lissée Selon un aspect particulier de l’invention, les milieux séparés par des parois sont anisotropes et l’étape de déterminer un trajet entre le point source et le point d’intérêt est réalisée en prenant en compte l’anisotropie de chaque milieu traversé. Selon un aspect particulier de l’invention, l’étape de déterminer un trajet entre le point source et le point d’intérêt est réalisé au moyen d’un algorithme de minimisation de la somme des segments composant le trajet, chaque segment correspondant à un milieu traversé, la somme étant pondérée par les lenteurs de chaque milieu. Selon un aspect particulier de l’invention, les modes de propagation de l’onde ultrasonore comprennent au moins un mode de propagation impliquant des réflexions de l’onde ultrasonore sur au moins certaines parois et la prise en compte du coefficient de réflexion du trajet issu du point central de chaque surface élémentaire dans le calcul de la réponse impulsionnelle du champ ultrasonore. Selon un aspect particulier de l’invention, la décomposition d’une surface d’un élément en surfaces élémentaires est réalisée à l’aide d’un diagramme de Voronoï. Selon un aspect particulier de l’invention, le calcul de la réponse impulsionnelle du champ ultrasonore comporte au moins le calcul de l’amplitude du champ et des bornes minimales et maximales des temps de vol correspondants aux trajets calculés. Selon un aspect particulier de l’invention, la zone d’intérêt est une région transcrânienne et les parois sont les parois d’un crâne. Selon un aspect particulier de l’invention, l’étape de déterminer un champ ultrasonore comprend au moins l’étape de convoluer le résultat de la somme avec un signal de référence associé à l’élément ultrasonore ou d’extraire l’amplitude et la phase d’une composante fréquentielle prédéfinie de chaque réponse impulsionnelle pour obtenir un champ ultrasonore transmis par ledit élément. L’invention a aussi pour objet un programme d’ordinateur comportant des instructions pour l'exécution de la méthode selon l’invention, lorsque le programme est exécuté par un processeur. L’invention a aussi pour objet un support d'enregistrement lisible par un processeur sur lequel est enregistré un programme comportant des instructions pour l'exécution de la méthode selon l’invention, lorsque le programme est exécuté par un processeur. L’invention a aussi pour objet un simulateur de champ d’une onde ultrasonore comprenant un processeur et une mémoire configurés pour exécuter les étapes de la méthode selon l’invention et un dispositif d’affichage pour afficher le champ calculé dans une zone d’intérêt. D’autres caractéristiques et avantages de la présente invention apparaîtront mieux à la lecture de la description qui suit en relation aux dessins annexés suivants. FIG. 1 représente un organigramme détaillant les étapes de mise en œuvre d’une méthode de détermination d’un champ ultrasonore émis à travers au moins une paroi selon un mode de réalisation de l’invention FIG. 2 représente un organigramme illustrant une étape de calcul de trajets ultrasonores et d’atténuation selon un mode de réalisation de l’invention, FIG. 3 représente un organigramme illustrant une étape de calcul d’une réponse impulsionnelle d’une onde ultrasonore émise, selon un mode de réalisation de l’invention, FIG. 4 représente un schéma illustrant un exemple de décomposition d’une sonde ultrasonore en plusieurs éléments, FIG. 5 représente un schéma illustrant le principe de la méthode de calcul de trajets de la FIG. 2 , FIG. 6 représente un schéma illustrant le principe de calcul d’un champ ultrasonore par la méthode des pinceaux, FIG. 7 représente un exemple de réponse impulsionnelle d’un champ ultrasonore, La FIG. 1 illustre, sur un organigramme, les principales étapes d’une méthode de détermination d’un champ ultrasonore émis depuis un point source d’une surface d’une sonde ultrasonore vers un point cible situé derrière une ou plusieurs parois. Dans la suite de la description, l’invention est présentée dans le contexte d’une paroi crânienne pour des applications de thérapie ou d’imagerie transcrânienne mais l’invention n’est pas limitée à ce type de paroi et peut s’appliquer à d’autres types de parois d’objets notamment pour des applications de contrôle non destructif. La méthode comporte une phase de prétraitement 100 suivie d’une phase de calcul de champ 110. La phase de prétraitement 100 prend en entrée un modèle d’un crane 101 qui est par exemple obtenu via une technique d’imagerie médicale par rayons X et un modèle de la géométrie d’une sonde ultrasonore 102 destinée à être utilisée pour générer un champ ultrasonore. A l’étape 103, on extrait les propriétés acoustiques du crâne à partir du modèle obtenu. Cette étape est par exemple réalisée au moyen des techniques d’extraction décrites dans les références [3] ou [4]. A l’étape 104, on extrait des maillages des parois du crâne, par exemple au moyen d’une technique de décomposition d’une surface en maillages. Les maillages sont par exemple constitués de nuages de points en 3D. Les éléments de description du crâne obtenus aux étapes 103 et 104 servent à modéliser les surfaces externes constituant les parois du crâne à l’étape 105. Afin de modéliser les différentes parois du crâne par des surfaces lisses, une méthode d’approximation par B-splines est utilisée à l’étape 106, par exemple à partir des techniques décrites dans la référence [5]. Un avantage à cette méthode est qu’elle consiste en une approximation et non une interpolation, ce qui permet d’être robuste à d’éventuels points aberrants du nuage de points obtenus à l’extraction 104. Cette méthode permet d’obtenir des surfaces plus lisses qu’avec une technique basée sur une interpolation, ce qui permet d’améliorer la convergence des algorithmes de calcul de trajets ultrasonores décrits par la suite. Elle présente également l’avantage de réduire le bruit dans le nuage de points. La méthode d’approximation est par exemple une méthode dite de multi-niveaux de B-spline ou « MBA ». A l’issue de l’étape 106, chaque paroi du crâne est modélisée via une surface 3D qui constitue une interface entre deux milieux d’indice i et i+1 et qui est décrite par une fonction où correspond aux coordonnées paramétriques de la surface décrite au moyen de splines. Dans le cas de la méthode MBA, on a . Les interfaces définissant les parois du crâne sont au moins au nombre de deux mais peuvent être en nombre supérieur à deux. La phase de prétraitement comporte aussi une phase de modélisation de la sonde ultrasonore. La sonde utilisée peut être un transducteur ultrasonore mono ou multi-éléments. L’invention est adaptée à n’importe quel type de géométrie de surface 2D ou 3D pour chaque élément du transducteur. A l’étape 106, chaque élément de la surface émettrice de la sonde est échantillonné spatialement pour obtenir un ensemble de points, répartis de manière uniforme sur la totalité de la surface de l’élément. Le nombre de points est un paramètre d’entrée de la méthode. Cet échantillonnage peut être obtenu par exemple à partir d’un diagramme de Voronoï. La somme des surfaces élémentaires associées à chaque point est égale à la surface totale de l’élément. La FIG. 4 illustre un exemple d’élément 400 d’un transducteur dont la forme est circulaire. L’élément 400 est décomposé en plusieurs surfaces élémentaires 401 de centres C. A l’étape 107, on détermine, pour chaque surface élémentaire, les coordonnées de 4 sommets d’un carré de centre égal au point C qui est le centre de la surface élémentaire. Le carré virtuel ainsi formé a une surface A de même aire que celle de l’élément 401. Dans le cas d’une surface conformée en trois dimensions, le carré virtuel est formé dans un plan ( , ) d’un repère orthonormé direct ( , , ) avec un vecteur normal à la surface élémentaire 401. Cette étape est illustrée sur la droite de la figure 4 qui montre, pour l’exemple de l’élément 401, le carré virtuel de centre C et dont les quatre sommets sont définis, dans un repère orthonormé direct ( , , ) par les coordonnées : Ces quatre points sont ainsi définis pour chaque surface élémentaire pour servir de points sources dans le calcul de trajets préalables aux calculs de champs. Cette modélisation des éléments de la sonde permet d’adresser n’importe quelle géométrie de sonde. A l’issue de la phase de prétraitement 100, on obtient donc un modèle des parois du crâne via des fonctions paramétriques en trois dimensions d’une part et d’autre part les coordonnées de cinq points (le centre et les quatre sommets du carré virtuel) caractérisant chaque surface élémentaire issue de la décomposition d’un élément de la sonde ultrasonore. L’ensemble de ces paramètres est fourni en entrée d’une seconde phase 110 de calcul de champs ultrasonores. L’exemple particulier de la FIG. 1 est donné pour le calcul du champ ultrasonore correspondant à un élément ultrasonore d’un transducteur. Si le transducteur comporte plusieurs éléments ultrasonores, une étape de combinaison supplémentaire doit être appliquée à partir des champs ultrasonores calculés pour chacun des éléments. Cette étape sera décrite plus en détail par la suite. La méthode de calcul de champ ultrasonore objet de la seconde phase 110 comporte trois boucles d’itérations. La première boucle vise à parcourir l’ensemble des points de la zone de calcul, c’est-à-dire les points du cerveau dans le cas d’une application transcrânienne ou les points d’une zone d’un objet à contrôler dans le cas d’une application de contrôle non destructif. Pour chaque point cible de la zone de calcul, on calcule un champ ultrasonore émis à travers les parois modélisées à l’étape 105. Ce calcul nécessite un premier calcul de trajets par minimisation de temps de vol 111 et un deuxième calcul de réponse impulsionnelle 112, ces deux étapes étant itérées sur une deuxième boucle d’itération parcourant plusieurs modes de propagation des ondes ultrasonores et une troisième boucle d’itération parcourant chaque point caractérisant chaque surface élémentaire de chaque élément de la sonde. Ces points correspondent à ceux déterminés à l’étape 107. La FIG. 2 détaille le calcul 111 de trajets ultrasonores et la FIG. 3 détaille le calcul 112 de la réponse impulsionnelle associée. Le but du calcul 111 de trajets ultrasonores consiste à rechercher le trajet le plus court entre un point A de la sonde et un point B cible situé de l’autre côté de la paroi crânienne modélisée par plusieurs surfaces. Ce principe est illustré à la FIG. 5 pour un exemple de paroi modélisée par deux surfaces f 1 et f 2 séparant trois milieux caractérisés par des vitesses de propagation de l’onde ultrasonore respectivement égales à c 1 , c 2 et c 3 . Le but de l’étape de calcul 111 est de trouver les points d’intersection I 1 et I 2 du trajet ultrasonore émis depuis A vers B qui permettent d’obtenir le temps de vol le plus court. Notons I1 = (x1,y1,z1) le point d’interface 1 et I2 = (x2,y2,z2) de l’interface 2 par lesquels passe le rayon, on a et car les interfaces sont décrites par des splines. Les 4 paramètres à faire varier sont donc : L’objectif du calcul de trajet 111 est de minimiser le temps de parcours donné par la relation : . L’étape de calcul de trajet 111 reçoit en entrée un point source 210, c’est-à-dire un point de la surface de la sonde émettrice et un point de calcul 211, c’est-à-dire un point de la zone de calcul. Il reçoit également en entrée les interfaces 213 modélisant la paroi à traverser ainsi que les vitesses 212 associées à chaque milieu séparé par une interface. L’étape 201 consiste en un tracé de rayon entre le point source et le point de calcul de sorte à définir les points d’intersection du rayon avec les différentes interfaces. Les points d’intersection appartiennent aux surfaces approximées obtenues à partir de splines. L’étape 201 peut être réalisée au moyen d’un algorithme de tracés de rayons tel que décrit par exemple dans la demande de brevet FR3114032 du Demandeur. L’étape 202 consiste ensuite à rechercher les coordonnées des points d’intersection qui minimisent le temps de vol donné par la relation précédente. Cette étape d’optimisation peut être réalisée par tout algorithme d’optimisation approprié, par exemple un algorithme de minimisation quasi-Newtonienne comme par exemple l’algorithme L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shannon), ou par la méthode décrite dans la demande de brevet FR3114032. L’objectif général de l’étape d’optimisation 202 est de rechercher le trajet entre le point source et le point de calcul au travers de interfaces tel que la fonction temps de vol soit minimale. Comme chaque milieu est homogène de lenteur le temps de vol est la somme des temps de vol dans chaque milieu traversé par l’onde. Pour trouver ce minimum, cela revient à trouver les points d’intersection du trajet ultrasonore sur chacune des interfaces. On utilise donc un algorithme de minimisation tel que l’algorithme L-BFGS pour résoudre le problème suivant : avec : est la surface paramétrée qui décrit l’interface entre les milieux d’indice i et i+1. Pour des milieux anisotropes, la lenteur s i dépend de la direction de propagation et la relation précédente s’écrit alors : La prise en compte de milieux anisotropes dans le calcul du trajet de l’onde ultrasonore permet de mieux caractériser l’objet traversé par l’onde, en particulier le crâne humain qui présente des propriétés spécifiques à chaque individu. Après optimisation, on obtient une description 203 du trajet ultrasonore entre le point source et le point de calcul. Selon un autre aspect de l’invention, l’atténuation de l’onde à travers chaque milieu traversé est également calculée à l’étape 204. Pour calculer l’atténuation totale sur un trajet, on somme les atténuations sur chacune des portions du trajet dans chaque milieu traversé. Si chaque milieu a un coefficient d’atténuation à la fréquence , et un exposant d’atténuation , l’atténuation totale est obtenue via la relation suivante: Les valeurs des atténuations de chaque milieu sont des données d’entrée de la méthode. Dans un mode de réalisation de l’invention, les calculs de trajet 111 sont réalisés pour différents modes de propagation. En particulier certains modes de propagation impliquent des réflexions sur les interfaces internes de la structure hétérogène traversée (par exemple le crâne). Dans ce cas, les réflexions sont considérées comme des segments de trajets supplémentaires. On applique alors la même méthode que celle décrite précédemment en ajoutant autant d’interfaces supplémentaires que de réflexions considérées qui correspondent à des propagations multiples de l’onde dans le même milieu. Selon d’autres modes de propagation, les valeurs des vitesses et/ou atténuation des milieux traversés peuvent différer d’un mode à l’autre. Les trajets ainsi calculés à l’étape 111 pour chaque couple de points (point source, point de calcul) et chaque mode de propagation sont ensuite transmis à une étape 112 de calcul de réponse impulsionnelle du champ ultrasonore émis Plus précisément, les trajets sont calculés pour chaque centre d’une surface élémentaire d’un élément de la sonde ainsi que pour les quatre sommets du carré virtuel formé à partir de ce centre (tel qu’explicité à l’étape 107). La FIG. 3 détaille la mise en œuvre de l’étape 112 de calcul de réponse impulsionnelle. Cette étape de calcul reçoit en entrée, pour chaque surface élémentaire d’un élément de la sonde, le trajet issu du point central 301 ainsi que les trajets 302 issus des quatre sommets du carré virtuel. Le calcul de la réponse impulsionnelle 112 est en partie basé sur la méthode dite des pinceaux décrite dans la référence [6]. On décrit tout d’abord le principe de la méthode des pinceaux. Pour reconstituer le front d’onde émis et donc calculer l’intégrale de Rayleigh associée au champ ultrasonore, il est possible d’utiliser la méthode des pinceaux. Physiquement, un pinceau représente un flux d’énergie acoustique en provenance de la source d’énergie ponctuelle O et se propageant en direction de la surface S. Puisque O est la seule source d’énergie considérée dans la modélisation, la loi de conservation de l’énergie nous indique que le flux d’énergie total traversant la surface du pinceau est nul. Ainsi, l’énergie surfacique traversant la surface S est inversement proportionnelle à son aire. Le rayon central du pinceau suit la trajectoire de Fermat, obtenue à partir de l’étape de calcul de trajet 111. Le rayon paraxial suit une direction différente de celle du rayon central, mais que l’on considère comme assez proche pour que l’approximation paraxiale soit valable (l’approximation paraxiale étant une approximation d’ordre 1, l’erreur augmente à mesure que l’on s’éloigne de l’axe). On le caractérise par quatre grandeurs différentielles exprimées dans un repère inclus dans le plan orthogonal à l’axe central du pinceau, permettant de définir un rayon paraxial du pinceau par un vecteur à 4 coordonnées: dx et dy caractérisent la position du point d’intersection entre le rayon paraxial et la surface S du pinceau. dSx et dSy donnent la projection du vecteur lenteur du rayon paraxial sur la surface S. Dans le cadre de l’approximation paraxiale, l’évolution d’un vecteur de pinceau dans un milieu homogène peut être caractérisée par une matrice 4 × 4 dépendante du milieu de propagation. On peut de même définir des matrices 4 × 4 caractéristiques d’interfaces entre deux matériaux. On appelle ces matrices les matrices de propagation et on les note usuellement avec une notation par blocs : où les matrices A, B, C et D sont des matrices de taille 2 × 2. En considérant un pinceau traversant une suite de milieux homogènes et d’interfaces caractérisés par les matrices de propagation (L i ) pour 1≤i≤n, le pinceau résultant peut s’écrire sous la forme : Puisqu’un milieu hétérogène peut être modélisé par une suite de milieux homogènes et d’interfaces les séparant, ces matrices suffisent à décrire l’évolution d’un pinceau lors de la traversée d’un milieu hétérogène. L’utilisation des pinceaux pour le calcul de champ ultrasonore permet de simuler la propagation spatiale d’une onde ultrasonore. Afin d’obtenir l’amplitude finale du pinceau, il faut tenir compte de sa divergence. En effet, la loi de conservation de l’énergie à l’intérieur du pinceau fait que le flux d’énergie diminue avec l’ouverture du pinceau, celui-ci étant, dans le cadre de l’approximation d’un pinceau à faible ouverture, réparti uniformément sur la surface de base S. Le rapport entre l’intensité surfacique finale I S et l’intensité par unité d’angle solide initiale IΩ s’écrit donc comme le rapport entre l’angle solide de départ dΩ et la surface de base finale du pinceau dS : En plus du facteur d’amplitude lié à la divergence du pinceau au cours de sa propagation, chaque interface traversée réduit l’énergie du pinceau puisque différents modes sont générés à chaque interface. La répartition de l’énergie entre les différents modes de propagation (transmis, réfléchis et évanescents pour chaque type de polarisation) est obtenue à partir des conditions de continuité aux interfaces et se traduit par le calcul d’un coefficient de Fresnel. Pour obtenir la valeur du coefficient de transmission en amplitude, applicable directement à la valeur du déplacement et donc de l’amplitude du champ ultrasonore, le rapport entre les sections de pinceau avant et après l’interaction est à considérer. Pour chaque pinceau émis depuis le point de champ considéré et dont la surface finale intersecte la surface du traducteur ultrasonore, il est possible d’extraire une réponse impulsionnelle élémentaire. La sommation de ces réponses élémentaires constitue la réponse impulsionnelle globale qui sert ensuite à déterminer la valeur d’amplitude du champ au cours du temps. La réponse impulsionnelle élémentaire d’un pinceau est assimilée à un échelon temporel dont l’amplitude, l’étalement temporel et le temps central sont extraits des caractéristiques finales du pinceau. La valeur d’amplitude correspond à une énergie. Par conséquent, elle est distribuée sur l’étalement temporel du pinceau, ce qui revient à ajouter un facteur inverse de l’étalement temporel à la valeur d’amplitude. L’amplitude A est donc donnée par la relation Avec A0 l’amplitude initiale, DF le facteur de divergence, T A le coefficient de transmission en amplitude, dS la surface du pinceau à l’arrivée, Δt l’étalement temporel. L’étalement temporel Δt est calculé comme la différence de marche maximale évaluée sur le pinceau, en fonction de l’angle d’incidence du pinceau sur le traducteur et de son étendue spatiale. Le temps central est égal au temps de vol du rayon central du pinceau. Pour un traducteur multiéléments utilisant des lois de retard ou d’amplitude, le temps central et l’amplitude de chaque pinceau sont modifiés en cohérence avec l’élément du traducteur intersecté. Les coefficients de Fresnel pouvant être complexes, l’amplitude du déplacement est elle-même potentiellement complexe. De plus, le déplacement étant un champ vectoriel, la réponse calculée est un signal à 3 dimensions complexes, soit 6 dimensions réelles. Une fois les réponses impulsionnelles élémentaires calculées pour un point de champ donné, nous les sommons sur une réponse impulsionnelle unique. Celle-ci correspond au champ rayonné par le traducteur vers le point de champ lorsque chaque élément du traducteur émet une impulsion de Dirac temporelle. Par conséquent, étant donné un signal d’entrée S 0 auquel est soumis le traducteur ultrasonore, le champ ultrasonore en un point de champ s’obtient par convolution des signaux complexes de réponse impulsionnelle avec le signal d’entrée. Pour chaque composante x, y et z du déplacement, la partie réelle est convoluée avec le signal de référence tandis que la partie complexe est convoluée avec la transformée de Hilbert du signal de référence. Le calcul de chaque composante scalaire si(t) du champ s’écrit donc, en notant RIi(t) la réponse impulsionnelle : si(t) = RIi(t) ∗ S0(t) Le modèle paraxial de pinceaux est basé sur les matrices de propagation pour proposer une approximation linéaire des rayons paraxiaux autour du rayon central d’un pinceau. Une approximation des rayons inclus dans le pinceau peut également être obtenue à l’aide d’une interpolation. Pour ce faire, on simule la propagation des quatre rayons r 0 , r 1 , r 2 et r 3 délimitant le pinceau en plus de propager le rayon central rc, comme cela est présenté sur la FIG. 6 . Comme explicité précédemment, les quatre rayons r 0 , r 1 , r 2 et r 3 correspondent aux quatre trajets issus des quatre sommets du carré virtuel centré sur le centre d’une surface élémentaire obtenue par décomposition d’un élément du transducteur ultrasonore. Cette approche permet de s’affranchir du calcul des matrices de propagation le long du trajet du pinceau tout en rendant accessible le calcul des grandeurs nécessaires à l’intégration de la réponse impulsionnelle élémentaire associée au pinceau. La fenêtre temporelle du créneau modélisant le pinceau est déterminée de sorte à encadrer l’ensemble des temps d’arrivée des rayons composant le pinceau. Ainsi, les bornes t min et t max du créneau temporel représentant la réponse impulsionnelle du pinceau sont déterminées respectivement comme l’instant minimum et l’instant maximum parmi l’ensemble des temps de vol calculés. La FIG. 7 représente le créneau temporel obtenu. En utilisant le modèle de pinceau à quatre rayons, la surface finale du pinceau est l’aire délimitée par les points d’arrivée des quatre rayons tel que représenté à la FIG. 6 . L’angle solide initial est quant à lui défini par la surface du pinceau de rayons à une distance unitaire. Le coefficient de divergence peut être calculé à partir de ces grandeurs et de la lenteur initiale S 0 du rayon axial via la relation suivante : De cette façon, en considérant le coefficient de transmission en amplitude cumulé T A du rayon central et l’amplitude initiale A 0 , on peut calculer de la même façon que précédemment l’amplitude du déplacement au point de champ considéré en la répartissant sur l’étalement temporel du pinceau via la relation . Ce modèle simplifié de pinceau permet de reproduire une réponse impulsionnelle élémentaire de pinceau sans avoir à calculer de matrices paraxiales, supprimant la nécessité de calculs matriciels coûteux (particulièrement en milieux anisotropes) et les remplaçant par des opérations de lancer de rayons ultrasonores. Pour un pinceau isolé, cette approche nécessite un total de 5 rayons ultrasonores pour produire une réponse impulsionnelle élémentaire. L’étape de calcul de réponse impulsionnelle 112 est donc réalisée sur la base de la méthode des pinceaux décrite ci-dessus. A l’étape 303, un calcul de polarisation est réalisé à partir du trajet central 301. A l’étape 304, on calcule un coefficient de transmission ou de réflexion cumulé T A (selon le mode de propagation) du trajet central. Dans le cas où le mode de propagation choisi implique à la fois des transmissions et réflexions sur les différentes interfaces, le coefficient cumulé est calculé à partir des coefficients de transmissions et réflexions sur chacune des interfaces considérées. A l’étape 305, on calcule les temps de vol correspondants aux 5 trajets 301,302 reçus en entrée. A l’étape 306, on calcule l’angle solide formé par les quatre trajets r 0 , r 1 , r 2 et r 3 correspondant aux sommets du carré virtuel, au point de calcul. A l’étape 307, on calcule le coefficient de divergence DF à partir de l’angle solide via la relation A l’étape 308, on calcule les valeurs de t min et t max de la façon décrite précédemment. Enfin, à l’étape 309, on calcule l’amplitude A via la formule donnée ci-dessus. La réponse impulsionnelle pour un couple (point source, point de calcul) et un mode de propagation est donnée par les bornes du créneau temporel, son amplitude et sa polarisation. On retourne ensuite à la FIG. 1 . Le calcul se poursuit à l’étape 113 par la détermination de la taille globale du signal qui est obtenue en relevant les instants temporels minimum et maximum parmi tous les créneaux temporels calculés pour les deux itérations (mode de propagation, point source). Cette étape est nécessaire pour allouer l’espace mémoire suffisant pour le signal de sortie. A l’étape 114, on somme l’ensemble des réponses impulsionnelles calculées. A l’étape 115, on effectue une convolution de la réponse impulsionnelle avec le signal de référence ou un calcul en monochromatique. Dans le cas le plus général, le signal de référence est celui utilisé pour générer le signal ultrasonore par la sonde. Dans le cas du calcul en monochromatique, une réalisation consiste à calculer uniquement l’amplitude et la phase du spectre mesuré à la fréquence souhaitée. Enfin, à l’étape 116 on obtient le champ de pression en tout point de la zone de calcul pour un élément ultrasonore. Dans le cas où le transducteur est multi-éléments, l’algorithme 110 de calcul du champ ultrasonore doit être itéré pour chaque élément du transducteur. Une étape supplémentaire (non représentée à la FIG. 1 ) doit également être appliquée afin de combiner les différents champs émis par chaque élément en un champ global associé au transducteur. Pour cela, une loi d’amplitude et une loi de retard prédéfinies selon l’application visée sont prises en compte dans cette combinaison. Plus précisément, cette étape supplémentaire consiste à calculer la somme des champs ultrasonores de chaque élément, pondérée à partir de la loi d’amplitude, chaque champ étant retardé à partir de la loi de retard. Autrement dit, la loi d’amplitude et la loi de retard donnent les amplitudes (respectivement retards) à appliquer à chaque champ correspondant à chaque élément. L’invention présente plusieurs avantages par rapport aux méthodes de l’art antérieur. Elle est adaptée à la prise en compte de sondes ultrasonores présentant des surfaces de géométrie quelconque, ce qui présente un avantage particulier dans le cadre des applications transcrâniennes pour lesquelles la surface émettrice de la sonde peut avoir un profil sphérique ou plus généralement conformé en 3D pour s’adapter à la forme du crâne. De même, dans le cas général, chaque élément a une surface projetée spécifique et de géométrie complexe, par exemple du type polygonal. En effet, la décomposition en surfaces élémentaires de chaque élément de géométrie quelconque couplé à la formation de quatre points d’intérêts associés au point central de l’élément permet de s’adapter à n’importe quelle forme géométrique. Par ailleurs, l’invention permet également de tenir compte de milieux anisotropes dans le calcul des trajets et de prendre en compte l’atténuation de chaque milieu dans le calcul du champ ultrasonore. Enfin une meilleure modélisation des surfaces définissant le crâne permet également une amélioration de la précision et du temps de convergence dans le calcul des trajets ultrasonores. L’invention peut être mise en œuvre en tant que programme d’ordinateur comportant des instructions pour son exécution de sorte à permettre une simulation des champs ultrasonores générés pour chaque point d’une zone ciblée. Le programme d’ordinateur peut être enregistré sur un support d’enregistrement lisible par un processeur. Le simulateur peut comprendre un processeur et une mémoire pour exécuter la méthode de génération d’un champ d’une onde ultrasonore selon l’invention ainsi qu’un dispositif d’affichage pour afficher le champ généré pour une zone d’intérêt, par exemple une zone transcrânienne. Un tel outil de simulation permet une assistance à un opérateur préalable à la génération du champ ultrasonore au moyen d’un transducteur de sorte à permettre une meilleure calibration du champ avant sa génération ou, dans le cas d’une sonde multiéléments, à calculer les lois de phase, ou de retard, permettant de corriger les aberrations subies par le champ lors de la traversée de la paroi crânienne Méthode, mise en œuvre par ordinateur, de détermination d’un champ d’une onde ultrasonore émis par un transducteur ultrasonore comprenant au moins un élément ultrasonore vers une zone d’intérêt, le transducteur étant séparé de la zone d’intérêt par plusieurs milieux séparés par des parois, la méthode comprenant les étapes de : Recevoir un modèle (105,213) de chaque paroi sous la forme d’une surface 3D paramétrée obtenue par approximation de B-splines multi-niveaux, Recevoir une décomposition (107) de la surface (400) de chaque élément ultrasonore en une pluralité de surfaces élémentaires (401) étant chacune décrite par un ensemble de points comprenant un point central et quatre points correspondants aux sommets d’un carré centré sur le point central, Déterminer un champ ultrasonore transmis par chaque élément ultrasonore et pour chaque point d’intérêt de la zone d’intérêt au moyen des sous-étapes de : Pour chaque mode de propagation de l’onde ultrasonore parmi plusieurs modes et pour chaque point source de l’ensemble de points caractérisant une surface élémentaire parmi l’ensemble des surfaces élémentaires décrivant chaque élément ultrasonore, déterminer (111) un trajet entre le point source et le point d’intérêt de sorte à minimiser le temps de vol de l’onde ultrasonore entre ces deux points, déterminer (204) l’atténuation de l’onde ultrasonore sur ce trajet, calculer (112) la réponse impulsionnelle du champ ultrasonore pour chaque surface élémentaire à partir des trajets et atténuations déterminés pour l’ensemble de points associés, Sommer (114) l’ensemble des réponses impulsionnelles calculées, Déterminer (115) le champ ultrasonore transmis par l’élément ultrasonore à partir du résultat de la somme. Méthode de détermination d’un champ d’une onde ultrasonore selon la revendication 1 dans laquelle le transducteur ultrasonore comprend plusieurs éléments ultrasonores, la méthode comprenant en outre une somme des champs ultrasonores de chaque élément pondérée par une loi d’amplitudes prédéfinie, chaque champ ultrasonore étant retardé via une loi de retards ou de phases prédéfinie. Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes comprenant une étape préalable de décomposition de la surface de chaque élément ultrasonore comprenant au moins: Décomposer (106) la surface en une pluralité de surfaces élémentaires contiguës, chaque surface élémentaire étant centrée sur un point central, Pour chacun des points centraux, former un carré centré sur le point central dans un plan tangent à la surface, Caractériser (107) la surface élémentaire par l’ensemble de points comprenant le point central et les quatre sommets du carré, Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes comprenant une étape préalable de détermination d’un modèle de chaque paroi comprenant au moins : Acquérir (101) une image des parois, Extraire (104) des maillages des parois sous la forme d’un nuage de points, Appliquer (105) une approximation par B-splines multi-niveau au nuage de points de sorte à obtenir une surface paramétrée 3D lissée Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes dans lequel les milieux séparés par des parois sont anisotropes et l’étape (111) de déterminer un trajet entre le point source et le point d’intérêt est réalisée en prenant en compte l’anisotropie de chaque milieu traversé. Méthode de détermination d’un champ d’une onde ultrasonore selon la revendication 5 dans laquelle l’étape (111) de déterminer un trajet entre le point source et le point d’intérêt est réalisé au moyen d’un algorithme de minimisation (202) de la somme des segments composant le trajet, chaque segment correspondant à un milieu traversé, la somme étant pondérée par les lenteurs de chaque milieu. Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes dans lequel les modes de propagation de l’onde ultrasonore comprennent au moins un mode de propagation impliquant des réflexions de l’onde ultrasonore sur au moins certaines parois et la prise en compte du coefficient de réflexion du trajet issu du point central de chaque surface élémentaire dans le calcul (112) de la réponse impulsionnelle du champ ultrasonore. Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes dans laquelle la décomposition (107) d’une surface d’un élément en surfaces élémentaires est réalisée à l’aide d’un diagramme de Voronoï. Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes dans laquelle le calcul (112) de la réponse impulsionnelle du champ ultrasonore comporte au moins le calcul de l’amplitude du champ et des bornes minimales et maximales des temps de vol correspondants aux trajets calculés. Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes dans laquelle la zone d’intérêt est une région transcrânienne et les parois sont les parois d’un crâne. Méthode de détermination d’un champ d’une onde ultrasonore selon l’une quelconque des revendications précédentes dans laquelle l’étape de déterminer un champ ultrasonore comprend au moins l’étape de convoluer (115) le résultat de la somme avec un signal de référence associé à l’élément ultrasonore ou d’extraire (115) l’amplitude et la phase d’une composante fréquentielle prédéfinie de chaque réponse impulsionnelle pour obtenir un champ ultrasonore transmis par ledit élément. Programme d’ordinateur comportant des instructions pour l'exécution de la méthode selon l’une quelconque des revendications 1 à 11, lorsque le programme est exécuté par un processeur. Support d'enregistrement lisible par un processeur sur lequel est enregistré un programme comportant des instructions pour l'exécution de la méthode selon l’une quelconque des revendications 1 à 11, lorsque le programme est exécuté par un processeur. Simulateur de champ d’une onde ultrasonore comprenant un processeur et une mémoire configurés pour exécuter les étapes de la méthode selon l’une quelconque des revendications 1 à 11 et un dispositif d’affichage pour afficher le champ calculé dans une zone d’intérêt.