FR2936889A1 - Procede pour le recalage d'un ensemble de points dans des images - Google Patents

Abstract

L'invention concerne un procédé pour le recalage d'un ensemble de points sources (S) dans une image de référence représentant un tissu organique et d'un ensemble de points destinations (D) dans une image courante représentant le tissu organique déformé, le procédé comprenant la détermination d'une fonction de recalage (R) par minimisation itérative d'une énergie de recalage Ereg, la fonction de recalage correspondant à une transformation permettant de passer des points sources (S) de l'image de référence aux points destinations (D) correspondants dans l'image courante, caractérisé en ce que la fonction de recalage (R) est une fonction de recalage élastique satisfaisant des conditions pour garantir le respect de critères physiques du tissu organique, les conditions comprenant : - le Jacobien de la fonction de recalage élastique est supérieur à 0, - la fonction de recalage est continûment différentiable, - la fonction de recalage est bijective.

Description

DOMAINE DE L'INVENTION

La présente invention concerne le domaine de la chirurgie assistée par ordinateur, et plus particulièrement un système et une méthode permettant le recalage d'un tissu 5 organique dans deux images.

PRESENTATION DE L'ART ANTERIEUR

La localisation précise d'une cible est essentielle au cours d'une intervention 10 chirurgicale pour réduire le taux de morbidité. Au cours d'une intervention chirurgicale, une déformation des tissus organiques peut se produire. Par exemple, au cours d'une intervention chirurgicale d'ablation d'une tumeur cérébrale, cette déformation peut intervenir après qu'une ouverture ait été réalisée dans le 15 crâne du patient. Cette déformation peut être liée à des phénomènes physiques (gravité, perte du liquide céphalorachidien, action du neurochirurgien, etc.) ou à des phénomènes physiologiques (tuméfaction due aux médicaments osmotiques, anesthésiques, etc.), dont certains ne sont pas connus à ce jour. Afin de compenser cette déformation des tissus organiques et de retrouver les 20 données préopératoires du patient pendant la procédure chirurgicale, une compensation de la déformation des tissus organiques peut être réalisée sur le principe suivant. Une première image d'une région des tissus organiques est acquise avant l'opération chirurgicale. Pendant l'opération chirurgicale, une seconde image de la même région est acquise. Des moyens de traitement du système traitent ces deux images pour 25 aligner les points dans les deux images afin de compenser la déformation des tissus organiques. Afin d'aligner les points dans les deux images, les moyens de traitement mettent en oeuvre un procédé de recalage d'images. Toutefois, les procédés de recalage ne permettent pas d'aligner avec précision les points dans deux images lorsque les caractéristiques des données présentes dans les 30 deux images sont hétérogènes. Un but de la présente invention est de proposer un procédé de recalage d'un ensemble de points sources dans une première image et d'un ensemble de points destinations dans une deuxième image, l'ensemble de points sources étant de nature distincte de l'ensemble de points destinations, du fait de la différence entre les zones de 35 recouvrement des données ou de la présence de bruit. 1 RESUME DE L'INVENTION

A cet effet l'invention prévoit un procédé pour le recalage d'un ensemble de points sources S dans une image de référence représentant un tissu organique et d'un ensemble de points destinations D dans une image courante représentant le tissu organique déformé, le procédé comprenant la détermination d'une fonction de recalage R par minimisation itérative d'une énergie de recalage Ereg, la fonction de recalage correspondant à une transformation permettant de passer des points sources S de l'image de référence aux points destinations D correspondants dans l'image courante, caractérisé en ce que la fonction de recalage R est une fonction de recalage élastique satisfaisant des conditions pour garantir le respect de critères physiques du tissu organique, les conditions comprenant : - le Jacobien de la fonction de recalage élastique est supérieur à 0, - la fonction de recalage est continûment différentiable, - la fonction de recalage est bijective. Ainsi, la fonction de recalage élastique est déterminée de sorte à garantir le respect de critères physiques associés au tissu organique que l'on souhaite recaler. Ces critères physiques à respecter sont les suivants : - non repliement du tissu organique sur lui-même (ce qui se traduit par le fait que le jacobien de la fonction de recalage doit être supérieur à 0), - régularité de la propagation de la déformation du tissu organique (ce qui se traduit par le fait que la fonction de recalage est continûment différentiable), - non-superposition de deux points distincts de l'espace initial en un point unique de l'espace déformé et portée de la fonction de déformation étendue à l'espace entier (ce qui se traduit par le fait que la fonction est bijective). Des aspects préférés mais non limitatifs de la présente invention sont les suivants : 30 - l'ensemble de points sources est de nature distincte de l'ensemble de points destinations, - l'énergie de recalage Ereg de la fonction de recalage (R) satisfait la relation suivante : Ereg = E d(R(s), D) ses 35 avec : - s un point de l'ensemble de points sources S, 2 20 25 - R(s) le transformé du point s par la fonction de recalage R, - d(R(s), D) la distance euclidienne entre le transformé du point s par la fonction de recalage et l'ensemble de points destinations, noté D, - le procédé comprend une étape de filtrage direct dans laquelle chaque point s de l'ensemble de points sources S dont l'association avec un point d de l'ensemble de points destinations D vérifie un critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points source S', - le critère de rejet est que la distance entre : 0 le transformé du point s de l'ensemble de points sources S par la fonction de recalage et o le point correspondant d de l'ensemble de points destinations D est supérieure à une valeur maximale prédéterminée, - Le procédé comprend en outre une étape de filtrage inverse dans laquelle chaque point d de l'ensemble de points destinations D dont l'association avec un point s de l'ensemble de points sources S vérifie un autre critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points destinations D', - l'autre critère de rejet est que la distance entre : o le transformé inverse du point d de l'ensemble de points destinations D par la fonction de recalage inverse et o le point correspondant s' de l'ensemble filtré de points sources S' est supérieure à la valeur maximale prédéterminée. L'invention concerne également un système pour le recalage d'un ensemble de points source dans une image de référence représentant un tissu organique et d'un ensemble de points destinations dans une image courante représentant le tissu organique déformé, caractérisé en ce qu'il comprend des moyens pour la mise en oeuvre du procédé décrit ci-dessus. L'invention concerne également un produit programme d'ordinateur comprenant des instructions de code programme enregistré sur un support utilisable dans un ordinateur, caractérisé en ce qu'il comprend des instructions pour la mise en oeuvre du procédé décrit ci-dessus.

BREVE DESCRIPTION DES DESSINS D'autres avantages et caractéristiques ressortiront mieux de la description qui va suivre de plusieurs variantes d'exécution, données à titre d'exemples non limitatifs, à partir des dessins annexés sur lesquels - la figure 1 illustre un mode de réalisation du procédé selon l'invention, - la figure 2 illustre un mode de réalisation du système pour la mise en oeuvre du procédé selon l'invention, - la figure 3 illustre une déformation d'un carré unitaire avec un repliement d'espace, - la figure 4 illustre en deux dimensions une approche itérative multi-échelle décrite en relation avec l'étape de recalage de tissu organique, - la figure 5 illustre des cellules insérées en marge d'une discrétisation pour deux niveaux de raffinement successifs, - la figure 6 illustre une fonction de forme w, - la figure 7 illustre un exemple d'approximation d'une courbe d'énergie de recalage par une parabole, - la figure 8 illustre trois étapes de l'approximation de l'inverse de la fonction de recalage, - la figure 9 illustre des étapes de filtrage direct et indirect.

DESCRIPTION DETAILLEE DE L'INVENTION

On va maintenant décrire plus en détail un mode de réalisation du procédé de recalage selon l'invention. Ce procédé de recalage sera présenté dans le cadre d'un procédé de traitement d'image pour estimer une déformation d'un cerveau.

Toutefois, il est bien entendu que ce procédé de recalage est totalement indépendant du procédé de traitement pour estimer une déformation du cerveau. Plus précisément, le procédé de recalage selon l'invention : - peut être mis en oeuvre pour le recalage d'autres tissus organiques que le cerveau, et/ou - peut être mis en oeuvre dans d'autres types de procédés de traitement d'image.

Principe général du procédé de traitement pour estimer une déformation du cerveau Comme expliqué précédemment, la localisation précise de la cible chirurgicale est essentielle pour réduire la morbidité au cours de l'exérèse chirurgicale d'une tumeur cérébrale. Lorsque les dimensions de la craniotomie sont importantes, une déformation des tissus mous du cerveau peut survenir au cours de l'intervention.

