DESCRIPTION TITRE : PROCEDE DE COMPENSATION LOCALE DES ABERRATIONS DANS UN MILIEU DYNAMIQUE EN IMAGERIE ULTRASONORE. Domaine technique [001] La présente invention concerne un procédé et un système permettant de construire une image confocale d’un milieu dynamique avec une compensation locale des aberrations. [002] L’invention s’applique avantageusement au domaine de l’imagerie médicale, mais peut être appliquée à tout domaine d’imagerie par ultrasons. Etat de la technique [003] Dans le domaine de l’imagerie acoustique, on cherche à caractériser un milieu inconnu en le sondant de manière active avec des ondes ultrasonores. C’est notamment le principe de l’échographe en imagerie médicale. [004] Cependant, du fait des inhomogénéités de vitesse du son entre les différents tissus du corps humain, les images échographiques souffrent d’aberrations à la fois transverses et axiales qui dégradent leur contraste et altèrent leur résolution. [005] Par ailleurs, en microscopie de localisation ultrasonore, on arrive à imager avec une haute résolution le réseau vasculaire d’un organe en détectant, localisant et trackant des bulles isolées. Toutefois, ce processus de détection et tracking des bulles repose sur une convolution par une fonction d’étalement gaussienne qui n’est valable que dans un cas idéal sans aberrations. [006] La figure 1 illustre un processus conventionnel de focalisation pour produire une image échographique d’un milieu. Le milieu comprend ici malheureusement une couche aberratrice avec une vitesse du son différente de la vitesse du son observée dans le reste du milieu. Il en résulte une distorsion spatiale et une dispersion temporelle (réverbération) du front d’onde acoustique qui conduit à des aberrations transverses et axiales de l’image échographique résultante. Ces phénomènes conduisent à une dégradation de sa résolution et
de son contraste ainsi qu’à l’apparition d’artefacts de réverbérations, particulièrement gênants lors d’un examen médical. [007] Sur le schéma de gauche de cette figure 1, le réseau de transducteurs, placé en vis-à-vis d’un milieu, permet d’insonifier et d’imager le milieu. La méthode conventionnelle consiste à insonifier le milieu à l’aide d’émissions focalisées par une technique dite de formation de voies (ou “beamforming” en langue anglaise). On applique aux signaux émis par chaque transducteur un jeu de retards appropriés τ basés sur un modèle de vitesse homogène c0, afin de faire interférer constructivement les ondes produites par chaque transducteur au point focal visé de position spatiale rin. = (xin, z). En raison des limites physiques de la diffraction, les ultrasons sont émis au travers de l’ouverture de la sonde échographique, concentrés dans une zone souvent nommée « tache focale ». En outre, les ondes traversant la couche aberratrice sont déformées, ce qui provoque une distorsion et un élargissement de la tache focale autour du point focal. Cet effet indésirable est illustré sur le schéma de gauche de la figure 1. [008] Les ondes réfléchies au point focal sont retournées vers le réseau de transducteurs et traversent de nouveau la couche aberratrice, ce qui déforme encore plus le front d’onde réfléchi mesuré par le réseau de transducteur. Un processus de formation de voie appliqué à ce type de signaux donne lieu à une image échographique présentant d’importantes distorsions latérales du fait de la présence de la couche aberratrice. Si cette couche est réverbérante, les échos de réflexion multiple peuvent donner lieu à une distorsion axiale de l’image échographique. Ces différents effets entrainent une perte de résolution et de contraste de l’image échographique. Une distribution hétérogène de vitesse du son dans les tissus traversés impacte donc la qualité de l’image reconstruite. [009] Le schéma de droite de la figure 1 explicite l’objet de l’invention : déterminer le front d’onde à émettre pour focaliser de manière optimale les ondes ultrasonores à la fois spatialement et temporellement vers chaque point du milieu. Des techniques de focalisation adaptative, ou plus récemment d’imagerie matricielle, ont été développées à cet effet. Toutefois, elles reposent sur une invariance de la tache focale sur une zone suffisante afin de pouvoir
combiner de manière intelligente les ondes réfléchies par différents points contigus. Cela permet en effet de s’affranchir du désordre et accéder à la loi d’aberration associée à la zone considérée, communément appelée zone d’isoplanétisme. Ces techniques de focalisation adaptatives sont bien connues mais restent limitées car elles ne permettent de compenser que des aberrations d’ordre relativement faibles associées à des patchs d’isoplanétisme suffisamment grands. Les aberrations d’ordre plus élevé ainsi que les réverbérations varient trop rapidement pour pouvoir être dressées par ces techniques de l’état de l’art. Il en résulte une distorsion spatio-temporelle du front d’onde acoustique qui conduit à des aberrations significatives de l'image échographique, et donc à une dégradation de sa résolution et de son contraste. Ces aberrations peuvent être telles qu’elles compromettent la caractérisation ultrasonore, et notamment dans le cas d’un examen médical. [0010] La présente invention a pour but un nouveau procédé d’imagerie ultrasonore dans lequel les aberrations d’ordre élevé et les réverbérations sont compensées de manière optimale pour chaque point de focalisation. L’objectif est d’obtenir une image échographique de résolution la plus fine possible et de contraste optimal. Exposé de l’invention [0011] On atteint au moins l’un des objectifs précités avec un procédé de construction ultrasonore d’une image confocale d’un milieu dynamique, le procédé comprenant les étapes suivantes : a) acquisition, au moyen d’un réseau de transducteurs, d’une série de matrices de réflexion canoniques Rui(t, #m)=[R(uout,iin,t,#m)] à différents instants, chaque matrice de réflexion canonique est définie entre une base d’émission d’ondes ultrasonores i en entrée et une base de réception u en sortie ; les coefficients de cette matrice de réflexion canonique correspondant aux signaux reçus par les transducteurs et induits par les ondes ultrasonores réfléchies dans le milieu ; t désignant le temps d’écho et #m désignant la mième matrice de réflexion canonique ; b) détermination d’une matrice de réflexion focalisée Rxx(z, #m) pour chaque matrice de réflexion canonique par focalisation en entrée-sortie pour tout point
d’au moins une région du milieu, les coefficients de cette matrice de réflexion focalisée sont obtenus en calculant un champ de pression acoustique entre tous les points de la région de positions latérales xin et xout, situés à une profondeur attendue z pour une vitesse du son supposée c0 ; c) détermination d’une composante dynamique pour chaque coefficient de chaque matrice de réflexion focalisée de sorte à constituer des matrices de réflexion focalisées dynamiques ^^ ^^ (^, #^) d) détermination d’une loi de correction ^(x, z) pour chaque point x et profondeur z du milieu à partir des matrices de réflexion focalisées dynamiques, e) détermination de matrices de réflexion focalisées corrigées ^ ^ ^^ (^, #^) par application de la loi de correction ^(x, z) en tout point du milieu, f) détermination d’un signal confocal dynamique Sc(x,z, #^) de tout point de position spatiale (x, z) à partir des coefficients diagonaux de la matrice de réflexion focalisée corrigée ^ ^ ^^ (z, #^), g) construction d’une image à partir des signaux confocaux dynamiques. [0012] Avec la présente invention, le procédé permet avantageusement de sonder localement le milieu pour obtenir une estimation locale d’une loi de correction des aberrations adaptée pour corriger le processus de formation de voies ultrasonore. Cette correction permet de réduire ou supprimer des aberrations par exemple dues à des variations de vitesse du son dans le milieu ou à réflexions multiples des ondes générées par une ou plusieurs zones d’aberration dans le milieu. [0013] L’étape d’acquisition peut consister en la réalisation des mesures permettant de traiter les matrices de réflexion ou bien à un traitement a posteriori des matrices à partir de données stockées dans un espace mémoire. La variable t représente le temps d’écho associé aux signaux enregistrés pendant la mesure et f la fréquence d’onde ultrasonore pendant la mesure. Pour chaque mesure, on acquière une amplitude et une phase pour chaque pixel. [0014] Ainsi les différents calculs de l’invention peuvent ainsi être réalisés indépendamment de la phase d’acquisition des mesures, en particulier en modifiant divers paramètres de calcul, ce qui permet de procéder à diverses analyses de construction ultrasonore soit en temps réel, soit a posteriori.
[0015] Ces calculs de correction bénéficient d’informations locales extraites grâce à la partie dynamique des échos réfléchis par les tissus. [0016] Alors que les méthodes de correction d’aberrations existantes (focalisation adaptative, imagerie matricielle) reposent sur une hypothèse d’isoplanétisme locale (invariance spatiale de la tache focale) afin de pouvoir moyenner les corrélations entre signaux ultrasonores provenant de plusieurs points de focalisation contigus, le procédé selon l’invention utilise la dynamique des signaux pour obtenir une loi de focalisation indépendante pour chaque pixel. On n’a pas besoin de moyenne spatiale et on peut donc accéder à des aberrations d’ordre très élevé présentant peu ou pas d’isoplanétisme. [0017] L’invention est notamment remarquable par le fait qu’on acquière une série donc plusieurs matrices de réflexion canoniques à des instants différents. Chaque acquisition peut comprendre plusieurs mesures avec application de plusieurs ondes planes différentes pour constituer une matrice. Une acquisition peut être nommée une trame, ou « frame » en anglais, symbolisée par #m. [0018] Le fait d’avoir plusieurs matrices de réflexion canoniques d’un même milieu acquises à des instants différents permet de tenir compte de l’aspect dynamique du milieu. On se sert de cet aspect dynamique, qui pourrait être considéré comme un inconvénient, pour améliorer la précision dans la détermination des lois de correction. [0019] [0020] Selon une mise en œuvre avantageuse de l’invention, l’étape d’acquisition d’une série de matrices de réflexion canoniques Rui(t, #m) peut comprendre l’émission d’une impulsion ultrasonore depuis chaque transducteur du réseau dont la position est repérée par la coordonnée uin ; cette impulsion donne lieu à une onde incidente cylindrique ou sphérique divergente qui est réfléchie par des diffuseurs du milieu ; ces échos réfléchis forment un champ rétrodiffusé qui est enregistré par chacun des transducteurs en fonction du temps ; la matrice de réflexion canonique Ruu(t,#m) exprimée dans la base des transducteurs étant composée d’un ensemble de réponses impulsionnelles R(uout,uin,t,#m) entre transducteurs. [0021]
[0022] Selon une variante, l’étape d’acquisition d’une série de matrices de réflexion canoniques Rui(t, #m) peut comprendre une insonification du milieu avec une série d’ondes planes avec un retard τ’ appliqué sur chaque signal à l’émission pour la formation d’un front d’onde incliné d’un angle θin par rapport au réseau de transducteurs, un champ rétrodiffusé par le milieu, R(uout, θin, t, est mesuré par tous les transducteurs de position uout pour chaque onde plane incidente θin, l’ensemble des réponses formant une matrice de réflexion canonique Ruθ(t, #m)=[ R(uout, θin, t, #m)]. [0023] [0024] Selon encore une variante, l’étape d’acquisition d’une série de matrices de réflexion canoniques Rui(t, #m) peut comprendre une insonification du milieu avec une série d’ondes divergentes. [0025] [0026] Selon l’invention, l’étape de détermination de la matrice de réflexion focalisée Rxx(z,#m) peut comprendre : [0027] - un processus de focalisation en entrée à partir de chaque matrice de réflexion canonique Rui(t,#m) qui utilise un temps de vol à l'aller des ondes entre la base d'émission d’ondes ultrasonores i et un transducteur virtuel en entrée TVin et qui crée une tache focale dite d’entrée autour d’un premier point P1 de position spatiale rin=(xin,z), ladite tache focale d’entrée correspondant au transducteur virtuel d’entrée TVin, et [0028] - un processus de focalisation en sortie à partir de la matrice de réflexion canonique Rui(t,#m) qui utilise un temps de vol de retour des ondes entre un transducteur virtuel de sortie TVout et les transducteurs de la base de réception u et qui crée une tache focale dite de sortie autour d’un deuxième point P2 de position spatiale rout=(xout,z), ladite tache focale de sortie correspondant au transducteur virtuel de sortie TVout. [0029] [0030] A titre d’exemple, la coordonnée xout du deuxième point P2 est situé à une distance de la coordonnée xin du premier point P1 qui est inférieure ou égale à une distance maximale Δxmax qui est fonction d’un nombre d’ondes ultrasonores générées lors de l’acquisition de la matrice de réflexion.
[0031] [0032] En particulier, la matrice de réflexion focalisée peut être déterminée dans le domaine temporel ou dans le domaine fréquentiel. [0033] Si on traite d’aberrations s’apparentant à de simples décalages temporels d’amplitude relativement faibles (inférieures typiquement à la résolution temporelle des signaux ultrasonores), alors on aura intérêt à calculer directement la matrice de réflexion dans le domaine temporel et au seul temps balistique. [0034] Si on traite d’aberrations d’amplitude plus importantes et/ou de réverbérations et/ou de dispersion fréquentielle, alors il est préférable d’adopter une approche polychromatique et calculer la matrice de réflexion focalisée dans le domaine fréquentiel. [0035] Selon une caractéristique avantageuse de l’invention, la composante dynamique peut être déterminée en soustrayant à chaque coefficient de chaque matrice de réflexion focalisée une moyenne glissante sur N matrices de réflexion focalisées. [0036] De plus, la composante dynamique peut être déterminée en appliquant un filtre passe haut ou passe-bande suivant une dimension des nombres #^ des matrices de réflexion. [0037] On utilise ici avantageusement le fait d’avoir acquis plusieurs trames de matrices de réflexion canoniques. [0038] [0039] La composante dynamique peut encore être déterminée en réalisant une décomposition en valeurs singulières des matrices de réflexion focalisées réarrangées dans une matrice bidimensionnelle dont l’une des dimensions est celle des nombres #^ des matrices de réflexion. En d’autres termes, on concatène l’ensemble des matrices de réflexion dans une matrice bidimensionnelle globale dont on fait la SVD. [0040] [0041] Selon l’invention, l’étape de détermination d’une loi de correction ^(^, ^) peut comprendre, pour chaque matrice de réflexion focalisée dynamique ^^ ^^ (^, #^), les étapes suivantes :
[0042] - détermination d’une matrice de réflexion duale Rcx(z, #^) par projection aller directement ou indirectement de la matrice de réflexion focalisée dynamique ^^ ^^ (^, #^) vers une base de correction (c), [0043] - calcul de la loi de correction ^(^, ^) à partir de la matrice de réflexion duale Rcx(z, #^), ladite loi de correction étant une loi de correction, ^ = [^(^, ^)] sur la base de correction (c), [0044] - détermination des matrices de réflexion focalisées corrigées par projection retour des matrices de réflexion duales corrigées ^^ ^^ (z, #^) vers la base focalisée (x), c’est-à-dire l’ensemble des points x à la profondeur z considérée. [0045] [0046] [0047] De préférence, les coefficients ^^ ^^ (^, #^) = [^^ ^ (^, ^, ^, #^)] de la matrice de réflexion duale corrigée ^ ^ ^^ peuvent être déterminés en effectuant un produit terme à terme entre la matrice de réflexion duale Rcx(z, #^) et le conjugué en phase de la loi de correction ^(^, ^), soit : [0048] ^ ^ ^^ = ^^^ ∘ ^∗ [0049] où le symbole * désigne une opération de conjugaison de phase, le symbole ∘ est le produit d’Hadamard, tel que :
[0051] L’étape de calcul de la loi de correction ^(x, z) peut comprendre les étapes suivantes : - construction d’une matrice de corrélation C(x,z) à partir des matrices de réflexion duale Rcx(z, , #^) pour chaque point (x,z) d’un champ de vision, - détermination de la loi de focalisation ^(x,z) pour chaque point du champ de vision en réalisant l’une des opérations suivantes : - décomposition en valeurs propres de la matrice de corrélation C(x,z), la loi de correction ^(x,z) étant le premier vecteur propre ^^ de la matrice de corrélation C(x,z) dans la base de correction (c), - décomposition en valeurs singulières de la matrice de réflexion duale réarrangée de la manière suivante : ^^#(x, z) = [^(^, #^ , x, z)]
la loi de correction ^(x, z) étant égale au premier vecteur singulier de la matrice de réflexion duale ^^#, i.e. ^(^, ^) = ^^ - résolution de l’équation suivante : ^(x, z) = exp(^ arg{^^^(x, z) × ^(x, z)}) de manière itérative par l’expression suivante, qui correspond à un calcul par retournement de phase itératif : ^^^^(x, z) = exp(^ arg{^^^ × ^^(x, z)}) Où × est le produit matriciel, avec ^^ un front d’onde arbitraire, la loi de correction ^(x, z) étant obtenue par : ^(x, z) = ^ ^^ →^ ^ ^^ (x, z). - résolution de l’équation suivante : ^(x, z) = exp(^ arg{^##(x, z) × ^(x, z)}) où × est le produit matriciel, de manière itérative par l’expression suivante : ^^^^(x, z) = exp(^ arg{^##(x, z) × ^^(x, z)}) où × est le produit matriciel, avec ^^ un front d’onde arbitraire, ce qui permet d’obtenir le vecteur W(x, z) suivant : ^(x, z) = ^ ^^ →^ ^ ^^(x, z) la loi de correction ^(x, z)étant obtenue par :
[0052] Avantageusement, la matrice de corrélation C(x,z) peut être déterminée dans la base de correction c et dans le domaine fréquentiel, par le calcul suivant des éléments de la matrice de corrélation C = Ccc :
où * est l’opérateur de conjugaison. c et c’ étant des points de la base de correction c [0053] Cette opération permet de corréler les champs réfléchis dans la base de correction pour chaque source virtuelle en (x,z). Cette corrélation est moyennée sur les différentes trames #^ , c’est à dire les différentes réalisations du speckle, afin de s’affranchir de la réflectivité aléatoire du milieu et ainsi
synthétiser une étoile guide cohérente à partir des différentes réalisations #^ du speckle. [0054] Autrement, la matrice de corrélation C(x,z) peut être déterminée dans la base des nombres #^, par le calcul suivant des éléments de la matrice de correction C = ^##:
* est l’opérateur de conjugaison. c et c’ étant des points de la base de correction c #m et #l désignant les mièmes et lièmes frames de la séquence de matrices de réflexion enregistrées. [0055] De préférence, l’étape de calcul de la loi de correction ^(^, ^) peut être itérée au moins deux fois avec à chaque itération, la projection aller utilise la matrice de réflexion focalisée corrigée ^ ^ ^^ (^, #^) obtenue durant la projection retour de l’itération précédente à la place de la matrice de réflexion focalisée Rxx(z,#^). [0056] L’étape de détermination d’une matrice de réflexion duale Rcx(z, #^) peut également être itérée au moins deux fois avec, à chaque itération, une utilisation d’une base de correction c différente. [0057] [0058] L’étape de détermination d’une matrice de réflexion duale Rcx(z, #^) peut autrement être itérée au moins deux fois avec, à chaque itération, une utilisation d’une projection aller soit vers une base de correction en entrée, soit vers une base de correction en sortie de la matrice de réflexion focalisée dynamique. [0059] [0060] Selon une caractéristique avantageuse de l’invention, la base de correction peut être l’une des bases suivantes : [0061] - une base des ondes planes ou base de Fourier spatiale, [0062] - une base des transducteurs u, [0063] - une base correspondant au lieu supposé d’aberrateurs dans le milieu, [0064] - une base correspondant à un plan déterminé par optimisation.
[0065] Cela permet de réaliser des projections au niveau des transducteurs, des aberrateurs ou autre. [0066] [0067] L’étape de détermination d’une matrice de réflexion duale Rcx(z, #^) peut être réalisée par projection aller de la matrice de réflexion focalisée dynamique ^^ ^^ (^, #^) vers la base de correction (c) en considérant un propagateur modèle décrivant la propagation des ondes depuis la base focalisée (x) vers la base de correction (c). [0068] Selon une caractéristique de l’invention, la matrice de réflexion considérée au départ du processus peut être la matrice de réflexion focalisée large bande, qui peut être obtenue par formation de voies numériques dans le domaine temporel ou dans le domaine de Fourier à partir des matrices fréquentielles:
[0069] Cette matrice de réflexion large bande correspond à une matrice de réflexion focalisée fenêtrée autour du temps balistique. [0070] L’étape de détermination d’une matrice de réflexion duale large bande Rcx(z, #^) peut être réalisée par projection aller de la matrice de réflexion focalisée large bande obtenue à partir de la matrice de réflexion focalisée dynamique ^^ ^^ (^, ^ = 0, #^) vers la base de correction (c) en considérant un propagateur à la seule fréquence centrale ^^ = (^^ + ^^)⁄ 2. La loi d’aberration obtenue ne présente alors pas de dépendance fréquentielle, ^(^, ^, ^, ^) = ^(^, ^, ^). Cette variante est par exemple à considérer si les aberrations sont assimilables à de simples décalages temporels d’amplitude inférieure à la résolution temporelle de la mesure ultrasonore, cette dernière variant comme l’inverse de la bande passante. Au-delà, une approche multi fréquence est à privilégier. [0071] [0072] Si les aberrations sont seulement axiales (variation de la Vitesse du son seulement suivant la profondeur), l’étape de détermination d’une matrice de réflexion duale Rcx(z, #^) n’est pas nécessaire. La loi de correction que l’on recherche est en effet seulement fréquentielle : ^(^, ^, ^, ^) = ^(^, ^, ^). Dans ce
cas, la matrice considérée au départ du processus est le signal confocal de la matrice de réflexion focalisée dynamique ^^ ^^ (^, #^), ce signal confocal étant tel que : [0073] ^(^, ^, ^, #^) = ^(^, ^, ^, ^, #^) [0074] Dans ce cas, la matrice de corrélation pour l’obtention de la loi de correction fréquentielle, ^(x, z) = [Φ(f, x, z)], peut être donnée par : [0075] ^(^, ^′, ^, ^) = ∑#^ ^(^, ^, ^, #^) ^∗(^, ^, ^^ , #^) où le symbole * désigne une opération de conjugaison de phase f et f’ désignent des fréquences de la bande passante du signal ultrasonore, la base de correction étant directement nommée ici (f), soit
, ^(^, ^, ^, #^) étant le signal confocal de la matrice de réflexion focalisée dynamique ^^ ^^ (^, #^). [0076] L’étape de calcul de la loi de correction ^(x, z), peut être réalisée : [0077] - en définissant une matrice de corrélation pour l’obtention de la loi de correction, ^(x, z) = [Φ(f, x, z)], par : [0078] ^(^, ^′, ^, ^) = ∑#^ ^(^, ^, ^, #^) ^∗(^, ^, ^^ , #^) [0079] - et en moyennant la matrice de corrélation sur les nombres #m mais également sur des pixels adjacents apparentant à un même patch
[0081] f et f’ désignent une fréquence du signal confocal dans la bande passante de la sonde [0082] xp est la coordonnée transverse du point central du patch d’isoplanétisme pour quel on cherche à estimer la loi d’aberration ^ [0083] zp est la coordonnée axiale du point central du patch d’isoplanétisme pour quel on cherche à estimer la loi d’aberration ^ [0084] Ω^: étant un patch d’isoplanétisme centré sur le point ^^^, ^^^. [0085] [0086] Selon l’invention, l’étape de détermination de matrices de réflexion focalisées corrigées ^ ^ ^^ (^, #^) peut être itérée au moins deux fois avec à chaque itération, la projection aller utilise la matrice de réflexion focalisée corrigée
^ ^ ^^ (^, #^) obtenue durant la projection retour de l’itération précédente à la place de la matrice de réflexion focalisée Rxx(z,#^). [0087] De préférence, l’étape g) de construction d’une image à partir des signaux confocaux dynamiques peut comprendre une étape de : - détermination d’une intensité de chaque signal confocal dynamique pour construire une image confocale du milieu, ou - détermination d’une intensité de chaque signal confocal dynamique pour construire une image Doppler de puissance en sommant, pour chaque point de position spatiale (x, z), les intensités de plusieurs signaux confocaux, ou - détermination d’une transformée de Fourier du signal confocal dynamique ^^(^, ^, #) suivant la dimension temporelle #^ :
où la fréquence ^ est la variable conjuguée au temps d’acquisition #, pour construire une image doppler directionnelle du milieu, - détermination de la dynamique locale pour chaque point d’image échographique de position spatiale (x, z), en mesurant la fréquence moyenne ^D(x,z) de la transformée de Fourier ^^(^, ^, ^) du signal confocal complexe:
pour construire une carte de la vitesse axiale des diffuseurs en chaque point de l’image. [0088] Selon une caractéristique avantageuse de l’invention, la matrice de réflexion focalisée Rxx(z,#m) peut dépendre d’une fréquence f des signaux ultrasonores, la loi de correction étant déterminée en fonction de cette fréquence f et des coordonnées du plan de correction. [0089] [0090] Selon une caractéristique avantageuse de l’invention, la matrice de réflexion focalisée Rxx(z,#m) peut être intégrée sur toute la bande passante, la loi de correction étant uniquement déterminée en fonction des coordonnées du plan de correction (c). La matrice de réflexion focalisée, ainsi que la dynamique, est considérée au seul temps balistique.
[0091] [0092] Selon un autre aspect de l’invention, il est prévu un système de construction ultrasonore d’une image confocale d’un milieu dynamique, le système comprenant : [0093] - un réseau de transducteurs adaptés pour générer une série d’ondes ultrasonores incidentes dans une zone d’intérêt du milieu, et pour mesurer en fonction du temps les ondes ultrasonores rétrodiffusées par ladite zone d’intérêt ; et [0094] - une unité de calcul reliée au réseau de transducteurs et adaptée pour mettre en œuvre le procédé décrit ci-dessus. [0095] [0096] On prévoit également un produit programme d'ordinateur comprenant des instructions qui, lorsque le programme est exécuté par un ordinateur, conduisent celui-ci à mettre en œuvre les étapes du procédé décrit ci-dessus. [0097] [0098] On prévoit également un support lisible par ordinateur comprenant des instructions qui, lorsqu'elles sont exécutées par un ordinateur, conduisent celui-ci à mettre en œuvre les étapes du procédé décrit ci-dessus. Brève description des dessins [0099] D’autres avantages et particularités de l’invention apparaîtront à la lecture de la description détaillée de mises en œuvre et de modes de réalisation nullement limitatifs, et des dessins annexés suivants. La figure 1 est une vue schématique illustrant la problématique des aberrations dans l’acquisition des images échographiques selon l’art antérieur ; La figure 2 est une vue schématique illustrant un exemple d’un système de construction ultrasonore pour la mise en œuvre du procédé selon la présente invention ; La figure 3 est un diagramme du procédé de construction d’une image ultrasonore selon la présente invention;
La figure 4 illustre plusieurs vues schématiques 2a à 2f illustrant des séquences d’émission/réception utilisées pour l’imagerie et la caractérisation ultrasonore d’un milieu; La figure 5 est une vue schématique présentant le principe de focalisation dans le procédé selon l’invention ; La figure 6 est une vue schématique illustrant l’opération d’extraction de la composante dynamique des signaux ultrasonores pour une expérience d’imagerie trans-crânienne d’un cerveau de mouton endormi ; La figure 7 est une vue schématique illustrant l’opération d’extraction des lois d’aberrations par retournement de phase itératif dans le speckle dynamique; La figure 8 représente plusieurs images illustrant l’effet de l’application de différentes lois d’aberrations transverses locales dans le cerveau du mouton; La figure 9 comportent plusieurs images confocales avant et après correction d’une partie du cerveau du mouton ; La figure 10 comportent deux visualisations de type « power Doppler » montrant le gain en contraste et résolution apporté par une correction ultra- locale des aberrations; La figure 11 comportent deux images de type cartes de fréquences dynamiques locales dans le cerveau du mouton, avant et après correction des aberrations transverses. Description détaillée des figures [00100] Il est bien entendu que les modes de réalisation qui seront décrits dans la suite ne sont nullement limitatifs. On pourra notamment imaginer des variantes de l’invention ne comprenant qu’une sélection de caractéristiques décrites par la suite isolées des autres caractéristiques décrites, si cette sélection de caractéristiques est suffisante pour conférer un avantage technique ou pour différencier l’invention par rapport à l’état de la technique antérieure. Cette sélection comprend au moins une caractéristique de préférence fonctionnelle sans détail structurel, ou avec seulement une partie des détails structurels si cette partie uniquement est suffisante pour conférer un avantage
technique ou pour différencier l’invention par rapport à l’état de la technique antérieure. [00101] Les différents modes de réalisation et aspects décrits dans la présente divulgation peuvent être combinés ou simplifiés de multiples manières. En particulier, les étapes des différents procédés peuvent être répétées, interverties, et/ou exécutées en parallèle, sauf précision contraire. [00102] La présente divulgation concerne des procédés et systèmes de caractérisation ultrasonore d'un milieu, et s'applique notamment à l'imagerie médicale de tissus vivants ou non. Le milieu est par exemple un milieu hétérogène que l'on cherche à caractériser pour par exemple identifier et/ou caractériser les hétérogénéités. Ces techniques de construction sont notoirement non invasives pour le milieu, qui est avantageusement préservé en particulier dans sa nature et son intégrité. [00103] [00104] Approche usuelle de l’imagerie ultrasonore [00105] Dans le domaine de l’imagerie par ultrasons, on cherche souvent à construire une image de la réflectivité d’un milieu à partir des échos rétrodiffusés par des hétérogénéités du milieu. C’est le principe de l’échographe utilisé en imagerie médicale, qui permet notamment de visualiser l’anatomie interne d’un individu ou d’un animal. A des fins de simplification, pour permettre la construction d’une image échographique, le milieu est considéré comme homogène, avec une vitesse de propagation du son c0 constante. [00106] La dynamique du milieu peut être naturelle (par exemple liée au passage du sang dans les vaisseaux sanguins et au cycle cardiaque) ou peut être induite par une vibration du milieu induit par pression de radiation ("push" acoustique) ou par vibration d’une sonde. [00107] Les méthodes d’échographie conventionnelle utilisent généralement un réseau de transducteurs piézo-électriques qui peuvent émettre et/ou recevoir des signaux ultrasonores de manière indépendante ou quasi indépendante, chaque transducteur étant à une position u dans la barrette supportant ledit réseau. Le réseau de transducteurs, placé en vis-à-vis d’un
milieu, permet d’insonifier et de construire une image représentative du milieu de différentes manières. Une méthode conventionnelle consiste à insonifier le milieu à l’aide d’émissions focalisées par une technique dite de formation de voies (ou “beamforming” en langue anglaise). Cette méthode consiste à appliquer aux signaux émis par chaque transducteur un jeu de retards appropriés τ(uin, xin, z, c0) basés sur un modèle de vitesse homogène c0, afin de faire interférer constructivement les ondelettes produites par chaque transducteur au point focal visé de position spatiale (xin, z). En raison des limites physiques de la diffraction, les ultrasons sont émis au travers de l’ouverture de la sonde échographique, concentrés dans une zone souvent nommée « tache focale », de largeur latérale δx. [00108] Afin de permettre de construire ensuite éventuellement une image échographique illustrant les caractéristiques du milieu étudié, une étape de focalisation numérique est également effectuée en réception. Les échos captés par les transducteurs du réseau, sont remis en phase en les décalant temporellement. Les délais τ(uout, xout, z, c0) sont identiques à ceux appliqués à l’émission, la variable uout désignant la position de chaque transducteur. Dans la phase d’émission, tous les signaux interfèrent au point de position (xin, z) au temps balistique t = z/c0 si le modèle de vitesse c0 utilisé correspond à la réalité du milieu étudié. En réception, les signaux provenant de ce même point (xout = xin) interfèrent par sommation des signaux au temps d’écho t = 2z/c0. Cette sommation permet d’obtenir le résultat final de la focalisation en réception. Cette méthode confocale à double focalisation à l’émission et à la réception permet d’imager directement la réflectivité du milieu avec une résolution latérale δx et un bon contraste. Toutefois, cette méthode est chronophage car elle nécessite de focaliser physiquement à l’émission en chacun des points du milieu ou au moins à une profondeur donnée, sur chacune des lignes de l’image construite représentative du milieu. [00109] [00110] La présente invention a pour objet de perfectionner les procédés de sondage par ultrasons connus, notamment afin de corriger des aberrations. [00111]
[00112] [00113] Système de construction ultrasonore [00114] La figure 2 illustre un exemple d’un système 1 d’imagerie ultrasonore pour la mise en œuvre du procédé d’imagerie ultrasonore d’un milieu tel qu'un milieu hétérogène M, selon la présente invention. Ce système et le procédé permettent la formation d’une image échographique par ultrasons d’au moins une partie (zone d’intérêt ou champ de vision) du milieu. [00115] Le système 1 comprend : [00116] - un dispositif de sondage 20 ou sonde 20, [00117] - une unité de calcul 30 pour calculer une image à partir des signaux reçus de la sonde 20, [00118] - un panneau de contrôle 40 relié à l’unité de calcul 30, ce panneau de contrôle comprenant par exemple des boutons 41 et un pavé tactile 42, [00119] - un dispositif d’affichage 50 pour visualiser une image et divers éléments ou mesures. [00120] La sonde 20 est reliée à l’unité de calcul 30 via un câble 21 ou via une connexion sans fil, et est capable d’émettre des ondes ultrasonores W dans le milieu M et de recevoir des ondes ultrasonores W depuis le milieu M, lesdites ondes ultrasonores étant résultantes de réflexions des ondes ultrasonores émises sur des particules diffusantes ou diffuseurs à l’intérieur du milieu. [00121] La sonde 20 peut comprendre un réseau 10 comprenant une pluralité de transducteurs 11. Le réseau 10 est par exemple un réseau linéaire ou courbé ou bidimensionnel ou matriciel. Les transducteurs 11 sont capables de convertir un signal électrique en une vibration et réciproquement. Les transducteurs 11 sont par exemple des transducteurs piézoélectriques ultrasonores pouvant se présenter sous la forme d'une barrette rigide mise en contact directement ou indirectement avec une surface externe du milieu M pour être couplé au milieu et pour faire vibrer et émettre et recevoir des ondes ultrasonores W. Le réseau 10 de transducteurs 11 de la sonde 20 est alors associé à l’unité de calcul 30. Le réseau 10 de transducteurs peut comprendre une centaine ou plus de transducteurs 11.
[00122] L’unité de calcul 30 peut comprendre un boitier 31 incluant des dispositifs de réception pour amplifier et/ou filtrer les signaux reçus de la sonde 20, et des convertisseurs (convertisseurs analogique vers digital, et convertisseurs digital vers analogique) pour transformer les signaux en données représentatives du signal. Les données peuvent être enregistrées dans une mémoire de l’unité de calcul 30 et/ou directement traités pour calculer des données intermédiaires (données de formation de voie ou autres). L’unité de calcul 30 peut implémenter tout procédé connu permettant de construire une image à partir des données des signaux reçus de la sonde 20, tel que la formation de voies. [00123] L’image calculée peut être : [00124] - une image du milieu (image B-mode) habituellement en niveau de gris pour visualiser des organes dans le milieu, et/ou [00125] - une image montrant une vitesse ou un flux dans le milieu (image couleur) par exemple utile pour visualiser des vaisseaux sanguins dans le milieu, et/ou [00126] - une image montrant une caractéristique mécanique du milieu (élasticité) par exemple utile pour identifier des tumeurs à l’intérieur du milieu. [00127] Par "connexion" ou "liaison" entre le dispositif de sondage 20, l'unité de calcul 30 et le dispositif d'affichage 50, on entend tout type de liaison filaire de type électrique ou optique, ou tout type de liaison sans fil utilisant tout protocole tel que le WiFiTM, BluetoothTM ou autres. Ces connexions ou liaisons sont à simple sens ou double sens. Le dispositif d’affichage 50 associé peut être de tout type, tel qu’un écran tactile ou non tactile, connecté ou pas. [00128] Le dispositif d’affichage 50 est un écran permettant de visualiser l’image calculée par l’unité de calcul 30. Le dispositif d’affichage 50 peut aussi visualiser d’autres informations telles que les échelles de l’image, ou des informations de configuration pour le calcul ou traitement ou toute information de mesure ou d’aide. L’écran 50 peut être articulé sur un bras support 51 pour un meilleur positionnement pour l’utilisateur. L’écran 50 est usuellement un écran de grande taille (au moins 20 pouces) pour une meilleure visualisation pour l’utilisateur.
[00129] Le panneau de contrôle 40 est par exemple une portion d’un boitier système, ladite portion comprenant un boitier de panneau disposant d’une surface sensiblement plane 40a inclinée vers l’utilisateur pour manipulation d’une seule main. Comme représenté en figure 2, le panneau de contrôle 40 peut comprendre un écran de contrôle 49 pour visualiser diverses informations de configuration. [00130] L’unité de calcul 30 est configurée pour la mise en œuvre d’étapes de calculs et/ou traitement, notamment pour la mise en œuvre d’étapes de procédés selon la présente divulgation. Par convention, comme représenté en figure 5 montrant un réseau 10 de transducteurs 11 sur une surface d’un milieu M, on définit un repère spatial du milieu M, en prenant un premier axe X et un deuxième axe Z perpendiculaire à celui-ci. Par simplification, le premier axe X correspond à la direction transversale dans laquelle les transducteurs 11 sont alignés dans l’exemple d’un réseau linéaire, et le deuxième axe Z correspond à la profondeur du milieu M par rapport à ce réseau 10 de transducteurs 11. Cette définition peut être adaptée au contexte et ainsi par exemple étendue à un repère spatial à trois axes dans le cas d’un réseau 10 bidimensionnel, ou à un repère polaire dans le cas d’un réseau 10 courbé, ou à tout autre repère adapté et/ou dépendant de la structure et forme du réseau 10 de transducteurs ultrasonores. Ainsi, dans la suite de la présente divulgation, nous utiliserons un repère cartésien XZ, correspondant à une sonde 20 linéaire, pour plus de simplicité dans les explications, mais un spécialiste du domaine généraliserait et appliquerait facilement les résultats à tout type de repère. [00131] Dans la suite de la divulgation, il est fait référence à un réseau 10 de transducteurs 11 pour l’émission et la réception, étant bien entendu que, dans un cas plus général, plusieurs réseaux de transducteurs pourront être utilisés simultanément. Les transducteurs 11 pourront être à la fois émetteur puis récepteur, ou bien seulement émetteur pour certains et seulement récepteur pour d’autres. De même, un réseau 10 peut être constitué d’un (1) à N transducteurs 11, de type identique ou de natures différentes. [00132] Le réseau 10 de transducteurs 11 sert par exemple à la fois comme émetteur et comme récepteur, ou est constitué de plusieurs sous-réseaux de
transducteurs, certains étant dédiés à l’émission, d’autres à la réception des ondes ultrasonores. Par réseau de transducteurs, on entend au moins un transducteur, une suite alignée ou non de transducteurs, ou une répartition bidimensionnelle de transducteurs (par exemple une matrice de transducteurs), ou toute répartition spatiale de transducteurs. [00133] Lorsque dans la présente divulgation, il est fait référence à des étapes de calcul ou traitement pour la mise en œuvre notamment d’étapes de procédés, il est entendu que chaque étape de calcul ou traitement peut être mis en œuvre par logiciel, hardware, firmware, microcode ou toute combinaison appropriée de ces technologies ou technologies avoisinantes. Lorsqu’un logiciel est utilisé, chaque étape de calcul ou traitement peut être mise en œuvre par des instructions de programme d’ordinateur ou du code qui peut être par exemple interprété, ou exécuté. Ces instructions peuvent être stockées ou transmises vers un support de stockage lisible par un ordinateur (ou unité de calcul) et/ou être exécutées par un ordinateur (ou unité de calcul) afin de mettre en œuvre ces étapes de calcul ou traitement. [00134] [00135] Sur la figure 3 est représenté un diagramme des principales étapes selon l’invention. On distingue une étape a) d’acquisition d’une série de matrices de réflexion canoniques Rui(t, #m). En b), on détermine un ensemble de matrices de réflexion focalisées Rxx(z, #m) du milieu par un processus de focalisation pour plusieurs points de position spatiale transverse (x, y) d’une région et de position axiale z=c0t/2 à partir de la matrice de réflexion canonique Rui(t, #m) pour un modèle de vitesse du son c0, [00136] A l’étape c) on construit des matrices de réflexion focalisées dynamiques ^^ ^^ (^, #^) en supprimant la composante statique. [00137] L’étape d) permet de calculer des lois de correction ^(x, z) pour chaque pixel de l’image à partir des matrices de réflexion duales, ces dernières étant obtenues par projection en sortie des matrices de réflexion focalisées dynamiques ^^ ^^ (^, #^) dans une base de correction c des aberrations, [00138] A l’étape e), on détermine des matrices de réflexion focalisées corrigées ^ ^ ^^ (^, #^).
[00139] A l’étape f), on détermine un signal confocal dynamique Sc(x,z, #^) pour tout point de position spatiale (x, z). [00140] Puis construction à l’étape g) d’une image à partir des signaux confocaux dynamiques. [00141] Ces étapes sont décrites plus en détails ci-après. [00142] [00143] Selon l’invention, le procédé de construction ultrasonore mis en œuvre par l'unité de calcul 30 du système 1 comprend une série d’acquisitions de matrices de réflexion à différents instants. Chaque réalisation de la matrice de réflexion sera appelé trame ou « frame » par la suite et la mième frame sera noté #m. Chaque matrice de réflexion peut être acquise de la manière suivante : [00144] - une étape de génération d’une série d’ondes ultrasonores incidentes USin dans une zone dudit milieu, au moyen d’un réseau 10 de transducteurs 11, ladite série d’ondes ultrasonores incidentes étant une base d’émission i ; et [00145] - pour chaque onde émise iin, le champ réfléchi par le milieu est mesuré par chaque transducteur et est noté R(uout,iin,t,#m), où t est le temps d’écho et le vecteur uout repère la position de chaque transducteur. Chaque champ est stocké dans une série de matrices de réflexion canonique Rui(t, #m)=[R(uout,iin,t,#m)] définie entre la base d’émission i en entrée et une base de réception u en sortie . [00146] Une première possibilité pour mesurer cette matrice de réflexion canonique est d’émettre successivement une impulsion ultrasonore depuis chaque transducteur du réseau dont la position est repérée par la coordonnée uin, comme cela est schématisé en figure 4(a). Cela donne lieu à une onde incidente cylindrique (ou sphérique) divergente. Cette onde est réfléchie par les diffuseurs du milieu et le champ rétrodiffusé est enregistré par chacun des transducteurs en fonction du temps comme représenté en figure 4(b). En répétant cette opération avec chaque transducteur utilisé successivement comme source, on détermine la matrice de réflexion canonique Ruu(t,#m) exprimée dans la base des transducteurs, composée de l’ensemble des réponses impulsionnelles R(uout,uin,t,#m) entre chaque transducteur. Cette matrice est alors riche de quantité d’information sur le milieu étudié. Toutefois, la méthode
suppose que le milieu reste fixe pendant toute la durée des mesures. En outre, les signaux enregistrés possèdent un mauvais rapport signal à bruit car le milieu est insonifié par un seul transducteur. [00147] Une deuxième manière de construire cette matrice de réflexion canonique consiste à insonifier le milieu avec une base de série d’ondes planes. Cette méthode permet de s’affranchir des problèmes précédents. La figure 4(c) illustre le principe de cette illumination en ondes planes. Une loi de retard τ’ est appliquée sur chaque signal à l’émission pour la formation d’un front d’onde incliné d’un angle θin par rapport au réseau de transducteurs. A la réception, illustré en figure 4(d), le champ rétrodiffusé par le milieu, R(uout, θin, t, #m), est mesuré par tous les capteurs de position uout pour chaque onde plane incidente θin. L’ensemble de ces réponses forment une matrice de réflexion canonique Ruθ(t, #m)=[ R(uout, θin, t, #m)]. Cette méthode a donné naissance à l’imagerie ultrarapide, et à l’élastographie, et elle est par exemple décrite dans le document : [00148] « Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography », G. Montaldo et al. (IEEE Trans. Ultrason., Ferroelect. Freq. Control 56 489-506, 2009). [00149] Une troisième manière pour créer cette matrice de réflexion canonique consiste à insonifier le milieu avec une base des ondes divergentes, comme représenté en figure 4(e) et figure 4(f), ce qui permet d’illuminer le champ acoustique de manière plus large que par l’utilisation des ondes planes. Cette base est repérée par la position sin de la source virtuelle associée à chaque onde divergente. Cette technique utilisée notamment en imagerie de super- résolution, est explicitée dans le document : [00150] « Ultrafast imaging of the heart using Circular Wave Synthetic Imaging with Phased Arrays», Couade et al., IEEE International Ultrasonics Symposium (2009). [00151] Chaque matrice de réflexion canonique Rui(t,#m) enregistrée peut être une matrice « réelle », c’est-à-dire composée de coefficients réels dans le domaine temporel, les signaux électriques enregistrés par chacun des transducteurs étant des nombres réels. En variante, cette matrice peut être une
matrice « complexe », c’est-à-dire composée de valeurs complexes, par exemple dans le cas d’une démodulation pour une formation de voies en phase et en quadrature (connu en langue anglaise sous la dénomination « beamforming IQ »). [00152] [00153] Focalisation de la matrice de réflexion [00154] Après la séquence d’acquisition ou en parallèle de celle-ci, un processus de formation de voies est appliqué de manière indépendante en entrée et en sortie des matrices de réflexion mesurées. Le résultat est une série de matrices de réflexion focalisées Rxx(z,#m) qui comprend des réponses R(xout, xin,z,#m) du milieu entre un transducteur virtuel d’entrée TVin de position spatiale (xin,z) et un transducteur virtuel de sortie TVout de position spatiale (xout,z). [00155] Les réponses de la matrice de réflexion focalisée Rxx(z,#m) correspondent à un champ de pression acoustique calculé entre tous les points du milieu de positions latérales xin et xout, situés à la profondeur attendue z=c0t/2 et à un temps d’écho t, et pour une vitesse du son supposée c0. [00156] Dans l’étape de détermination de la matrice de réflexion focalisée Rxx(z,#m), on applique en référence à la figure 5 : [00157] - un processus de focalisation en entrée à partir de chaque matrice de réflexion canonique Rui(t,#m) qui utilise un temps de vol à l'aller des ondes entre la base d'émission i et le transducteur virtuel en entrée TVin et qui crée une tache focale dite d’entrée autour du premier point P1 de position spatiale rin = (xin,z), ladite tache focale d’entrée correspondant au transducteur virtuel d’entrée TVin, [00158] - un processus de focalisation en sortie à partir de la matrice de réflexion canonique Rui(t,#m) qui utilise un temps de vol de retour des ondes entre le transducteur virtuel de sortie TVout et les transducteurs de la base de réception u et qui crée une tache focale dite de sortie autour du deuxième point P2 de position spatiale rout=(xout,z), ladite tache focale de sortie correspondant au transducteur virtuel de sortie TVout.
[00159] Ces processus de focalisation en entrée et en sortie forment un processus de focalisation en entrée-sortie, dénommé dans la suite de cette divulgation, processus de focalisation ou plus simplement focalisation. [00160] Autrement dit, dans ce procédé de construction ultrasonore, le transducteur virtuel d'entrée TVin correspond à une "source virtuelle" ultrasonore située à la position spatiale rin dans le milieu et le transducteur virtuel de sortie TVout correspond à un "capteur virtuel" ultrasonore situé à la position spatiale rout. Cette source et ce capteur virtuels sont spatialement séparés de la différence de leurs positions spatiales Δx = xout - xin. [00161] La distance maximale Δxmax entre xout et xin est dictée par le nombre d’illuminations utilisée pour acquérir la matrice de réflexion. Dans la base des ondes planes, par exemple, l’échantillonnage angulaire ^^ de la séquence d’illuminations impose Δxmax~^/(2^^) pour éviter que la matrice de réflexion focalisée soit polluée par des lobes de réseau. Dans l’idéal, Δxmax est choisi en fonction du niveau d’aberrations afin de limiter le nombre de données à acquérir et conserver en mémoire. Il est typiquement de l’ordre de quelques mm en imagerie ultrasonore. [00162] La profondeur attendue de transducteurs virtuels est le paramètre z utilisé dans la loi de focalisation pour un modèle de vitesse du son c0. Leur profondeur réelle est dictée par la position axiale (en profondeur) du volume isochrone, c’est à dire par le temps d’écho t et par la distribution de vitesse du son c(r) dans le milieu. La dimension latérale des transducteurs virtuels est dictée par la tache focale produite par la focalisation à cette profondeur réelle. [00163] Chaque matrice de réflexion focalisée Rxx(z,#m), peut être déterminée ou calculée : [00164] - soit dans le domaine temporel, auquel cas elle peut être notée explicitement avec le paramètre de temps t, c’est-à-dire notée Rxx(z,t,#m), et [00165] en pratique les données des matrice de réflexion focalisée Rxx(z,t,#m) sont calculées entre deux instants temporels prédéterminés ; et [00166] en pratique, elle peut être seulement considérée au temps balistique attendu (i.e à t=0 dans la base focalisée)
[00167] - soit dans le domaine fréquentiel, auquel cas elle peut être notée explicitement avec le paramètre de pulsation ^ qui correspond à une fréquence f par ^ = 2^^, c’est-à-dire notée Rxx(z, ^,#m), et [00168] en pratique les données de la matrice de réflexion focalisée Rxx(z, ^,#m) sont calculées entre deux pulsations ^, d’une bande passante fréquentielle, par exemple entre une pulsation inférieure ^^ et une pulsation supérieure ^^ , pour une pulsation centrale ^^. [00169] Ainsi, la suite des calculs du procédé peut être effectué dans le domaine temporel ou dans le domaine fréquentiel. [00170] [00171] Dans le premier cas de calcul dans le domaine temporel, chaque matrice de réflexion focalisée Rxx
, du milieu entre le transducteur virtuel d’entrée TVin et le transducteur virtuel de sortie TVout est obtenue par focalisation par un calcul de formation de voies en entrées et en sortie. Les coefficients de cette matrice de réflexion focalisée Rxx(z,t,#m), peuvent être déterminés par :
[00173] dans laquelle : [00174] Nin est un premier coefficient de normalisation, [00175] Nout est un deuxième coefficient de normalisation, [00176] Rui(t, #^) est la matrice de réflexion, dont chaque coefficient R(uout, iin, t, #^) est le champ enregistré par le transducteur de position spatiale uou consécutif à l’émission d’indice iin dans la base d’émission (i) et au temps t, auquel on a appliqué les temps de retard τin et τout ; [00177] ^^^(^^^, ^^^^) et ^^^^(^^^^ , ^^^^ , ^) sont des coefficients d’apodisation qui sont prédéfinis, par exemple pour garder une ouverture numérique constante aussi bien à l’émission qu’à la réception ;
[00178] Le premier coefficient de normalisation Nin peut par exemple être défini par : ^^^(^^^, ^) = ∑^^^ ^^^(^^^, ^^^, ^) De même, le deuxième coefficient de normalisation Nout peut être défini par : ^^^^(^^^^ , ^) = ∑^^^^ ^^^^(^^^^ , ^^^^ , ^) [00179] ^^^(^^^, ^^^, ^) est le temps de vol attendu pour chaque onde incidente ^^^ pour atteindre le premier point focal de position spatiale (^^^ , ^) dans un milieu modèle de vitesse du son c0 [00180] ^^^^(^^^^ , ^^^^ , ^) est le temps de vol attendu pour une onde réfléchie depuis le deuxième point focal de position spatiale (^^^^ , ^) au transducteur de position ^^^^. [00181] Ces temps de retard ^^^ et ^^^^ sont habituellement calculés par un homme du métier à partir d’un modèle établi de vitesse du son. Une hypothèse relativement simplificatrice consiste à faire l’hypothèse d’un milieu homogène avec une vitesse du son constante c0 dans le milieu. Dans ce cas, les temps de vols sont directement obtenus à partir des distances entre les transducteurs de la sonde et les transducteurs virtuels. Ainsi, ces calculs de temps de retard sont fonction du type d’onde, de la vitesse du son supposée, et de la géométrie du réseau de transducteurs. [00182] Par exemple, dans le cas particulier d’une onde plane avec un angle d’émission ^ in : [00183] - le temps de retard à l’émission ^ in peut être obtenu par :
[00185] - le temps de retard à la réception ^ out peut être obtenu par :
[00187] Ainsi, ces exemples de calculs de temps de retard montrent bien qu’ils sont fonction du type d’onde et de la vitesse du son, supposée ici constante dans le milieu. [00188] [00189] Le nombre d'éléments de la base d'émission Nin est par exemple supérieur ou égal à un (1), et avantageusement supérieur ou égal à deux (2). Le nombre d'éléments de la base de réception Nout est par exemple supérieur ou égal à deux (2).
[00190] [00191] Enfin, chaque matrice de réflexion focalisée Rxx(z,t,#m) exprimée dans le domaine temporel, peut être transformée dans le domaine fréquentiel en une matrice de réflexion focalisée Rxx(z, ^, #m) par une transformée de Fourier, c’est-à-dire par :
[00193] Cette transformée de Fourier peut être implémentée par tout type de transformée de Fourier discrète, normalisée ou non. [00194] [00195] Dans le deuxième cas de calcul dans le domaine fréquentiel, chaque matrice de réflexion canonique Rui(t, #m) exprimée dans le domaine temporel, puisqu’elle est constituée des signaux reçus par les transducteurs, peut être transformée dans le domaine fréquentiel en une matrice de réflexion canonique #^)) par une transformée de Fourier, c’est-à-dire par
[00197] [00198] Cette transformée de Fourier peut être implémentée par tout type de transformée de Fourier discrète, normalisée ou non. [00199] [00200] Ainsi, la matrice de réflexion focalisée Rxx(z, ^, #m) du milieu peut être obtenue par focalisation par le calcul matriciel ci-dessous, sensiblement équivalent à la focalisation par la formation de voie temporelle explicitée précédemment, c’est-à-dire par le produit matriciel suivant :
[00202] dans lequel : [00203] la matrice ^^^(^, #^) est la transformée de Fourier de chaque matrice de réflexion canonique ^^^(^, #^), [00204] la matrice ^^^(^, ^) est la matrice de passage de réception adaptée pour le passage de la base de réception (u) à la base focalisée (x) à la profondeur z et à la pulsation ^,
[00205] la matrice ^^^(^, ^) est la matrice de passage d’émission adaptée pour le passage de la base d’émission (i) à la base focalisée (x) à la profondeur z et à la pulsation ^, [00206] Les symboles ∗ et † désignent respectivement les opérations matricielles de conjugaison et de transposition-conjugaison. [00207] Le symbole × désigne un produit matriciel. [00208] [00209] Extraction de la composante dynamique [00210] Afin d’extraire les signaux des diffuseurs en mouvement, une étape de filtrage de la composante statique de la matrice de réflexion est effectuée. La séparation des composantes dynamique et statique peut être réalisée de diverses manières. [00211] Suivant un premier mode de réalisation, la composante dynamique est isolée en soustrayant au signal une moyenne glissante sur N frames (typiquement N=5). Les coefficients de la matrice de réflexion focalisée dynamique R sont ainsi obtenues : [00212] ^^(^, ^, ^, ^, #^) = ^(^, ^, ^, ^, #^) − 〈^(^, ^, ^, ^, #^)〉^^^^^^^^^ [00213] où ^^ ^^ (^, ^, #^) = [^^(^, ^, ^, ^, #^)] est la composante dynamique de la matrice de réflexion. Notez que le même filtre peut être appliqué en amont à la matrice mesurée Rui(t, #m), ce qui serait computationnellement plus avantageux. Toutefois, l’intérêt de l’appliquer sur la matrice de réflexion focalisée est de pouvoir adapter le nombre N de frames sur lequel la moyenne glissante est réalisée. [00214] Suivant un deuxième mode de réalisation, un filtre passe haut (ou passe bande) plus sophistiqué est appliqué suivant la dimension #^ de la matrice de réflexion afin de filtrer la composante statique de la matrice de réflexion. [00215] Suivant un troisième mode de réalisation, la composante dynamique est isolée en réalisant une décomposition en valeurs singulières de la matrice de réflexion focalisée réarrangée sous une forme bidimensionnelle de la manière suivante, ^# = [^({^, ^, ^, ^}, #^)]. La décomposition en valeurs singulières de cette matrice s’écrit :
[00216] ^ = ∑ ^ ^ # ^ ^^^ ^^ [00217] avec ^^ = ^^^({^, ^, ^, ^}, ^)^ correspondant aux vecteurs singuliers de la matrice de réflexion dans l’espace associé {^, ^, ^, ^}. [00218] ^^ = ^^^(#^)^ correspondant aux vecteurs singuliers de la matrice de réflexion focalisée ^# dans la base des frames, [00219] ^^ correspondant aux valeurs singulières réelles et positives de la matrice de réflexion focalisée ^# rangées dans un ordre décroissant : ^^ > ^^ > ⋯ > ^^ [00220] Les P premiers espaces propres associés aux valeurs singulières les plus hautes (^^^^^^^) sont associées à la composante statique. Les Q derniers espaces propres associés aux valeurs singulières les plus faibles (^^ < ^^^^) sont associées au bruit. Les espaces propres associés aux valeurs singulières intermédiaires (^^^^< ^^ < ^^^^) sont associées au signal dynamique d’intérêt :
[00222] [00223] Les seuils ^^^^ et ^^^^ peuvent être déterminés comme étant des points d’inflexion de la distribution des valeurs singulières. [00224] [00225] A des fins d’allégement des notations, la composante dynamique de la matrice de réflexion sera notée R et non ^^ dans la suite. [00226] La figure 6 illustre l’extraction de la composante dynamique des signaux à partir de données acquises sur un cerveau de mouton. Les différentes représentations de cette figure 6 permet de comparer l’évolution de la réflectivité d’un point d’une zone de speckle du milieu au fil du temps avant filtrage, sur l’image A, et à celle du même point après filtrage, image B. L’évolution temporelle de la réflectivité complexe associée au pixel désigné par la croix blanche [A et B] est représentée à l’aide d’un nuage de points dans le plan complexe [C] avant et après filtrage de la partie dynamique des signaux, respectivement en périphérie et au centre. L’image [D] est un agrandissement de la partie centrale du plan complexe permettant de mettre en évidence le caractère quasi aléatoire de la réflectivité complexe après filtrage.
[00227] [00228] Procédé de correction des aberrations [00229] Le procédé a pour but la correction des aberrations, ces aberrations étant par exemple dues à des variations de structures dans le milieu qui induisent des variations de vitesse du son et des variations de réflectivité. [00230] Le procédé comprend un traitement de correction comprenant des étapes de : [00231] - détermination d’un ensemble de lois de focalisation ^ à partir des réponses du milieu obtenues pour les différentes réalisations du speckle mesurées aux instants #^. [00232] - détermination de matrices de réflexion corrigées ^′^^(^, ^, #^) du milieu par application de la loi de correction ^ aux matrices de réflexion mesurées ^^^(^, ^, #^). [00233] Grâce à ces dispositions, le procédé permet avantageusement de sonder localement le milieu et de corriger la matrice de réflexion focalisée par rapport aux aberrations, notamment en déterminant une loi de correction pour chaque point du milieu et, optionnellement, pour chaque fréquence f de l’onde ultrasonore. [00234] En outre, le procédé peut comprendre une étape de : [00235] - détermination d’une intensité ^^ d’un point d’image échographique correspondant à un point de position spatiale r = (x, z) à partir des diagonales des matrices de réflexion corrigées ^′^^(^, ^, #^) [00236] - détermination d’une intensité d’un point d’une image power Doppler en sommant l’intensité de l’image échographique au même point sur l’ensemble ou une partie des frames #^ [00237] - détermination d’une image Doppler colorée permettant d’accéder à la directionnalité du mouvement dans les tissus en effectuant une transformée de Fourier suivant la dimension temporelle #^ de l’amplitude complexe de l’image échographique comme décrit dans l’article : [00238] E. Mace, G. Montaldo, B. -F. Osmanski, I. Cohen, M. Fink and M. Tanter, "Functional ultrasound imaging of the brain: theory and basic
principles," in IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 60, no. 3, pp. 492-506, March 2013, [00239] [00240] [00241] Détermination de la loi de focalisation
[00242] L’étape de détermination d’une loi de correction ^, comprend alors des sous-étapes effectuées à chaque profondeur z et chaque fréquence f, de : [00243] - détermination d’un ensemble de matrices de réflexion duale Rcx(z,f, #^) par projection aller des matrice de réflexion focalisée Rxx(z,f,#^) vers une base de correction (c), [00244] - calcul de la loi de correction ^(x, z) à partir de la matrice de réflexion duale Rcx(z,f, #^), ladite loi de correction étant déterminée sur la base de correction (c), ^ = [^(^, ^, ^)], de sorte que ladite loi de correction ^(x, z) est une loi de correction spatio-fréquentielle, [00245] - détermination d’une matrice de réflexion duale corrigée ^ ^ ^^ autour du point de référence et dont les coefficients s’écrivent selon ^ ^ ^^ (^, ^, #^) = [ ^^ ^ (^, ^, ^, ^, #^)], déterminée en effectuant le produit terme à terme entre la matrice de réflexion duale Rcx(z, f, #^) et le conjugué en phase de la loi de correction ^(x, z), c’est-à-dire par : [00246] ^ ^ ^^ = ^^^ ∘ ^ ∗ [00247] où [00248] le symbole * désigne une opération de conjugaison de phase [00249] le symbole ∘ est le produit d’Hadamard, tel que : [00250] ^^^ ^ (^, ^, ^, ^, #^) = ^^^(^, ^, ^, ^, #^)^∗(^, ^, ^, ^). [00251] L’étape finale du procédé comprend la détermination d’un ensemble de matrices de réflexion focalisées corrigées ^ ^ ^^ (^, #^) par projection retour des matrices de réflexion duales corrigées ^ ^ ^^ (z,f, #^) vers la base focalisée (x). [00252] Grâce à ces dispositions, le procédé permet avantageusement de sonder localement le milieu et de corriger la matrice de réflexion focalisée par rapport aux aberrations, notamment en déterminant une loi de focalisation pour chaque point du milieu, et pour chaque fréquence f de l’onde ultrasonore.
[00253] Cette correction est effectuée dans une base de correction c adaptée aux aberrations à corriger. La base de correction c est une base de correction en entrée ou une base de correction en sortie. [00254] Des exemples de base de correction sont : [00255] - une base des ondes planes ou base de Fourier spatiale, [00256] - une base des transducteurs u, [00257] - une base correspondant au lieu supposé des aberrateurs dans le milieu, c’est-à-dire par exemple un plan entre le plan des transducteurs (base des transducteurs u) et le plan de focalisation (base de focalisation x), [00258] - une base correspondant à un plan déterminé par optimisation, par exemple par une matrice de corrélation dont la première valeur propre est maximale. [00259] [00260] Matrice de réflexion duale Rcx(z,f, #^) [00261] Selon un mode de réalisation du procédé de la présente divulgation, la projection aller permet de déterminer une matrice de réflexion duale Rc(z,f, #^). Cette projection aller peut être effectuée par : [00262] un produit matriciel entre une matrice de passage P et la matrice de réflexion focalisée Rxx(z, f, #^), c’est-à-dire : [00263] ^^^(^, ^, #^) = ^(^, ^) × ^^^(^, ^, #^) [00264] où : ^(^, ^) = [^(c, x, ^, ^)] est la matrice de passage à chaque fréquence f entre la base focalisée (x) à la profondeur z et la base de correction (c). [00265] [00266] La matrice de passage P dépend de la base de correction c utilisée. [00267] [00268] Dans le cas d’une base de correction correspondant à une base des ondes planes (c=k), la matrice de passage P est l’opérateur de transformée de Fourier. [00269] Dans le cas d’un réseau 10 de transducteurs 11 de type linéaire pour générer une image à deux dimensions, les coefficients de cette matrice de passage P peuvent s’écrire selon :
[00270] ^(^^ , ^, ^, ^) = ^(^^ , ^) = ^^^(−^^^^) [00271] où ^^, la composante transverse du vecteur d’onde k associé à chaque onde plane. [00272] Dans le cas d’un réseau 10 de transducteurs 11 de type matriciel pour générer une image en trois dimensions, les coefficients de cette matrice de passage P peuvent s’écrire :
[00275] ^|| est la composante transverse du vecteur d’onde k associé à chaque onde plane, et [00276] ^ = (^, ^), le vecteur position transverse. [00277] [00278] Dans le cas d’une base de correction correspondant à une base des transducteurs (c=u), les coefficients de la matrice de passage P correspondent à la dérivée normale de la fonction de Green reliant chaque point focal de position spatiale (x, z) et chaque transducteur de position spatiale (u, 0). [00279] Dans le cas d’un réseau de transducteurs de type linéaire pour générer une image à deux dimensions, les coefficients de la matrice de passage P peuvent s’écrire selon : [00280] ^(^, ^, ^, ^) = ∇^^^^(^, ^, ^) [00281] où ∇^ est le gradient projeté suivant la direction de profondeur z, et [00282] ^^^(^, ^) est la fonction de Green 2D qui relie chaque transducteur ^ = (^, 0) à chaque point ^ du milieu M , avec :
[00284] où ^^ = 2^^/^^ est le nombre d’onde, [00285]
est la fonction de Hankel du 1 er ordre dont l’expression asymptotique est la suivante : ℋ^(2^^|^ − ^|/^^) =
[00286] Dans le cas d’un réseau de transducteurs de type matriciel pour générer une image en trois dimensions, les coefficients de cette matrice de passage P peuvent s’écrire : [00287] ^(^, ^, ^, ^) = ∇^^^^(^, ^, ^)
[00288] où ^^^(^, ^) est la fonction de Green 2D qui relie chaque transducteur ^ = (^^ , ^^ , 0) à chaque point ^ = (^, ^) du milieu M, avec :
[00290] Les coefficients de la matrice de passage P s’écrivent donc dans ce cas de la manière suivante : [00291]
[00292] [00293] Matrice de corrélation [00294] Selon un premier mode de réalisation, cette étape de calcul comprend : [00295] - la construction d’une matrice de corrélation C(x,z) à partir des matrice de réflexion duale Rcx(z,f, , #^) pour chaque point (x,z) du champ de vision [00296] - l’analyse de cette matrice de corrélation C(x,z) pour déterminer la loi de focalisation spatio-fréquentielle ^(x,z) pour chaque point du champ de vision. [00297] [00298] Selon une première variante, la matrice de corrélation C(x,z) est déterminée dans la base de correction c et dans le domaine fréquentiel, par le calcul suivant des éléments de la matrice de corrélation C = Ccc : [00299] ^({^, ^}, {^’, ^’}, ^, ^) = ∑#^ ^(^, ^, ^, ^, #^) ^∗(^, ^^, ^, ^^ , #^) [00300] Où * est l’opérateur de conjugaison. [00301] Cette opération permet de corréler les champs réfléchis dans la base de correction pour chaque source virtuelle en (x,z). Cette corrélation est moyennée sur les différentes frames #^ (i.e différentes réalisations du speckle) afin de s’affranchir de la réflectivité aléatoire du milieu et ainsi synthétiser une étoile guide cohérente à partir des différentes réalisations #^ du speckle. [00302] [00303] Selon une deuxième variante, la matrice de corrélation C est déterminée dans la base des frames #^, par le calcul suivant des éléments de la matrice de correction C = ^##:
[00305] * est l’opérateur de conjugaison [00306] [00307] Estimation de la loi de focalisation [00308] Selon une première variante, l’analyse de la matrice de corrélation C(x,z) est effectuée par une décomposition en valeurs propres de la matrice de corrélation C(x,z), et la loi de correction spatio-fréquentielle ^(x,z) est le premier vecteur propre ^^ de la matrice de corrélation C(x,z) dans la base de correction (c), c’est-à-dire C = Ccc. [00309] La matrice de corrélation étant hermitienne (^ = ^^), ses valeurs propres sont réelles et positives. [00310] La matrice de corrélation Ccc(x,z) peut ainsi s’écrire :
[00312] ou en termes de coefficients matriciels : [00313] ^({^, ^}, {^’, ^’}, ^, ^) = ∑^ ^^^^(^, ^)^^∗(^′, ^′) [00314] avec ^^ correspondant aux vecteurs propres de la matrice de corrélation C [00315] ^^ correspondant aux valeurs propres réelles et positives de la matrice de corrélation Ccc(x,z) rangées dans un ordre décroissant :
> ^^ > ⋯ > ^^ [00316] On a alors la loi de correction spatio-fréquentielle ^(x, z) qui est égale au premier vecteur propre, i.e. ^(x, z) = ^^ ; ou à sa version normalisée, ^(x, z) = exp(^arg{^^}), i.e. une loi de correction spatio-fréquentielle dont les coefficients sont d’amplitude unité mais dont la phase est égale à celle de ^^ (le symbole arg{^} désigne la phase du vecteur X ; ou à une correction de type filtre inverse, ^(x, z) = exp(^arg{^^})/|^^|. La première option est à privilégier si on est en présence d’un mauvais rapport signal à bruit (filtre adapté). En général, la seconde option sera toutefois privilégiée afin que la correction n’agisse pas comme un filtre en amplitude mais permette la correction des seules distorsions de phase. Enfin la troisième option est pertinente quand le milieu aberrateur atténue de manière inhomogène certaines composantes et/ou fréquences du
champ que l’on souhaite réhausser afin d’avoir un estimateur plus fidèle de la réflectivité in fine. [00317] [00318] Selon une deuxième variante, l’analyse de la matrice de corrélation Ccc(x,z) est effectuée par une décomposition en valeurs singulières de la matrice de réflexion duale réarrangée de la manière suivante : [00319] ^^#(x, z) = [^({^, ^}, #^ , x, z)] [00320] La décomposition en valeurs propres de la matrice de corrélation Ccc(x,z) réalisée dans la première variante est en effet équivalente à la décomposition en valeurs singulières (SVD) de chaque matrice de réflexion duale ^^# [00321] La décomposition en valeurs singulières s’applique sur des matrices de forme rectangulaire, et appliquée à la matrice de réflexion duale ^^# en chaque point (x,z), elle s’écrit ainsi :
[00323] ou en termes de coefficients matriciels :
[00325] avec ^^ = ^^^(^, ^)^ correspondant aux vecteurs singuliers de la matrice de fréflexion duale ^^#(^, ^) dans la base de correction, ou de manière équivalente, aux vecteurs propres de la matrice Ccc tels que définis dans la première variante. [00326] ^^ = ^^^(#^)^ correspondant aux vecteurs singuliers de la matrice de réflexion duale ^^# dans la base des frames, [00327] [00328] ^^ correspondant aux valeurs singulières de la matrice de réflexion duale ^^# qui sont, par définition, égales à la racine carrée des valeurs propres ^^ de la matrice de corrélation C telles que définies dans la première variante: ^^ = ^ ^ ^. [00329] [00330] On a alors la loi de correction spatio-fréquentielle ^(x, z) qui est égale au premier vecteur singulier de la matrice de distorsion duale ^^#, i.e. ^(^, ^) =
^^ ; ou à sa version normalisée, ^(^, ^) = exp(^arg{^^}), i.e. une loi de correction spatio-fréquentielle dont les coefficients sont d’amplitude unité mais dont la phase est égale à celle de ^^ (le symbole arg{^} désigne la phase du vecteur X) ; ou à une correction de type filtre inverse, ^(^, ^) = exp(^arg{^^})/|^^| . [00331] L’intérêt de la décomposition en valeurs singulières de la matrice de réflexion duale^^#, par rapport à une décomposition en valeurs propres de la matrice de corrélation Ccc est la rapidité du calcul des algorithmes numériques de la décomposition en valeurs singulières [00332] [00333] Cette recherche de la loi de correction spatio-fréquentielle ^(x, z) est aussi équivalente à résoudre l’équation suivante : [00334] a^(x, z, ) = ^^^(x, z) × ^(x, z) [00335] où × est le produit matriciel et a est une constante [00336] de manière itérative par l’expression suivante, qui correspond à un calcul par retournement temporel itératif : [00337] ^^^^(x, z) = ^^^(x, z) × ^^(x, z), [00338] avec ^^ un front d’onde arbitraire, par exemple ^^ = [1 ⋯ 1]^ [00339] [00340] Alors, la loi de correction spatio-fréquentielle ^(x, z) est obtenue par : [00341] ^(x, z) = ^ ^^ →^ ^ ^^ (x, z), [00342] ou sa version normalisée :
[00344] ou sa version filtre inverse :
[00346] Pour ^ → ∞, l’algorithme de retournement temporel itératif converge vers le même premier vecteur propre ^^ de la matrice ^^^. En pratique, il peut y avoir un intérêt à passer par un algorithme de retournement temporel itératif plutôt que par une SVD car il peut converger au bout de quelques itérations, d’où une plus grande rapidité de calcul. [00347]
[00348] Selon une troisième variante, l’analyse de la matrice de corrélation Ccc est effectuée par la résolution de l’équation suivante : [00349] ^(x, z) = exp(^ arg{^^^(x, z) × ^(x, z)}) [00350] de manière itérative par l’expression suivante, qui correspond à un calcul par retournement de phase itératif : [00351] ^^^^(x, z) = exp(^ arg{^^^(x, z) × ^^(x, z)}) [00352] où × est le produit matriciel, [00353] avec : ^^ un front d’onde arbitraire, par exemple ^^ = [1 ⋯ 1] ^. [00354] Alors, la loi de correction spatio-fréquentielle ^(x, z) est obtenue par :
[00358] L’intérêt d’un algorithme de retournement de phase itératif par rapport aux alternatives précédentes est d’être un estimateur plus fiable de la phase de la loi de correction ^(x, z) et donc d’accéder in fine à une meilleure compensation des distorsions de phase induites par l’aberrateur. [00359] [00360] Selon une quatrième variante, l’analyse de la matrice de corrélation C## est effectuée par la résolution de l’équation suivante : [00361] ^(x, z) = exp(^ arg{^##(x, z) × ^(x, z)}) [00362] où × est le produit matriciel, [00363] de manière itérative par l’expression suivante : [00364] ^^^^(x, z) = exp(^ arg{^##(x, z) × ^^(x, z)}) [00365] où × est le produit matriciel, [00366] avec ^^ un front d’onde arbitraire, par exemple ^^ = [1 ⋯ 1] ^ [00367] ce qui permet d’obtenir le vecteur W(x, z) suivant : [00368] ^(x, z) = ^ ^^ →^ ^ ^^(x, z) [00369] Ce vecteur ^(x, z) = [^(#^ , ^, ^)] défini dans la base des frames contient la phase de chaque étoile guide incohérente synthétisée par focalisation au point (x,z) pour chaque frame #^ .
[00370] Le conjugué en phase de ce vecteur ^(x, z) peut alors être exploité pour rephaser chaque étoile virtuelle incohérente de sorte à pouvoir les recombiner de manière cohérente et obtenir ainsi un estimateur de la loi de correction spatio-fréquentielle ^(x, z)non biaisé par la réflectivité aléatoire du milieu. Mathématiquement, cette opération s’écrit de la manière suivante :
[00372] L’intérêt de cette approche par rapport à une SVD de la matrice distorsion Rc# (deuxième variante) ou de l’algorithme de retournement phase itératif (troisième variante) est de converger vers une loi de correction non biaisée par l’amplitude plus importante du signal ultrasonore sur certaine frame de la matrice de réflexion (amplitude importante provoquée par le passage d’un diffuseur brillant tel qu’une bulle ou un problème expérimental). [00373] La figure 7 illustre un exemple d’extraction des lois d'aberrations par retournement de phase itératif dans le speckle dynamique. [A] La matrice de réflexion focalisée aberrée ^^^(^, #) [B] est projetée en sortie dans une base de correction (c), ici illustrée par la base des transducteurs (u). ^ = (^, ^) est le vecteur position transverse. Chaque frame de la matrice duale ^^^(^, #) constitue une réalisation dynamique du désordre. [C] Cette matrice est réarrangée suivant les dimensions (c) et (#) pour construire une matrice ^^#(^^^ , ^) en chaque point (^^^, ^) et calculer [D] la matrice de corrélation associée ^^^(^^^, ^). Une analyse par RPI (retournement de phase itératif) de la matrice ^^^(^^^, ^) permet d'extraire de ces fronts d'ondes une estimation ^(^^^, ^) des lois d'aberration en chaque point du champ de vision. Le conjugué en phase des lois d'aberrations ^∗(^^^, ^) est utilisé comme loi de focalisation adaptative en réception pour chaque frame permettant de retrouver une tache focale de sortie corrigée. Bien qu'illustrée en base (c) = (u) ici, le raisonnement peut être étendu à une autre base de correction, en particulier celle des ondes planes (k). Cette approche peut être décrite en construisant une matrice distorsion ou schématisée comme sur la figure 7 à l’aide de la matrice de réflexion duale. Les deux approches sont équivalentes : dans le premier cas on
extrait une loi d’aberration, tandis que dans le second, on extrait une loi de focalisation (aberration et courbure géométrique). [00374] La figure 8 illustres des exemples de lois d'aberrations transverses locales extraites dans le speckle dynamique. [A] Plusieurs régions du champ de vision choisies à titre illustratif et leurs [B] lois d'aberrations transverses correspondantes. L’association se fait deux à deux par code couleur. Le contour de chaque zone en transparence sur l’image [A] définit la zone sur laquelle la loi est estimée tandis que celle opaque au centre sur l’image [A] délimite la zone sur laquelle cette loi sera appliquée pendant la correction. Ces zones ont également une extension selon l'axe y qui n'est pas représentée ici. Le facteur de Strehl « S » associé à chacune des lois d'aberration sur les images [B] est indiqué. [00375] [00376] Correction des aberrations dans les matrices de réflexion [00377] A partir de la loi de focalisation obtenue, on va compenser en post- traitement les aberrations en recalculant des matrices de réflexion duales et focalisées à partir de la loi de focalisation estimée. [00378] Matrice de réflexion duale corrigée [00379] La matrice de réflexion duale corrigée ^^ ^^ (^, ^, #^) = [^′(^, ^, ^, ^, #^)] est déterminée en effectuant le produit terme à terme entre la matrice de réflexion duale Rc(z, f) et le conjugué en phase de la loi de focalisation ^, c’est- à-dire par : [00380] ^^ ^^ = ^^^ ∘ ^∗ [00381] où l’exposant * représente l’opération de conjugaison de phase [00382] le symbole ∘ est le produit d’Hadamard, c’est-à-dire le produit matriciel terme à terme de coefficients des matrices, tel que [00383] ^′(^, ^, ^, ^, , #^) = ^(^, ^, ^, ^, #^)^∗(^, ^, ^, ^). [00384] [00385] Matrice de réflexion focalisée corrigée [00386]
[00387] Une matrice de réflexion focalisée corrigée ^ ^ ^^ (^, ^) est alors déterminée par projection retour de la matrice de réflexion duale corrigée ^^ ^^ (z,f) vers la base focalisée (x). [00388] La projection retour est effectuée par un produit matriciel entre la matrice de passage P définie plus haut et la matrice de réflexion focalisée ^ ^ ^^ (z, f), c’est-à-dire : [00389] ^ ^ ^^ (^, ^, #^) = ^^(^, ^) × ^^ ^^ (^, ^, #^) [00390] Où l’exposant † désigne l’opération matricielle de trans-conjugaison. [00391] [00392] Itérations du traitement de correction [00393] Selon le mode de réalisation du procédé, les étapes du traitement de correction, c’est-à-dire les étapes de : [00394] - de détermination de la loi de correction ^(x, z) ladite étape comprenant éventuellement les étapes de détermination de la matrice de réflexion duale Rcx(z, f, #^), S160 de calcul de la loi de correction ^(x, z)et de détermination de la matrice de réflexion duale corrigée ^^ ^^ (^, ^, #^), et [00395] - de détermination de la matrice de réflexion focalisée corrigée
[00396] sont itérées plusieurs fois (deux fois ou plus de deux fois). [00397] A chaque itération, la projection aller utilise la matrice de réflexion focalisée corrigée ^ ^ ^^ (^, #^) obtenue durant la projection retour de l’itération précédente à la place de la matrice de réflexion focalisée Rxx(z, #^). [00398] Ainsi, à chaque itération, la loi de correction de spatio-fréquentielle est améliorée pour prendre en compte de mieux en mieux une ou plusieurs aberrations du milieu. [00399] [00400] Selon une première variante de ce processus itératif, à chaque itération de l’étape de détermination de la matrice de réflexion duale Rcx(z, on utilise une base de correction c différente, par exemple pour corriger des aberrations différentes localisées en des lieux différents du milieu. [00401] Par exemple, le milieu peut être discrétisé ou modélisé par une succession de couches selon la direction de la profondeur z, et les bases de
correction c des itérations correspondent à des plans de ces couches successives. Autrement dit, on applique au cours des itérations des corrections correspondant à une pluralité d’aberrations du milieu. [00402] [00403] Selon une deuxième variante de ce processus itératif, à chaque itération de l’étape de détermination de la matrice de réflexion duale, on utilise une projection aller soit vers une base de correction en entrée, soit vers une base de correction en sortie de la matrice de réflexion. Dans ce dernier cas, la projection de la matrice Rxx(z,f, #^) vers la base de correction est effectuée de la manière suivante : [00404] ^^^(^, ^, #^) = ^(^, ^) × ^⊺ ^^ (^, ^, #^) [00405] où le symbole ⊺ désigne l’opération matricielle de transposition. [00406] Dans la succession des itérations, on peut alterner entre l’utilisation d’une base de correction en entrée et une base de correction en sortie. Ainsi, la loi de correction de spatio-temporelle ^ est améliorée à chaque itération, et la correction des aberrations est améliorée. [00407] [00408] Variantes du processus de correction des aberrations [00409] Suivant les conditions expérimentales, le processus d’aberrations pourra être simplifié ou complexifié suivant plusieurs variantes décrites ci- dessous. [00410] [00411] Compensation spatiale des aberrations [00412] Si le milieu en amont du plan focal n’induit pas de dispersion temporelle des échos (absence de réverbérations, diffusion multiple négligeable), il n’y a pas forcément d’intérêt à considérer de manière indépendante les composantes fréquentielles des signaux confocaux. Dans ce cas, on peut considérer la matrice de réflexion focalisée large bande,
[00414] comme point de départ du processus de correction des aberrations. La projection de la matrice de réflexion dans la base de correction se fait alors en considérant le propagateur à la fréquence centrale. Tout le reste du
processus est identique à ce qui a été décrit ci-dessus. Seule la dépendance fréquentielle des différentes quantités exprimées n’a plus lieu d’être. [00415] [00416] Compensation temporelle des aberrations [00417] Si, au contraire, le milieu induit une dispersion principalement temporelle des échos ultrasonores et peu d’aberrations spatiales, il n’y a pas d’intérêt à considérer la matrice de réflexion focalisée mais seulement son signal confocal : [00418] ^(^, ^, ^, #^) = ^(^, ^, ^, ^, #^) [00419] La matrice de corrélation à considérer dans ce cas pour l’obtention de la loi de correction, ^(x, z) = [Φ(f, x, z)], est donnée par :
[00421] [00422] Exploitation de l’effet mémoire spatial [00423] Le nombre de frames peut être insuffisant pour une estimation correcte de la loi de focalisation. Dans ce cas on pourra considérer, à la place de la matrice de réflexion duale, la matrice distorsion associée, ^^^(^, ^, #^) = [ ^^^(^, ^, ^, ^, #^)], et moyenner la matrice de corrélation non plus seulement sur les différentes frames mais également sur des pixels adjacents apparentant à
[00426] et ^^^^ = ^^^^^(^, ^, ^, ^)^, une matrice de référence d’un milieu modèle dans lequel la vitesse du son est c0 vitesse du son attendue du milieu et dans lequel un réflecteur plan est positionné à la profondeur z. Une moyenne de la matrice de corrélation sur plusieurs grains de speckle peut accélérer la convergence de la matrice de corrélation vers la matrice de covariance associée. Il en résulte une estimation plus précise vers la loi de focalisation. [00427] [00428] Construction d’images du milieu [00429] Image confocale dynamique
[00430] Selon un mode de réalisation du procédé de construction ultrasonore de la présente divulgation, le procédé comprend en outre une étape de : [00431] - détermination du signal confocal dynamique Sc(x,z, #^) d’un point de position spatiale (x, z), à partir des coefficients diagonaux de la matrice de réflexion focalisée corrigée ^ ^ ^^ (z,#^) intégrés sur la bande passante des signaux ultrasonores, c’est-à-dire par la combinaison des réponses confocales du point à plusieurs fréquences f. On a ainsi par exemple le calcul suivant : [00432] ^^(^, ^, #^) = ∑^ ^^^ ^ (^, ^, ^, ^, #^) [00433] L’intensité de cette quantité, [00434] ^^(^, ^, #^) = |^^(^, ^, #^)| ^, [00435] déterminée en une pluralité de points (x,z) permet de construire une image confocale du milieu corrigée, qui correspond à une image échographique classique débarrassée des problèmes d’aberrations, de réverbérations et de dispersion fréquentielle de la vitesse du son dans le milieu étudié. [00436] Sur la figure 9 est illustrée une comparaison des images confocales dynamiques avant et après correction. Plusieurs exemples d'images confocales originales (à gauche) ainsi que leurs équivalents après correction (à droite) sont montrés pour différentes frames, avec microbulles [A-C] ou sans [D]. Les bienfaits de la méthode sont appréciables notamment sur l’image des microbulles (points brillants) qui apparaissent mieux résolus et plus contrastés après correction. [00437] [00438] Image Doppler de puissance (« Power Doppler » en anglais) [00439] Selon un mode de réalisation du procédé de construction ultrasonore de la présente divulgation, le procédé comprend en outre une étape de : [00440] - détermination d’une image Power Doppler IP(x,z) d’un point de position spatiale (x, z), à partir de la somme sur les différentes frames de l’intensité confocale dynamique ^^(^, ^, #^) mesurée en chaque point (x,z). On a ainsi par exemple le calcul suivant :
[00442] L’intensité précédente déterminée en une pluralité de points permet de construire une image power Doppler corrigée du milieu, qui quantifie le mouvement des tissus au cours de la séquence d’acquisitions. [00443] La figure 10 illustre deux images de type Doppler de puissance. On distingue le gain en contraste apporté par la correction. [A] est un affichage PWD des 400 frames du bloc d’application avant et [B] après correction. [00444] [00445] Image Doppler Directionnelle [00446] Selon un mode de réalisation du procédé de construction ultrasonore de la présente divulgation, le procédé comprend en outre une étape de : [00447] - détermination d’une image Doppler directionnelle ID(x,z) d’un point d’image échographique de position spatiale (x, z), à partir de la transformée de Fourier du signal confocal complexe ^^(^, ^, #) suivant la dimension temporelle #^ des frames:
[00449] Où la fréquence ^ est la variable conjuguée au temps d’acquisition #. [00450] L’image power doppler initiale peut ainsi être séparée en une composante positive et une composante négative du déplacement:
[00453] [00454] Ces deux images peuvent ensuite être combinées pour former une image Doppler colorée ID(x,z) dont le bleu encode la composante négative et le route la composante positive du mouvement des tissus en chaque point (x,z). [00455] [00456] Carte de la dynamique locale [00457] Selon un mode de réalisation du procédé de construction ultrasonore de la présente divulgation, le procédé comprend en outre une étape de : [00458] - détermination d’une carte de dynamique locale d’un point d’image échographique de position spatiale (x, z), à partir de la fréquence moyenne ^D(x,z) de la transformée de Fourier ^^(^, ^, ^) du signal confocal complexe:
Cette fréquence moyenne est reliée à la vitesse moyenne des diffuseurs en chaque point du milieu. Sur la figure 11 sont représentées ces cartes de la dynamique locale permettant une comparaison avant [A] et après [B] correction des aberrations transverses. La correction des aberrations permet une estimation plus fine de la vitesse de déplacement des diffuseurs dans le milieu. [00460] [00461] Image confocale statique [00462] Au-delà de la composante dynamique des données acquises, les lois de focalisation obtenues peuvent être exploitées pour également corriger la matrice de réflexion brute, i.e la matrice de réflexion considérée avant le filtrage de la composante statique des données. [00463] [00464] Le domaine d'application préférentiel est l'imagerie ultrasonore dans son ensemble et notamment ses modes Doppler et ULM (« Ultrasound Localization Microscopy »). L'invention peut aussi bien servir à imager le cerveau humain à travers un crâne que le réseau vasculaire du foie pour l’imagerie médicale. Elle peut servir à imager des processus dynamiques en contrôle non destructif. L’invention a également des applications potentielles dans différents domaines de la physique des ondes, comme la tomographie à cohérence optique pour la lumière ou la sismologie pour les ondes sismiques. [00465] [00466] Bien entendu, l’invention n’est pas limitée aux exemples qui viennent d’être décrits. De nombreuses modifications peuvent être apportées à ces exemples sans sortir du cadre de la présente invention telle que décrite.