Du fait de cette déformation du cerveau, les données préopératoires ne correspondent plus à la réalité, et la neuronavigation s'en trouve fortement compromise. La présente invention permet de prendre en compte cette déformation et de calculer une position rectifiée des structures anatomiques du cerveau afin de localiser des ancillaires. Avant une intervention chirurgicale sur le cerveau d'un patient, une angiographie par résonance magnétique (ARM) du cerveau du patient peut être acquise. Au cours de l'intervention, suite à une déformation importante du cerveau, le chirurgien effectue un balayage échographique Doppler de la région d'intérêt.

Le procédé de traitement pour estimer une déformation du cerveau permet la mise à jour des données préopératoire en fonction de la déformation du cerveau au cours de l'intervention. Pour ce faire, on propose de suivre au cours de l'intervention la déformation de la structure anatomique particulière qu'est l'arbre vasculaire cérébral.

Du fait de sa répartition dans le volume du parenchyme, les déplacements de l'arbre vasculaire sont représentatifs des déformations globales des tissus. En effet, les artères et artérioles irriguant les tissus nerveux sont présentes à la fois dans les parties profondes de l'encéphale, mais aussi en surface, et plus particulièrement au voisinage de la tumeur dont le développement passe par une angiogenèse et donc la formation de nouveaux vaisseaux. De la localisation de l'arbre vasculaire par rapport à sa position initiale en début d'intervention il est possible de déduire les déplacements de ses vaisseaux, provoqués par la déformation du cerveau. Cette information, loin d'être dense dans le volume global, reste localisée à 25 certaines régions. On propose alors d'utiliser une modélisation biomécanique du cerveau du patient afin d'extrapoler une déformation globale à partir du déplacement de l'arbre vasculaire cérébral. Le champ de déplacements global ainsi obtenu est ensuite répercuté à l'ensemble 30 des données préopératoires afin de les mettre à jour.

Description d'un mode de réalisation du procédé pour estimer une déformation du cerveau En référence à la figure 1, on a illustré différentes étapes du procédé pour estimer 35 la déformation du cerveau.

Le procédé de traitement d'image pour estimer une déformation du cerveau comprend les étapes suivantes : - le traitement 100 d'une image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient, - le traitement 200 d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient, - la détermination 300, à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, d'un champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, - l'application 400 du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour estimer le déplacement global du cerveau du patient ; ce modèle sert d'interpolateur mécanique afin de calculer au cours de l'intervention l'ensemble des déplacements volumiques à l'intérieur du cerveau - la génération 500, à partir du déplacement du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé.

25 Détermination de l'arbre vasculaire cérébral de référence Comme décrit précédemment, le procédé comprend une étape 100 consistant à traiter une image du cerveau du patient acquise avant l'intervention chirurgicale. Ce traitement a pour but de localiser l'arbre vasculaire cérébral dans l'image acquise avant l'intervention chirurgicale. 30 Cet arbre vasculaire cérébral dans l'image acquise avant l'intervention chirurgicale correspond à l'arbre vasculaire avant la déformation du cerveau du patient, et est appelé dans la suite arbre vasculaire cérébral de référence . Préférentiellement, l'image acquise avant l'intervention chirurgicale est une image tridimensionnelle. 35 Préférentiellement encore, cette image tridimensionnelle est une angiographie par résonnance magnétique. 10 15 20 L'angiographie par résonance magnétique permet de révéler des contrastes importants entre les tissus. L'analyse de l'angiographie par résonance magnétique permet de distinguer avec précision la nature des tissus observés. En particulier, les données qu'elle contient permettent de localiser l'arbre vasculaire cérébral de référence. Par ailleurs, la technique permettant l'acquisition d'une angiographie par résonance magnétique présente l'avantage d'être non invasive et donc sans dangers pour le patient, contrairement par exemple à des techniques émettant un rayonnement ionisant comme la tomographie. Une fois le patient installé sur la table opératoire, le traitement de l'image tridimensionnelle pour obtenir l'arborescence artérielle cérébrale de référence comprend une étape d'orientation de l'image tridimensionnelle par rapport à la tête du patient. En d'autres termes, la position spatiale de l'arbre vasculaire de référence est déterminée par rapport à un référentiel de l'espace physique du patient. Plus précisément, l'arbre vasculaire cérébral de référence est transposé dans le référentiel physique du patient, à l'aide d'une technique dite de recalage rigide qui s'appuie sur la mise en correspondance de points de calibrages avec leurs équivalents segmentés dans l'angiographie par résonance magnétique acquise avant l'opération. On entend, dans le cadre de la présente invention, par recalage rigide , une transformation géométrique composée d'une rotation et d'une translation.

Un calcul permet alors de trouver la transformation entre les deux systèmes de coordonnées. Les points de calibrage utilisés peuvent être des points ou des surfaces anatomiques remarquables comme la pointe du nez, les arcades sourcilières, la surface des joues ou du front.

Les points de calibrage peuvent aussi être définis par une pluralité de pastilles adhésives (au moins trois et préférentiellement dix) collées sur la peau du patient avant l'acquisition de l'angiographie par résonance magnétique. Ainsi, ces pastilles sont acquises dans l'angiographie en même temps que le cerveau du patient. Les homologues radiométriques des points de calibrage sont identifiables dans l'angiographie par résonance magnétique. Une fois le patient installé sur la table d'opération, le chirurgien indique à un système de localisation la position de chaque point de calibrage au moyen d'un palpeur. Le système de localisation peut être un système optique de localisation (comprenant par exemple une caméra stéréoscopique), ou un système magnétique de localisation, ou tout autre système de localisation connu de l'homme du métier.

Connaissant la position de l'arbre vasculaire cérébral de référence par rapport aux points de calibrage, et connaissant la position des points de calibrage une fois le patient installé sur la table d'opération, le système est apte à calculer la position spatiale de l'arbre vasculaire de référence par rapport à la position spatiale du patient.

Le recalage entre les images et le patient, réalisé en début d'intervention, permet par la suite de localiser dans le champ opératoire tous les éléments identifiés dans l'image par résonance magnétique et, réciproquement, de transposer la position des ancillaires dans le volume de données afin de suivre et contrôler la bonne réalisation de l'intervention.

Détermination de l'arbre vasculaire cérébral courant Le procédé comprend également une étape 200 consistant à traiter une pluralité d'images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement une arborescence artérielle cérébrale du patient dite arborescence artérielle cérébrale courante . Durant l'opération, suite à une déformation importante du cerveau du patient, le chirurgien effectue une acquisition d'images bidimensionnelles. On entend, dans le cadre de la présente invention, par pluralité d'images bidimensionnelles , des images 2D dans un volume (3D) déterminé dont les plans d'acquisition sont sécants ou parallèles. Ces images 2D peuvent être obtenues à partir de tout dispositif médical d'acquisition d'images 2D ou 3D telles que des images échographiques 3D par exemple. Préférentiellement, ces images bidimensionnelles sont des images échographiques acquises en mode Doppler, à l'aide d'une sonde localisée par le système de localisation. Les images échographiques en mode Doppler permettent de visualiser les flux sanguins. Par ailleurs, la technique d'échographie Doppler présente l'avantage d'être bénigne pour le patient. L'échographie Doppler per-opératoire est donc un bon outil pour l'exploration des déplacements du cerveau au cours de l'intervention. Ce mode de visualisation repose sur l'effet Doppler et permet de localiser, en analysant les variations de fréquence des ultrasons émis, des déplacements au sein des tissus, tels que les écoulements de sang. La segmentation de ce type d'images est rendu plus aisée du fait du codage en couleur utilisé pour représenter les vélocités des écoulements. Cette modalité permet donc de visualiser dans le volume de l'encéphale les flux sanguins qui s'effectuent le long de l'arbre vasculaire et d'identifier ainsi la position des vaisseaux qui irriguent le parenchyme.

L'acquisition des images bidimensionnelles est réalisée par exemple comme suit. La sonde localisée est mise en contact avec le cerveau du patient à travers la craniotomie effectuée par le chirurgien. Un mouvement de rotation est manuellement imprimé à la sonde localisée de sorte que le plan échographique balaye la partie du cerveau la plus étendue possible autour de la tumeur. La sonde localisée comprend un marqueur permettant au système de localisation de définir sa position et son orientation par rapport au référentiel physique du patient, et donc de repositionner spatialement les images bidimensionnelles dans le référentiel physique du patient.

Après l'acquisition d'une série d'images bidimensionnelles, celles-ci sont traitées pour localiser l'arbre vasculaire cérébral courant du patient. Plus précisément, les centres des vaisseaux contenus dans chaque image bidimensionnelle sont sélectionnés de sorte à obtenir un nuage de points. Ce nuage de points représente au moins partiellement l'arbre vasculaire cérébral 15 courant du patient.

Recalaqe des arbres vasculaires de référence et courant Dans une autre étape 300 du procédé, l'arbre vasculaire de référence et courant sont mis en correspondance pour déterminer un champ de déformation représentant le 20 déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence. Plus spécifiquement, l'arbre vasculaire de référence est recalé par rapport à l'arbre vasculaire courant. On entend par recalage élastique , le fait de déterminer, à partir de jeux de données qui représentent deux images numériques d'une même zone, une transformation 25 permettant de passer de l'image de référence à l'image courante. En d'autres termes, il s'agit de superposer "au mieux" des structures contenues dans deux images comparables. Les données issues de l'image tridimensionnelle et des images bidimensionnelles sont de natures distinctes. 30 D'un côté, l'arbre vasculaire de référence obtenu à partir des données préopératoires est presque continu, l'ensemble de la tête du patient est contenu dans le volume de données, et la segmentation réalisée n'engendre pour ainsi dire pas d'artefacts dans la reconstruction de l'arbre vasculaire cérébral. D'un autre coté, l'arbre vasculaire courant obtenu à partir des images 35 bidimensionnelles per-opératoires est composé d'un nuage de points, ces images bidimensionnelles ne couvrent qu'une fraction du volume du cerveau. Selon l'habileté de l'opérateur, les images bidimensionnelles peuvent être plus ou moins régulièrement réparties dans l'espace, certaines régions comportant des acquisitions multiples tandis que d'autres ne comportent aucune acquisition. La tâche qui incombe à l'étape 300 de recalage des arbres vasculaires cérébraux 5 est donc d'estimer une transformation élastique entre : - la configuration courante, per-opératoire, de l'arbre vasculaire définie par un nuage de points épars, redondants, incomplet et bruités, et - la configuration de référence définie par les données complètes et 10 propres obtenues à partir de l'imagerie préopératoire. Le champ de déformation est alors calculé à partir de la correspondance établie entre l'arbre vasculaire de référence et l'arbre vasculaire courant. Chaque vecteur du champ de déformation a pour origine un point de l'arbre vasculaire de référence et désigne un point de l'arbre vasculaire courant, déformé. 15 Le recalage élastique est guidé par la minimisation d'une énergie de recalage. Cette énergie de recalage quantifie la similarité entre l'arbre vasculaire courant et l'arbre vasculaire de référence. Le recalage optimal est atteint pour une énergie de recalage nulle. L'objectif de l'étape de recalage élastique est de minimiser itérativement la valeur de cette énergie de 20 recalage en faisant évoluer les points de l'arbre vasculaire courant vers les points de l'arbre vasculaire de référence. La rapidité du calcul est obtenue grâce : - au calcul préalable d'une carte de distances - entre chaque point de l'arbre vasculaire courant et son homologue au sens radiométrique 25 (i.e. le point qui désigne la même structure) dans l'arbre vasculaire de référence - sur laquelle s'appuie l'évaluation de l'énergie de recalage et - la détermination de la direction de plus forte descente, calculée à chaque itération du processus de minimisation. 30 Etant donnée la nature éparse et irrégulière des données, la recherche d'une fonction de recalage élastique doit être contrainte afin de ne pas aboutir à un recalage aberrant. La fonction de recalage élastique est une transformation non-linéaire permettant de passer de l'arbre vasculaire courant à l'arbre vasculaire de référence. Plus précisément, des conditions liées à la nature des données que l'on cherche à 35 recaler sont déterminées pour garantir le respect de critères physiques essentiels.

S'agissant de données organiques, une première condition que doit satisfaire la fonction de recalage élastique concerne le non-repliement de l'espace. En effet, la déformation du cerveau du patient ne peut induire le repliement du cerveau du patient sur lui-même.

Une deuxième condition que doit satisfaire la fonction de recalage élastique concerne la continuité de la propagation de la déformation. En effet, les arbres vasculaires de référence et courant sont des observations du phénomène physique de déformation du cerveau. De façon similaire au processus de déformation du cerveau du patient, sans déchirement ou résection de tissu, la fonction de recalage doit être continue et même continument différentiable i.e. différentiable en tout point du domaine et ne présentant pas de discontinuité de différentielle. Une troisième condition que doit satisfaire la fonction de recalage est que deux points distincts de l'arbre vasculaire courant ne peuvent avoir la même image dans l'arbre vasculaire de référence, et vice versa.

La fonction de recalage élastique telle qu'elle est proposée dans l'invention présente l'avantage d'être inversible avec le degré de précision souhaité. Une fois cette fonction de recalage calculée, un post-traitement du résultat produit par le recalage permet d'éliminer les aberrations dues à la présence de bruit dans les données des images bi- ou tridimensionnelles per-opératoires.

L'ensemble de déplacements de l'arbre vasculaire finalement retenu forme un champ de déformation représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence. Le lecteur appréciera que l'étape de recalage peut être appliquée au recalage de structures autre qu'un arbre vasculaire cérébral. Notamment, le recalage élastique 25 présenté ci-dessus peut s'appliquer au recalage de tout type de tissu organique.

Extrapolation du champ de déplacements de l'arbre vasculaire cérébral à l'ensemble du volume de l'orqane par le biais d'un modèle biomécanique Dans une autre étape 500 du procédé, le champ de déformation de l'arbre 30 vasculaire déterminé est appliqué à un modèle biomécanique des tissus mous du cerveau du patient pour estimer le déplacement global du cerveau du patient. Ce modèle mathématique formalise une connaissance a priori du comportement de l'organe et permet d'en inférer le comportement global sous les contraintes de déformation locales déduites de l'observation per-opératoire, appelées conditions aux 35 limites.

Le modèle biomécanique (aussi connu sous le nom de modèle déformable) joue un double rôle. En premier lieu, le modèle biomécanique permet d'extrapoler les déplacements de l'arbre vasculaire à l'ensemble du volume du cerveau afin de simuler l'effet de la 5 déformation du cerveau. En deuxième lieu, le modèle biomécanique permet de réguler les données. En effet, l'estimation du champ de déplacements de l'arbre vasculaire est entachée d'erreurs dues au mode d'imagerie utilisé, aux limitations des algorithmes de segmentation ainsi qu'à celles du calcul du recalage élastique effectué pour mettre en correspondance l'arbre 10 vasculaire courant et l'arbre vasculaire de référence. Le modèle biomécanique permet d'harmoniser le champ de déplacements en supprimant les erreurs dans les données d'entrée du modèle biomécanique. Ceci permet de rendre le procédé résistant aux artefacts présents dans les données bruitées recueillies au cours de l'intervention (i.e. les images échographiques bi- ou tridimensionnelles). 15 L'originalité du procédé concerne ainsi la complémentarité entre un modèle biomécanique approximant les comportements des tissus cérébraux du patient et un suivi de structure anatomique potentiellement bruité mais représentatif du déplacement du cerveau grâce à l'arbre vasculaire cérébral. A l'issu de l'étape 500 du procédé, le déplacement global du cerveau du patient 20 est estimé. Une dernière étape 600 du procédé concerne la génération, à partir du déplacement global du cerveau estimé, d'au moins une image du cerveau du patient dans laquelle le déplacement du cerveau est compensé. Cette image est par exemple l'image tridimensionnelle (ou des images bidimensionnelles en coupe de cette image 25 tridimensionnelle) acquise avant l'opération à laquelle est appliquée la déformation globale estimée du cerveau du patient de sorte à mettre à jour les données préopératoire en fonction de la déformation globale estimée du cerveau du patient.

Vérification de la cohérence de la déformation calculée 30 Avantageusement, le procédé peut comprendre une étape de calcul des disparités entre des points de l'image du cerveau générée et une ou plusieurs images bidimensionnelles de contrôle acquises pendant l'opération. Ceci permet à l'utilisateur de vérifier la cohérence et la précision du calcul réalisé par le procédé décrit ci-dessus. Pour ce faire, l'utilisateur peut par exemple acquérir une 35 ou plusieurs images bidimensionnelles de contrôle du cerveau du patient pour calculer l'erreur de position entre les points des images bidimensionnelles de contrôle et leurs homologues dans l'image générée à partir de la déformation globale estimée du cerveau du patient. Par exemple, l'utilisateur peut scruter les tissus cérébraux dans la région d'intérêt tout en projetant dans les images bidimensionnelles de contrôle la position de l'arbre vasculaire déformé par le modèle. La superposition des deux données i.e. Doppler en temps réel obtenu par balayage de l'organe, d'une part et signal Doppler simulé, calculé à partir de l'incidence des coupes échographiques localisées par rapport à l'arbre vasculaire déformé par le système, d'autre part, permet à l'utilisateur d'évaluer la coïncidence entre le modèle et la réalité du patient. L'appréciation de la qualité de la déformation permet à l'utilisateur de prendre la décision de suivre ou non les indications du système sur la base des données préopératoires mises à jour. Le déplacement du cerveau étant un processus évolutif, cette validation est réalisée aussitôt après l'acquisition des données sur lesquelles se base l'algorithme de déformation afin d'en garantir la pertinence.

Le procédé décrit ci-dessus est parfaitement adapté à un suivi continu de la déformation des tissus mous du cerveau, tant par ses temps de réponse qui sont de l'ordre de quelques secondes, que par la répétabilité de la procédure d'acquisition de données qui ne nécessite pas une interruption majeure de l'intervention ou la mise en place d'équipements encombrants.

Théorie associée au recalaqe élastique On va maintenant décrire plus en détail la théorie associée à l'étape de recalage élastique. Les explications qui vont suivre sont faites en relation avec le recalage élastique de l'arbre vasculaire cérébral, mais il est bien entendu que ce recalage peut s'appliquer à d'autres types de tissus organiques. Un recalage élastique est une applications qui minimise une certaine erreur de recalage, ou plutôt, une énergie de recalage définie entre un ensemble de points sources, S et un autre ensemble de points destination, D. L'énergie du recalage R est la mesure d'une similarité entre l'ensemble de points transformés R(S) et l'ensemble cible D. Lorsque les deux ensembles se confondent alors cette énergie devient nulle et le recalage est accompli. Si l'énergie de recalage et l'espace de fonctions dans lequel évolue la transformation R ne sont pas définis correctement, alors le problème consistant à trouver un R optimal est mal posé, i.e. l'existence et l'unicité de R ne sont pas garanties.

Nous présentons dans cette partie, le cadre de calcul de la déformation élastique qui permet d'établir de façon robuste une correspondance élastique entre les données considérées. Les spécificités de cette approche résident dans les propriétés mathématiques de la transformation R, démontrées plus loin, qui en garantissent le réalisme physique.

Enerqie de recalage L'énergie de recalage peut être librement définie afin de refléter la nature du problème de recalage traité. Par exemple, dans le cas de recalage entre deux images, cette mesure peut être basée sur les différences entre les valeurs des pixels, ou voxels, des images source et destination. Dans le cas présent, la similarité entre les jeux de données S et D est définie par la proximité spatiale des deux nuages de points. Sans autre connaissance de la nature des données manipulées, il est possible, par exemple, définir une énergie de recalage Ereg comme la moyenne des distances euclidiennes minimales entre R(S) et D. 1 )+ s Où : - R représente la fonction de recalage élastique, - Card(S) est le cardinal de l'ensemble source S, - Card(D) est le cardinal de l'ensemble destination D, - d(R(s),D) représente la distance entre l'ensemble D et le transformé par R d'un point s de l'ensemble S, - d(d, R(S) représente la distance entre un point d de l'ensemble D et la transformée par R de l'ensemble S. Dans ce cas un recalage optimal est réalisé lorsque la valeur de Ereg s'annule. L'énergie symétrique présentée ci-dessus suppose que chaque structure présente dans S peut trouver, au moyen d'une transformation R appropriée, une correspondance parmi D, et réciproquement. Les désignations S et D des deux jeux de données sont interchangeables. Dans le cas présent, toutefois, les structures présentes dans S et D peuvent ne pas être équivalentes et l'hypothèse de symétrie ne tient pas puisque les structures dans S ne représentent qu'un sous-ensemble de D. Une définition alternative doit être trouvée pour en tenir compte. Dorénavant, on désignera par S l'ensemble de points définissant la configuration courante, per-opératoire, de l'arbre vasculaire, et D ceux qui en définissent la configuration de référence, préopératoire. La fonction de recalage élastique R est donc celle qui amène les points de S sur D et une définition asymétrique de l'énergie de ce recalage peut être formulée en tronquant l'expression symétrique ci-dessus : D 1i La signification des appellations source et destination peut être mieux saisie en voyant le recalage élastique comme un processus itératif où les points source S sont déplacés vers l'ensemble de destination selon une trajectoire optimale sur une surface d'énergie potentielle dont le relief est défini par la configuration spatiale des points de D. L'évaluation d'une configuration candidate R(S) est faite sur l'ensemble contenant moins de données tandis que la fonction d'énergie est plus strictement définie par D, lequel contient l'information de référence sur l'ensemble de la structure de l'arbre vasculaire. Des solutions alternatives au problème de l'asymétrie figurent dans la littérature. La définition ci-dessus de Ereg peut aboutir à des situations où un groupe de points sources est recalé vers une partie confinée de D, en d'autres termes, cette énergie ne favorise pas la répartition des points sources sur la structure destination. Pour éviter cela, il est possible d'en étendre l'expression en rajoutant un terme qui quantifie cette répartition mais qui opère uniquement dans un voisinage restreint de points de D qui contiennent des points de R(S). Finalement, étant donné que le cardinal de l'ensemble S n'évolue pas au cours du recalage, il est encore possible de simplifier l'expression précédente de l'énergie Ereg en éliminant la constante multiplicative qui la précède, ce qui nous amène à la définition suivante de Ereg : L[Z '.{; :.L) Espace des fonctions de déformation Etant donné la nature éparse et irrégulière des données, la recherche d'une fonction de recalage élastique doit être contrainte afin de ne pas aboutir à un recalage aberrant. Cet écueil est principalement dû aux nombreux minima locaux qui jalonnent le produit de la fonction d'énergie pour une configuration D donnée. Les deux ensembles S et D sont des observations du phénomène physique étudié i.e. la déformation du cerveau.

On suppose donc que de façon similaire au processus du brain-shift , sans déchirement ou résection de tissu, la transformation non-linéaire R cherchée est continue et même continuement différentiable i.e. différentiable en tout point du domaine et ne présentant pas de discontinuité de différentielle. Ceci définit une première restriction sur l'espace de recherche des fonctions de recalage : Une autre contrainte forte sur la nature de la transformation est qu'elle ne doit pas replier l'espace. La contrainte de non-repliement peut s'exprimer mathématiquement comme suit. Considérons un trièdre orienté positivement constitué de trois vecteurs infinitésimaux, dX1, dX2 and dX3, placés en un point p de l'espace. Soit dV le volume, positif, défini par ces trois vecteurs et dont la valeur est donnée par le déterminant de la matrice formée par les trois vecteurs. dxi Après application de la transformation R au point p les trois vecteurs se transforment en dxl, dx2 et dx3. La nature infinitésimale de ces vecteurs nous autorise l'approximation suivante : Le volume défini par le trièdre transformé, noté dv est alors obtenu de façon similaire. C"1'?+ r ... La contrainte de non-repliement s'énonce alors tout simplement : dv > 0. Ce qui signifie que le volume infinitésimal déformé n'a pas été retourné (dv < 0) par R ou totalement écrasé (dv = 0). Ainsi l'application R ne replie pas localement l'espace si et seulement si son jacobien en p est strictement positif. Dans notre cas le recalage élastique ne doit replier l'espace en aucun point, ce qui nous amène à l'expression de la contrainte de non-repliement : Si en un point donné, le jacobien est plus grand que 1 cela signifie que le recalage élastique étire localement le volume, si au contraire il est plus petit que 1, alors le volume est contracté. Enfin le produit de la matrice jacobienne de R calculée en un point par un vecteur unitaire nous donne le taux d'allongement ou de contraction de l'espace dans une direction donnée au voisinage de ce point. La figure 3 montre une déformation d'un carré unitaire avec un repliement d'espace. Finalement pour les mêmes motivations physiques, l'application R doit être injective. Ce qui se traduit par le fait que deux points distincts ne peuvent pas avoir la même image par R et, physiquement, R ne doit pas amener deux atomes en un même point de l'espace. Si R est une injection, alors pour chaque point atteint P ), un prédécesseur peut être calculé. Dans le cas présent où S représente une fraction de la structure présente dans D, la transformation cherchée in fine est celle qui amène D sur S, soit R-1, afin de définir le champ de déplacements de l'arbre vasculaire de la configuration de référence D vers la configuration courante S, comme décrit ci-dessus. Pour que l'expression R-' ait un sens, il faut que 1 _ :În' Y'• ('" . Afin de simplifier ces considérations, on considère les fonctions R surjectives, i.e. telles que i .4' i(' . Il est donc possible de conclure que la fonction de recalage R cherchée est bijective.

La contrainte de non-repliement garantit que le jacobien est strictement positif en tout point de -ceci est cependant insuffisant pour garantir la bijectivité globale. L'espace de fonctions doit alors être convenablement défini afin d'assurer la bijectivité de la fonction R. Les sections suivantes décrivent la stratégie de recalage adoptée afin d'aboutir à une minimisation de l'énergie de recalage tout en respectant les contraintes énoncées plus haut de différentiabilité, non-repliement et bijection.

Calcul du recalage L'approche présentée ici est inspirée de la modélisation mécanique des solides basée sur une discrétisation de type éléments finis. Afin d'approximer le comportement mécanique des données non-structurées, on suppose que le nuage de points soumis à la déformation est englobé dans un solide virtuel élastique. Puisque le déplacement des points traduit les propriétés physiques d'une structure sous-jacente dont la géométrie est inconnue, on attribue arbitrairement à ce solide virtuel la forme d'une boîte englobante des points. On définit alors l'ensemble des déformations admissibles comme les combinaisons successives de déformations élémentaires du solide englobant, dont les effets se propagent à travers son volume et modifient la position des points internes. Les déformations élémentaires mentionnées ici sont définies plus loin. La technique de recalage 3D décrite ici peut aisément être dérivée en une technique de recalage 2D en remplaçant la notion de solide virtuel par plan virtuelle . Les illustrations présentées ci-après se basent sur une implémentation en 2D du recalage pour en clarifier les concepts. La fonction de recalage élastique est assemblée itérativement de manière à ce que l'énergie de recalage Ereg diminue de façon optimale à chaque pas. Soit So := S, le nuage de points source initial, et soit Si le nuage de points source obtenu suite à l'itération i. Le domaine du solide virtuel est discrétisé avec un certain niveau de raffinement en hexaèdres réguliers (ou carrés, en 2D) appelés cellules . Les déformations élémentaires du solide sont obtenues en déplaçant individuellement les noeuds de la grille ainsi formée. L'union des cellules entourant un noeud n définit son voisinage , noté V (n). Une fonction de recalage élémentaire r est définie en propageant le déplacement du noeud n à l'ensemble des points de sous les contraintes suivantes : -La Front i( hl ii oisiElage le Tl, f 'f(ii) est. 1 n ilaIi ;Ë[-' Le. 1F2ti7i1 '(' Chi 1Y FC'yff? (.~silÏf?Illi°tif, 111 Ilr:i11 tir' 2,('. 1 nts !('y 1)i)ilir s f I tFi'2`]cur (' vois 1 T'(sf (Tu iFlcii=117g(.' 1 ~t11). )(p = j] Chaque noeud du solide discrétisé est considéré et le gradient de l'énergie de recalage, en tant que fonction du déplacement du noeud, est calculé. L'opposé de ce vecteur définit le déplacement préférentiel du noeud, c'est-à-dire celui qui, appliqué au noeud, entraîne la diminution la plus importante de l'énergie de recalage. De tous les noeuds, celui qui permet la diminution la plus importante du gradient est choisi et son déplacement préférentiel lui est appliqué tandis que les autres noeuds demeurent immobiles. La déformation locale du solide est calculée en propageant le déplacement du noeud à travers son voisinage et la position des points source qu'il contient est modifiée en conséquence, donnant ainsi la nouvelle configuration du nuage de points Si+1.

Une fois la déformation élémentaire appliquée, le solide retourne dans sa configuration initiale, au repos, et les points source Si+1 sont distribués dans son volume. Une nouvelle itération peut alors être calculée en prenant la nouvelle configuration Si+1 comme point de départ. Les itérations s'arrêtent, enfin, lorsqu'aucune diminution significative de l'énergie ne peut être obtenue en déplaçant les noeuds de la grille du solide. La boucle itérative décrite ci-dessus suppose un degré de raffinement du domaine fixe. Afin de maintenir la consistance spatiale de S au cours de l'assemblage de R et en éviter le morcellement, une approche multi-échelle a été adoptée. L'assemblage itératif de R commence au niveau de raffinement le plus grossier et lorsque l'énergie de recalage ne peut être améliorée davantage à un niveau donné, les cellules sont subdivisées en 8 sous-cellules (ou 4, en 2D) à la manière d'un octree.

Les itérations de raffinement s'arrêtent lorsque le degré maximal de subdivision spécifié manuellement a été atteint. La dimension de la plus petite cellule atteinte au cours du processus de raffinement définit les dimensions de la plus petite structure présente dans S et qui ne sera pas recalée correctement si une structure équivalente est absente de D.

La figure 4 illustre en 2D l'approche itérative multi-échelle décrite plus haut. Les points source S sont représentés par les points référencés 1. L'ensemble D n'est pas représenté, pour plus de clarté. En (a) l'ensemble initial SO := S est placé au sein du solide virtuel (référencé 2) discrétisé par une cellule (niveau 1) et les gradients de l'énergie de recalage sont calculés en ses quatre noeuds 3. En (b) le meilleur déplacement nodal est appliqué au solide et les points source sont déplacés en S1 (référencé 4). Le solide retourne alors à sa configuration initiale (a). Après plusieurs itérations, si aucune amélioration significative de l'énergie ne peut être obtenue, le niveau de raffinement de la grille augmente (c). Les gradients de l'énergie sont calculés aux neuf noeuds 5 de la nouvelle grille. En (d) le meilleur déplacement est appliqué au solide et la nouvelle configuration de points source S2 (référencé 6) est calculée. Grâce à la technique de remise à zéro de la déformation du solide après chaque itération, de grandes déformations de S peuvent être envisagées sans avoir à maintenir une grille irrégulière, soumise à une somme de déformations importantes ce qui aurait été le cas si le solide suivait strictement les déformations de S. La régularité de la grille avant chaque itération permet également un gain substantiel d'effort calculatoire, étant donné que l'interpolation des déplacements nodaux peut être calculée sur la base d'une cellule type cubique. Avant de détailler davantage le calcul de R, nous devons étendre le concept intuitif du solide virtuel de manière à ce que chaque déformation élémentaire soit définie comme une application R., }, afin de répondre à la restriction énoncée plus haut. Considérons pour cela un point source placé à la frontière du solide virtuel par exemple à proximité du bord supérieur du carré sur la figure 4. Si le noeud est déplacé comme cela est illustré en (d), ce point va tomber en dehors du solide dans sa configuration au repos (carré) et deviendra un point mort , dont la position ne pourra être modifiée par les itérations suivantes, lesquelles n'ont d'effet que sur les cellules de 0.

Pour pallier cet inconvénient, on ajoute autour du domaine, une marge de cellules de dimensions identiques à celles utilisées pour la discrétisation. Cet ajout est réalisé à chaque niveau de raffinement. Les noeuds extérieurs de ces cellules sont fixes et assurent ainsi un rôle d'interface entre la zone déformable et le reste de ' , qui demeure inchangé. La figure 5 montre les cellules supplémentaires 7 insérées en marge de la discrétisation du domaine pour deux niveaux de raffinement successifs.

Propriétés du recalaqe élastique Soient G1, ...,GJ les niveaux de raffinement successifs, G1 étant le niveau le plus grossier et GJ le plus fin. Soit l'ensemble des déformations élémentaires calculées à chaque niveau de raffinement j (où lj est le nombre total de déformations élémentaires au niveau j). L'expression de la fonction de recalage élastique R assemblée itérativement est : p Il est clair que les propriétés de différentiabilité, non-repliement et bijection de R dépendent des propriétés des fonctions de recalage élémentaires r'j. Si chaque r'; est une 5lé alors la fonction de recalage R est également une bijection et son bijection inverse est donnée par : Si chaque r'; est une application continuement différentiable alors la fonction de recalage R l'est également. Soitil un point quelconque et pk le point transformé après application de k déformations élémentaires. Soit K = 11 + 12 + ... + IJ , alors 1. ' 1. est l'ensemble de toutes les positions successives de p0 et R(p) = pK. La différentielle de R en p0 est donnée par la règle de composition suivante : I)I' j) , ï..Y (:)p Finalement, si tous les r', vérifient la condition de non-repliement, alors cela est vrai également de l'application R. Ceci est la conséquence de l'expression du jacobien de R : la: J r iep .

Ainsi pour prouver que toute fonction R de la famille de fonctions de recalage élastique générée à l'aide de la procédure décrite plus haut vérifie les contraintes de différentiabilité, bijection et non-repliement, il suffit de prouver que ces propriétés sont bien vérifiées pour chacune des fonctions de recalage élémentaire r', qui composent R.

Etude des fonctions de recalage élémentaire Expression générale de r, A ce point, la famille de fonctions de recalage élémentaires doit être spécifiée. Dans la méthode des Eléments Finis, une valeur nodale (déplacement ou autre quantité) est interpolée à travers le volume discrétisé au moyen de fonctions de forme '1' ., 1]. La valeur de la fonction de forme associée à un noeud et évaluée en un point du domaine donne la proportion de la valeur nodale transférée en ce point. Ainsi la valeur de w en son noeud est 1. Le même mécanisme pour propager le déplacement d'un noeud à travers son voisinage va être utilisé. On considère la ième itération de recalage effectuée au niveau de raffinement j et mettant en oeuvre la fonction de recalage élémentaire r',. Les indices d'ordre et de niveau de raffinement i et j seront omis dans la suite de ce paragraphe. Soit le déplacement et w la fonction de forme associés au noeud dont le déplacement est réalisé lors de cette itération. L'expression de r est donc : p+?{1~r, Les conditions sur la zone d'influence de r se traduisent par des contraintes sur le support de w de la façon suivante : 21 éfl in).,t; , i;,p) _ 0 Remarquons que le support de w est i: n ', qui est un compact Différentiabilité La différentiabilité de r découle de celle de w. 5 alors et la matrice jacobienne de r s'écrit : i u.). (ir 1 10 où Id3x3 est la matrice identité 3 x 3. Donc pour que r soit , {.. il suffit que w le soit aussi.

Non-repliement En réécrivant la contrainte de non-repliement à l'aide de l'expression de la matrice jacobienne de r et en développant le déterminant de la matrice, une nouvelle expression de la contrainte de non-repliement est obtenue : 1) M .0 le point où I (I' Ilatteint son maximum M - Il P s. i) Alors l'équation précédente est équivalente à : Le non-repliement de r peut donc être garanti en limitant l'amplitude des 20 déplacements nodaux. La norme maximale du déplacement nodal est donnée par l'inverse de la norme maximale du gradient de la fonction de forme considérée w.

Bijection Soit V := V (n) l'union fermée des cellules entourant le noeud n. Dans cette partie 25 nous allons montrer que si une application élémentaire r, ne replie pas l'espace, alors elle est également bijective de V dans lui-même mais également de R3 ! R3, puisque les points à l'extérieur de V ne sont pas affectés par r.

Implémentation , Puisque p ~ 'y' u.> l est une application continue à support compact, soit t 15 Choix d'une expression de w Afin de simplifier le calcul du recalage élastique, une expression polynomiale de w est choisie. Soit ° un polynôme de degré 3, défini par : .( = t.) `=1 -7- Son expression est Une fonction de forme w définie sur la cellule canonique [0' I1 et prenant la valeur 1 au noeud n = 1 1_ . 1 peut ainsi être obtenue comme : 'äTr ) 77. p2 La figure 6a montre la fonction de forme restreinte à une dimension, i.e. lf, ù sur l'intervalle [0, 1]. La fonction de forme w est définie sur l'union des cellules voisines autour du noeud n par un changement de variable et une mise à l'échelle afin d'adapter son domaine de définition canonique [t h1.3 aux dimensions des cellules dans la grille de discrétisation du solide virtuel. La figure 6b montre les changements de variable nécessaires pour assembler à partir de la fonction de forme canonique, la fonction de forme w sur un voisinage de 2x2 cellules du plan. La figure 6c montre le tracé de la fonction de forme ainsi obtenue w R2 ù .

Contrainte de non-repliement et de bijection L'expression de w étant choisie, il est possible de calculer l'amplitude du 20 déplacement nodal maximal autorisé pour une fonction de recalage élémentaire r. Pour définir cette amplitude, une borne supérieure M est calculée pour la quantité 1 'Tl sur le domaine R La dérivée de ù Elle atteint son maximum sur [0, 1] en tmax = 1/2 où ~ -r . On obtient l'expression de : V SurL``' on a . Cette borne supérieure est valide pour chaque coordonnée, ainsi et donc 11-7.' Il Etant donnée cette borne supérieure pour la norme du gradient, une condition suffisante qui garantit la bijection et par là même, le non repliement de l'espace pour des fonctions de déformation élémentaires définies pour une grille composée de cellules canoniques [0, 1] est donc : Ce résultat peut être étendu à n'importe quelle dimension de cellule de la façon suivante. Soit + ] + Li ..r [ + Li une cellule cubique Ti ., d'arête L. Soit la fonction R' qui transforme cette cellule en cellule canonique , définie par : La fonction de forme w définie sur Oeil est ,3.i .1 ) ;(] et: De la même façon que pour la cellule canonique, il est possible de conclure qu'une borne supérieure pour la norme du gradient est : Ce qui conduit à la condition suffisante garantissant pour une cellule de dimension L la bijection et le non-repliement de la déformation élastique élémentaire r : Etant donné que l'approche présentée ici est basée sur la modélisation mécanique des solides, on choisit de restreindre l'amplitude des déplacements individuels des noeuds au domaine dit des petites déformations . Pour ce faire on limite les amplitudes des déformations imposées aux cellules du solide à 10% de leur dimension, soit : Une conséquence de cette restriction est que la compression et l'étirement de alors : l'espace appliqués au solide à chaque itération sont bornés. Puisque .. y } 20 L'équation ci-dessus fournit une approximation i.-., y1/20 0.25980 .. des bornes inférieure et supérieure du jacobien de r valable pour toutes les dimensions de cellules : 0.71 . ,: 1..2( Minimisation de l'énergie On considère deux ensembles S et D ayant les cardinaux nS = Card(S) et nD = Card(D). La minimisation de l'énergie de recalage Ereg nécessite le calcul de la distance minimale entre chaque point de S et l'ensemble D. La complexité du calcul brutal, considérant pour chaque point de S sa distance à chaque point de D est s est nS x nD. De nombreuses techniques ont été proposées pour accélérer le calcul de la distance la plus proche entre deux ensembles de points, principalement dans le cadre du recalage rigide basé sur le principe de l'ICP pour Iterative Closest Point .

Afin d'accélérer le calcul de Ereg à chaque itération, une carte de distance euclidienne est calculée à partir des points de D. La carte de distance est calculée sur une grille 3D discrétisant l'espace autour de D. A chaque sommet de la grille gi, la distance minimale à D, d(gi,D) est mémorisée. Cette structure de données est initialisée lors de la phase préopératoire, juste après la segmentation de l'arbre vasculaire D dans les images ARM.

Le contenu de la carte de distance est ensuite enregistré sous forme de fichier binaire pour être lu au début de l'intervention par le système, avant tout calcul de recalage élastique. Pour un point source donné s, la valeur de d(s,D) est calculée par interpolation 5 trilinéaire parmi les 8 valeurs lues dans la grille de la carte de distance, associées aux sommets qui entourent s, {gi}i=sl ,...,s8. Lors du recalage, à chaque itération, le meilleur déplacement pour chaque noeud de la grille qui discrétise le solide virtuel est calculé à l'aide du gradient de Ereg. Le calcul du déplacement nodal est réalisé en deux étapes : 10 1. Le gradient de Ereg par rapport aux déplacements nodaux est estimé par différences finies ; 2. Une recherche le long de la direction de plus forte descente est effectuée afin d'estimer le pas optimal ; 3. L'amplitude du déplacement ainsi obtenu est limitée de manière à satisfaire une 15 contrainte. Etant donné que l'on ne considère que les déplacements d'un noeud à la fois, n, seules les trois coordonnées du vecteur de déplacement nodal sont impliquées dans le calcul de la minimisation. La fonction minimisée est '' définie par :

20 où r P H I + WW,` est la fonction de recalage élémentaire associée au noeud n, et R est la fonction de recalage global assemblée avant l'itération courante. La stratégie consistant à revenir à la forme initiale du solide virtuel à l'issue de chaque itération permet de supposer que l'ensemble des points source au début de l'itération a déjà été transformé par R selon l'algorithme décrit ci-dessus Ceci nous permet 25 de laisser de coté l'historique du recalage que constitue R. Soit `~ !, l'ensemble de points source dont la position appartient au support de la fonction de forme w associée au noeud du maillage n. Puisque ce sont les seuls points concernés par le déplacement de n à l'itération de recalage courante, le calcul de l'énergie de recalage peut être restreint à cet ensemble. Ceci nous permet 30 d'exprimer la fonction à minimiser f comme : f ( t; = La fonction f n'est pas systématiquement différentiable par rapport à ?' , du fait de la racine carrée portée par la distance, et de l'opérateur min(x, y) qui opère sur les distances entre les nuages de points. La direction de plus grande descente doit donc être estimée localement par différences finies centrées. On obtient ainsi : où est une valeur très petite par rapport à l'amplitude maximale du déplacement possible de n, qui est inférieure à 10% de la taille de la cellule L, d'où ). Les vecteurs " F '.= 1-.2 quant à eux, forment la base canonique de Le vecteur de plus grande descente sur f peut s'écrire : La longueur du pas maximal autorisé par la contrainte de non-repliement est 10% de L, ce qui permet de définir le vecteur maximal de déplacement du noeud n dans la direction de descente surf : Le meilleur déplacement nodal au sens de la minimisation de l'énergie de recalage correspond donc au minimum de la fonction t i ` + f zt. ,,.r, (t ,: ;:< pour i [0.1] Une estimation convenable de ce minimum est obtenue en approximant sur l'intervalle [0, 1] la fonction t ` + par une parabole 1 ° t I . _Ï + + Les trois coefficients de la parabole A, B et C sont déterminés par les informations 20 locales sur le comportement de f ainsi qu'une évaluation de sa valeur au pas maximal : nia. Ainsi . L r..3.r_E = t. t. (_F B = { Ilv Pi( l) A + t_3 + f: % ù i'4 = ! \ Y! i:lzl"\r + L110 Il ,t{,1 Deux situations sont alors possibles. Si A > 0 alors P est une parabole convexe et son minimum est atteint en tmin = -B/(2A). Dans ce cas, si tmin > 1 alors la valeur du pas est limité à 1. Si au contraire, A < 0 alors P est une parabole concave et son minimum sur est atteint en tmin = 1. La figure 7 montre un exemple représentatif d'approximation de la courbe d'énergie de recalage 11 par une parabole P 12 sur l'intervalle [0, 1]. Des mesures d'écart entre les deux courbes, effectuées sur un grand nombre de jeux de données lors du déroulement de la procédure de recalage montre les qualités de cette approximation. La courbe du bas 11 est obtenue en échantillonnant l'énergie de recalage sur l'intervalle [0, 1], celle du haut 12 est le tracé de la parabole P. Malgré un léger écart entre les valeurs des deux courbes, on remarque que les valeurs de t pour lesquelles ces courbes atteignent leur minimum sont quasiment identiques. Pour un niveau de raffinement de la discrétisation du solide virtuel donné, la minimisation de l'énergie est ainsi menée, pas à pas, jusqu'à ce qu'aucune diminution significative ne puisse être obtenue par application d'un déplacement à un des noeuds du maillage.

Dans ce cas, si le niveau de raffinement maximal n'a pas été atteint, la grille est raffinée et la même procédure est répétée. Inversion Le calcul de l'inverse de la fonction de recalage élastique-Iiva maintenant être décrit. Ce calcul nécessite l'inversion successive de toutes les fonctions de recalage élémentaires : (r';)-'. Etant donné qu'aucune expression analytique ne peut être donnée pour (r';)-1, le calcul de (r';)-' (q) est réalisé individuellement pour chaque point considéré q. Pour simplifier les écritures, les indices i et j ne seront pas précisés dans la suite, et on considéra plus généralement une fonction élémentaire r, associée au déplacement "' 25 du noeud n. Pour un point donnée 1' tel', son antécédent p = r-1(q) est calculé comme suit. On a vu précédemment que p peut s'écrire comme P 1 +Ç ? où tp est la racine de la fonction polynomiale f sur l'intervalle [ta, tb]. En fait, puisque (1 ii'tp l'intervalle de recherche peut être restreint à [max{-1, ta}, 0]. 30 Le choix de l'expression de la fonction de forme c1.' ( . - :T ' amène à l'expression suivante de f : tI- t+ I+ t, )z_ 3) est un polynôme de degré 3 en t, donc f est de degré 9. Il n'y a pas de formule explicite permettant de trouver les racines d'un polynôme de degré 9 aussi doit on faire appel à une technique itérative de type Newton. En supposant que la valeur de ta est déjà initialisée à max{-1, ta}, on commence 5 avec l'estimation de la solution suivante : Les itérations sont ensuite menées de manière classique : L'approximation de ù,k à chaque itération est donnée par l'erreur correspondante dans l'espace déformé est 1171. La précision de l'approximation de l'inverse est mesurée dans l'espace initial, non déformé, en évaluant la quantité M. p n'est pas connu, ce qui rend difficile l'évaluation cette erreur. Afin d'évaluer 11P P , on s'appuie sur le fait que l'application 15 inverse r-' est k-lipschitzienne, en d'autres termes : 10 q 'ù ù ; .-1~q - Cette propriété permet alors de majorer l'erreur I I par : IIP- pII=Il ll _ Ill t,ll = 20 On désire poursuivre les itérations de l'algorithme de Newton jusqu'à ce que la

précision voulue soit atteinte i.e. lorsque 11f)' r •. Au vu de ce qui précède, sf'; cette condition d'arrêt est vérifiée lorsque Y p., t' '# 3 Cette inégalité constitue un critère d'arrêt fiable défini dans l'espace non-déformé. Considérons maintenant le cas d'une fonction de recalage R composée de 25 plusieurs fonctions élémentaires rj . Dans ce cas, la précision du calcul de l'inverse

' doit être contrôlée à chaque étape d'inversion d'une fonction élémentaire.

Voyons comment nous pouvons évaluer la précision du calcul de R-1 à l'aide des -1 constantes de Lipschitz des fonctions :{ Pour plus de clarté, on suppose que On veut évaluerR ,, en ayant, pour chacune des , une constante de Lipschitz kj . Soit P l'antécédent cherché de p3 = R(p) ayant pour images intermédiaires p1 = r1 (p) et p2 = r2(r1(p)) = r2(p1). La figure 8 illustre, de droite à gauche, les trois étapes de l'approximation de'~` # ainsi que l'accumulation successive des erreurs commises à chaque pas : 1. q2 approxime '' r > 'P3'. L'erreur dans l'espace déformé est _3 := kq3-p3k. Cette erreur se propage vers l'espace non-déformé comme et finalement 11(1 )11 2. r1 approxime ~i - avec une erreur dans l'espace déformé de(, .- llr-- .rll Cette erreur se propage vers la gauche, i.e. vers l'espace non-déformé, comme puis r - q 3. Finalement s approxime avec une erreur dans l'espace déformé := ù 115, . _, de . Cette erreur se propage comme 118 L'erreur finale peut ainsi être majorée par la somme des erreurs propagées : =Ils ù II .- Ils ù III+III {III11(- + L'expression ci-dessus peut être immédiatement étendue à toute fonction de recalage élastique R. Si R est composée de N fonctions élémentaires et la précision voulue dans l'espace non-déformé est (7 alors l'effort d'approximation peut être réparti parmi les N fonctions élémentaires en ajustant le seuil de précision de l'inversion de 25 Newton à :

où chaque constante de Lipschitz ki ne doit être évaluée qu'une seule fois étant donné qu'elle ne dépend que de l'amplitude du déplacement appliqué au noeud associé à la fonction élémentaire ri et de la taille de la cellule au niveau de raffinement considéré.

Filtraqe du recalaqe élastique et calcul du champ de déplacements Comme décrit ci-dessus, les images per-opératoires sont bruitées et le recalage élastique peut produire occasionnellement des artefacts de recalage. Il est donc important de mettre en place une stratégie permettant de réduire le nombre d'association entre points source et points destination erronées, chacune de ses associations se traduisant par la suite par un déplacement, leur effet risque de se propager jusqu'au calcul de la déformation biomécanique du modèle avec pour conséquence, une perte de précision importante. Une technique robuste et simple a donc été implémentée pour éliminer les erreurs de recalage grossières, dues à la présence de faux vaisseaux dans les données issues de la segmentation. Cette méthode est divisée en deux phases : un filtrage direct, qui s'appuie sur le calcul de la fonction de recalage R, suivi par un filtrage inverse, lequel fait intervenir l'inverse du recalage : R-'. Le rôle du filtrage est d'éliminer les associations de points qui ne vérifient pas un critère de précision donné. Soit dmax la distance maximale acceptable entre un point source et son homologue dans D. Le filtrage direct permet d'éliminer les points de S qui ne remplissent pas la condition d(R(s),D) < dmax. L'ensemble de points source filtré est Sf , défini par :

La figure 9a illustre l'étape de filtrage direct. Comme on peut le constater dans 25 cette figure, les points sources n'ayant pas atteint D sont éliminés. L'objectif de la procédure de recalage est de générer un champ de déplacements correspondant à la déformation de la structure de l'arbre vasculaire (dans la présentation faite ici, mais plus généralement du tissus organique déformé) définie à partir des images préopératoires par D, et qui l'amènent dans la configuration courante S, obtenue par 30 analyse des images US Doppler per-opératoires. Le filtrage inverse permet de vérifier que le champ de déplacements obtenu par l'inversion de la fonction de recalage R transforme effectivement les points de D en S. Afin de ne pas tenir compte des points source éliminés précédemment, nous restreignons l'ensemble S à Sf et définissons l'ensemble de points destination filtré Db par : Ta.y. La figure 9b illustre l'étape de filtrage inverse. Comme l'illustre cette figure, les déplacements de D qui n'atteignent pas S sont éliminés.

Finalement, l'ensemble de déplacements de l'arbre vasculaire, de sa configuration de référence vers sa configuration courante, est assemblé en s'appuyant sur ce dernier ensemble de points. Le champ de déplacements destiné au modèle biomécanique est ainsi donné par : Le recalage élastique décrit en détail ci-dessus permet de construire une association entre deux observations de l'arbre vasculaire cérébral pré- et per-opératoire, de natures très distinctes. Le calcul est guidé par la minimisation de l'énergie de recalage qui quantifie la similarité entre les deux structures. Le recalage optimal étant atteint pour une énergie nulle, l'objectif de l'algorithme de recalage est d'en minimiser itérativement la valeur en faisant évoluer la structure courante, dite 'source', vers la structure de référence, dite 'de destination'. La rapidité du calcul est obtenue grâce au calcule préalable d'une carte de distance sur laquelle s'appuient l'évaluation de l'énergie de recalage ainsi que la détermination de la direction de plus forte descente, calculées à chaque itération du processus de minimisation. Les propriétés mathématiques de la famille de fonction de recalage engendrée par la méthode décrite ci-dessus garantissent le respect de critères physiques essentiels, tels que le non-repliement de l'espace, le non-recouvrement de matière et la continuité de la propagation de la déformation. Ce recalage élastique présente également l'avantage d'être inversible avec le degré de précision souhaité. Sur la base de cette dernière caractéristique, un post-traitement du résultat produit par le recalage est proposé. Ce post-traitement permet d'éliminer les aberrations dues à la présence de bruit dans les données d'imagerie per-opératoire. L'ensemble de déplacements de l'arbre vasculaire finalement retenu constitue la donnée d'entrée du modèle biomécanique dont le rôle est d'en étendre la portée à l'ensemble du volume du cerveau du patient afin de simuler l'effet du brain-shift. On a illustré à la figure 2 un mode de réalisation du système pour la mise en oeuvre du procédé décrit ci-dessus.

Le système comprend des moyens de localisation 10, des moyens d'acquisition 20, des moyens de traitement 30 et des moyens d'affichage 40. Les moyens de localisation permettent la localisation dans l'espace des instruments utilisés par l'utilisateur lors de l'intervention chirurgicale. Ces moyens de localisation peuvent comprend par exemple une caméra stéréoscopique et des marqueurs associés aux différents objets à localiser. Les moyens de localisation 10 permettent également la localisation dans l'espace de la tête du patient lors de l'intervention, notamment pour le repositionnement spatial de l'arbre vasculaire de référence calculé à partir de l'angiographie par résonance magnétique.

Les moyens d'acquisition 20 permettent l'acquisition des images bidimensionnelles pendant l'opération. Ces moyens d'acquisition 20 comprennent par exemple une sonde à ultrasons localisable par les moyens de localisation 10. Les moyens d'affichage permettent par exemple l'affichage des données préopératoire mise à jour. Ces moyens d'affichage peuvent comprend un écran.

Les moyens de traitement permettent de mettre en oeuvre les différentes étapes de calcul du procédé, et notamment : - l'étape de traitement 100 de l'image tridimensionnelle du cerveau du patient, acquise avant une intervention chirurgicale, pour obtenir une arborescence artérielle cérébrale de référence du patient, - l'étape de traitement 200 des images bidimensionnelles du cerveau du patient, acquises pendant l'opération, pour reconstituer au moins partiellement l'arborescence artérielle cérébrale courante du patient, - l'étape de détermination 300, à partir de la mise en correspondance des arborescences artérielles cérébrales de référence et courante, du champ de déformation de l'arbre vasculaire représentant le déplacement de l'arbre vasculaire courant par rapport à l'arbre vasculaire de référence, et - l'application 400 du champ de déformation de l'arbre vasculaire déterminé à un modèle biomécanique du cerveau du patient pour 30 estimer le déplacement du cerveau du patient. Ces moyens de traitement peuvent être un ordinateur et comprendre des moyens de commande pour le contrôle des différents éléments matériels du système. L'homme du métier aura compris que de nombreuses modifications peuvent être apportées au procédé décrit sans sortir matériellement des nouveaux enseignements 35 présentés ici. Notamment, le procédé de recalage décrit précédemment en relation avec le recalage de l'arbre vasculaire cérébral peut être mis en oeuvre pour permettre le recalage de tout type de tissu organique ayant subit une déformation. Par ailleurs, même si cette étape de recalage a été présentée en relation avec des ensembles de points sources S et destinations D présentant des nombres de points différents, l'étape de recalage peut être mise en oeuvre sur des jeux de données S et D comportant le même nombre de points. Il est donc bien évident que les exemples qui viennent d'être donnés ne sont que des illustrations particulières en aucun cas limitatives.

Claims (7)

1. REVENDICATIONS1. Procédé pour le recalage d'un ensemble de points sources (S) dans une image de référence représentant un tissu organique et d'un ensemble de points destinations (D) dans une image courante représentant le tissu organique déformé, le procédé comprenant la détermination d'une fonction de recalage (R) par minimisation itérative d'une énergie de recalage Ereg, la fonction de recalage correspondant à une transformation permettant de passer des points sources (S) de l'image de référence aux points destinations (D) correspondants dans l'image courante, caractérisé en ce que la fonction de recalage (R) est une fonction de recalage élastique satisfaisant des conditions pour garantir le respect de critères physiques du tissu organique, les conditions comprenant : - le Jacobien de la fonction de recalage élastique est supérieur à 0, - la fonction de recalage est continûment différentiable, - la fonction de recalage est bijective.
2. 2. Procédé selon la revendication 1, caractérisé en ce que l'ensemble de points sources est de nature distincte de l'ensemble de points destinations.
3. 3. Procédé selon l'une des revendications 1 ou 2, caractérisé en ce que l'énergie de recalage Ereg de la fonction de recalage (R) satisfait la relation suivante : Ereg = E d (R(s), D) ses avec : - s un point de l'ensemble de points sources S, - R(s) le transformé du point s par la fonction de recalage R, - d(R(s), D) la distance euclidienne entre le transformé du point s par la fonction de recalage et l'ensemble de points destinations, noté D.
4. 4. Procédé selon l'une des revendications 1 à 3, caractérisé en ce qu'il comprend une étape de filtrage direct dans laquelle chaque point s de l'ensemble de points sources S dont l'association avec un point d de l'ensemble de points destinations D vérifie un critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points source S'.
5. 5. Procédé selon la revendication 4, caractérisé en ce que le critère de rejet est que la distance entre :o le transformé du point s de l'ensemble de points sources S par la fonction de recalage et o le point correspondant d de l'ensemble de points destinations D est supérieure à une valeur maximale prédéterminée.
6. 6. Procédé selon l'une des revendications 4 ou 5, caractérisé en ce qu'il comprend en outre une étape de filtrage inverse dans laquelle chaque point d de l'ensemble de points destinations D dont l'association avec un point s de l'ensemble de points sources S vérifie un autre critère de rejet est supprimé de sorte à obtenir un ensemble filtré de points destinations D'..
7. 7. Procédé selon la revendication 6, caractérisé en ce que l'autre critère de rejet est que la distance entre o le transformé inverse du point d de l'ensemble de points destinations D par la fonction de recalage inverse et o le point correspondant s' de l'ensemble filtré de points sources S' est supérieure à la valeur maximale prédéterminée. 9. Système pour le recalage d'un ensemble de points source dans une image de référence représentant un tissu organique et d'un ensemble de points destinations dans une image courante représentant le tissu organique déformé, caractérisé en ce qu'il comprend des moyens pour la mise en oeuvre du procédé selon l'une des revendications 1 à 7. 10. Produit programme d'ordinateur comprenant des instructions de code programme enregistré sur un support utilisable dans un ordinateur, caractérisé en ce qu'il comprend des instructions pour la mise en oeuvre du procédé selon l'une des revendications 1 à 7.