FR2936890A1 - Procede de traitement d'image pour estimer une deformation d'un cerveau d'un patient - Google Patents

Abstract

L'invention concerne un procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient, le procédé comprenant 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 du cerveau du patient et 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é.

Description

DOMAINE DE L'INVENTION

La présente invention concerne le domaine de la neurochirurgie assistée par ordinateur, et plus particulièrement un système et une méthode permettant la localisation précise d'une cible.

PRESENTATION DE L'ART ANTERIEUR

La localisation précise d'une cible est essentielle pour réduire le taux de morbidité 10 au cours d'une intervention chirurgicale d'ablation d'une tumeur cérébrale. Au cours d'une intervention chirurgicale sur le cerveau, une déformation des tissus cérébraux (ou brain shift selon la terminologie anglo-saxonne) peut se produire après qu'une ouverture ait été réalisée dans le crâne du patient. Cette déformation peut être liée à des phénomènes physiques (gravité, perte du 15 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 du cerveau du patient et de retrouver les données préopératoires du patient pendant la procédure chirurgicale, des systèmes de 20 neuro-navigation ont été développés. Ces systèmes de neuro-navigation permettent d'afficher au chirurgien la position des objets chirurgicaux par rapport à la position des caractéristiques anatomiques cérébrales du patient. Dans ce type de systèmes la compensation de la déformation du cerveau du 25 patient est en général réalisée sur le principe suivant. Une première image de la surface corticale du patient est acquise avant l'opération chirurgicale. Pendant l'opération chirurgicale, une seconde image de la surface corticale du patient est acquise. Des moyens de traitement du système comparent ces deux images et estiment un champ volumique de déformation qui permet d'aligner les points 30 dans les deux images afin de compenser la déformation du cerveau du patient. De tels systèmes sont notamment décrits dans les documents US 7 072 705 et US 2007/0021669. Toutefois, un inconvénient de cette technique de compensation de la déformation du cerveau du patient à partir de points de la surface corticale est qu'elle laisse une zone 35 importante du cerveau (sous la surface corticale) dans laquelle la déformation réel du cerveau est inconnue. 1 Un but de la présente invention est de proposer un système et un procédé de compensation de la déformation du cerveau du patient permettant de pallier au moins l'un des inconvénients précités.

RESUME DE L'INVENTION

A cet effet l'invention prévoit un procédé de navigation passive pour estimer une déformation d'un cerveau d'un patient, le procédé comprenant : - le traitement 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 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, à 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 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 du cerveau du patient . - la génération, à 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é. L'arbre vasculaire initial est recalé élastiquement vers sa configuration déformée, segmentée dans les images US Doppler localisées, et le champ de déplacements résultant est ensuite étendu à l'ensemble du volume de l'organe par le biais d'un modèle biomécanique spécifique au patient. La déformation globale ainsi calculée permet de mettre à jour les données préopératoires et de rectifier le planning chirurgical initial.

L'invention fournit une réponse rapide, robuste et précise au chirurgien et lui donne la possibilité de valider la pertinence de la déformation calculée avant toute modification du planning. Des aspects préférés mais non limitatifs de la présente invention sont les suivants : - le procédé peut également comprendre le 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 acquise pendant l'opération ; - l'image tridimensionnelle acquise avant l'opération peut être une angiographie par résonance magnétique, - les images bidimensionnelles acquises pendant l'opération peuvent être des images ultrasonores en mode Doppler, - 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 après l'installation de celui-ci sur une table d'opération, - l'étape d'orientation de l'image tridimensionnelle comprend : o la localisation d'un jeu de points de calibrage sur la tête du patient après l'installation de celui-ci sur la table d'opération, o la détection de leurs homologues radiométriques dans l'image tridimensionnelle, et o la détermination de transformations rigides permettant de passer des positions des points du jeu de points de calibrage aux positions de leurs homologues dans l'image tridimensionnelle, - les points de calibrage peuvent comprendre des points de la surface de la tête du patient après l'installation de celui-ci sur la table d'opération, - le traitement des images bidimensionnelles pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient comprend une segmentation des images bidimensionnelles pour obtenir la position médiane des vaisseaux et/ou artères du cerveau de sorte à reconstituer au moins partiellement l'arborescence artérielle cérébrale courante du patient, - l'étape de détermination d'un champ de déformation comprend la détermination de transformations non-linéaires, dites transformations élastiques, permettant le passage des positions des points de l'arborescence artérielle cérébrale de référence aux positions des points correspondants de l'arborescence artérielle courante. L'invention concerne également un système pour la détermination d'une déformation du cerveau d'un patient, le système comprenant des moyens pour la mise en oeuvre du procédé décrit ci-dessus.35 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, à 5 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, 10 - 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, 15 - 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. 20 DESCRIPTION DETAILLEE DE L'INVENTION

On va maintenant décrire plus en détail un mode de réalisation du procédé selon l'invention. 25

Principe général de l'invention Comme expliqué précédemment, la localisation précise de la cible chirurgicale est 30 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. L'invention propose une technique permettant 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, l'invention 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 15 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 20 par la déformation du cerveau. Cette information, loin d'être dense dans le volume global, reste localisée à certaines régions. L'invention 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 25 vasculaire cérébral. Le champ de déplacements global ainsi obtenu est ensuite répercuté à l'ensemble des données préopératoires afin de les mettre à jour.

Description d'un mode de réalisation du procédé selon l'invention 30 En référence à la figure 1, on a illustré différentes étapes du procédé selon l'invention. Le procédé de traitement d'image pour estimer un déplacement d'un cerveau d'un patient comprend les étapes suivantes : - le traitement 100 d'une image tridimensionnelle du cerveau du patient, acquise 35 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é.

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. 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. Préférentiellement encore, cette image tridimensionnelle est une angiographie par résonnance magnétique. L'angiographie par résonance magnétique permet de révéler des contrastes 30 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 35 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 15 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 20 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 25 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 30 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 35 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 15 é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 20 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 25 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 30 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 35 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 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 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 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 20 sont de natures distinctes. 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. 25 D'un autre coté, l'arbre vasculaire courant obtenu à partir des images 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 30 que d'autres ne comportent aucune acquisition. La tâche qui incombe à l'étape 300 de recalage des arbres vasculaires cérébraux 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 35 - la configuration de référence définie par les données complètes et 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é.

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 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 (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. 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 à 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 5 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 10 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 présenté ci-dessus peut s'appliquer au recalage de tout type de tissu organique. 15 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 vasculaire déterminé est appliqué à un modèle biomécanique des tissus mous du cerveau 20 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 limites. 25 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 déformation du cerveau. 30 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 vasculaire courant et l'arbre vasculaire de référence. Le modèle biomécanique permet 35 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). 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 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 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 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 20 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 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 25 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 30 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 35 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 application L 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 recalaqe 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. r 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 : fs) ti 5 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 : 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. ("XI 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. 4.,,f:7ti 1 (.IX i t' 3 = 1,(P)I (1 itX2 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 r )` atteint P 11 t .fl''" ' 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 Afin de simplifier ces considérations, on considère les fonctions R surjectives, i.e. telles que llrr_',rr.'.l,` 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 1F, , 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 : <~. frontière (j 1"1+Cs1E1r1 E' it( %1. ,_C '1 est IIt"Lat.,.çé' l.Ei. l F a - (n), 1,E' )4E1 ï11a e ilt'! Il rut Yll0h< le#1,E;t#{ MI bail, 4 Ta' ;I1 3 = (i11 1_Efi#ti h'-4 pOff1l', à [`0x1{'1'li`11C' 4h'' f' 1'i#14li#Clfaç' r(st'(?Ul' #11f'ÎÏ IT1 1'4 L. T'f(ll) `{.p) = p.

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 20 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 i)ç 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 25 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. 30 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. 35 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 : 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';. Si chaque r'; est une ,)3 bijection alors la fonction de recalage R est également une bijection et son 10 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 + ... + .1 Pi k 15 IJ , alors` est l'ensemble de toutes les positions successives de N- p0 et R(p) = pK. La différentielle de R en p0 est donnée par la règle de composition suivante : Finalement, si tous les r'; vérifient la condition de non-repliement, alors cela est vrai 20 également de l'application R. Ceci est la conséquence de l'expression du jacobien de R : 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 `f . [tt' 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 le 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 p + Les conditions sur la zone d'influence de r se traduisent par des contraintes sur le support de w de la façon suivante : { _p .?u p qui est un compact Remarquons que le support de w est Différentiabilité La différentiabilité de r découle de celle de w. alors ' 3 ) et la matrice jacobienne de r s'écrit : 1 ùP. rhr. i)r où Id3x3 est la matrice identité 3 x 3. -?1#~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 : ) Puisque P ti Iest une application continue à support compact, soit le point où I~ ., "`1I atteint son maximum II:ù ' Vii%'{).rir;::F O. Alors l'équation précédente est équivalente à : Le non-repliement de r peut donc être garanti en limitant l'amplitude des 15 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 20 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 25 Choix d'une expression de w Afin de simplifier le calcul du recalage élastique, une expression polynomiale de w est choisie. Soit 7 un polynôme de degré 3, défini par : Son expression est 2t Une fonction de forme w définie sur la cellule canonique et prenant la valeur 1 au noeud ii =' . peut ainsi être obtenue comme : ii.,tt k = ,u'i ,~1-:~ , P2 P3 23 La figure 6a montre la fonction de forme restreinte à une dimension, i.e. '1' ù 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 [0' Il 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 u' ; R2 k.

Contrainte de non-repliement et de bijection L'expression de w étant choisie, il est possible de calculer l'amplitude du 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é IIlsur le domaine [O. 1~.e La dérivée de {fl~1 ù . Elle atteint son maximum sur [0, 1] en tmax = 1/2 où `. . On obtient l'expression de ar? f r r on a supérieure est valide pour chaque coordonnée, ainsi Il J.. 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 . Cette borne et donc canoniques [0, 1] est donc : Il ' Ce résultat peut être étendu à n'importe quelle dimension de cellule de la façon s s suivante. Soit une cellule cubique d'arête L. Soit E.- la fonction qui transforme cette cellule en cellule canonique [tt'~ , définie par : {,t?J _ 3,? i , .11.1#1 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 : La fonction de forme w définie sur Oeil est ,3.i',f ) ;(] et: Il 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 : (..3$1 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 5 l'espace appliqués au solide à chaque itération sont bornés. Puisque r 10 et , alors : 20 - 20 L'équation ci-dessus fournit une approximation 20 ù 0.29`l) ) des bornes inférieure et supérieure du jacobien de r valable pour toutes les dimensions de cellules : 10 0.7, L26 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 15 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 . 20 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 25 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 30 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 : 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 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 *II sont impliquées dans le calcul de la minimisation. La fonction minimisée est définie par :

où P P + fl''UU) 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 de laisser de coté l'historique du recalage que constitue R. il' l'ensemble de points source dont la position appartient au Soit i,.El 30 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 d'exprimer la fonction à minimiser f comme : iEJ 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 : 5 où f-, 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ù c- ''ti` `? ] . Les vecteurs t-{ k 5.k=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 : 10 Le meilleur déplacement nodal au sens de la minimisation de l'énergie de recalage correspond donc au minimum de la fonction f' () + t pour t [(,?. ]. Une estimation convenable de ce minimum est obtenue en approximant sur l'intervalle [0, 1 ] la fonction 3` par une parabole 1> 4/ + 1,1 + -.' Les trois coefficients de la parabole A, B et C sont déterminés par les informations 15 locales sur le comportement de f ainsi qu'une évaluation de sa valeur au pas maximal : 'i1+(-. +.4 Ir "II t i) . 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 20 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 25 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 ? -'va 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 t.f du noeud n. Pour un point donné t1 ]` n k 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 ù q + f 7 où tp est la racine de la fonction polynomiale f sur l'intervalle [ta, tb]. En fait, puisque w l'intervalle de recherche peut être restreint à [max{-1, ta}, 0]. Le choix de l'expression de la fonction de forme amène à I expression suivante de f : 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 30 avec l'estimation de la solution suivante : Les itérations sont ensuite menées de manière classique : L'approximation de RI) à chaque itération est donnée par l'erreur correspondante dans l'espace déformé est 'i 1111 La précision de l'approximation de l'inverse est mesurée dans l'espace initial, non déformé, en évaluant la quantité 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 inverse r-' est k-lipschitzienne, en d'autres termes : lEÏ "'7'p, q ù' ; .-1 ù (q 10 ù Cette propriété permet alors de majorer l'erreur par - .- .l hall ù(111= i;'I ~I . On désire poursuivre les itérations de l'algorithme de Newton jusqu'à ce que la 15 précision voulue soit atteinte i.e. lorsque I ,, {.. Au vu de ce qui précède, Ysf'; cette condition d'arrêt est vérifiée lorsque ' 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 plusieurs fonctions élémentaires rj . Dans ce cas, la précision du calcul de l'inverse

20 ' 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 constantes de Lipschitz des fonctions

Pour plus de clarté, on suppose que _ , On veut évaluer R _ en ayant, pour chacune des , une 25 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 1'~`' ainsi que l'accumulation successive des erreurs commises à chaque pas : 1. q2 approxime '7 ù '' L'erreur dans l'espace déformé est _3 := kq3-p3k. Cette erreur se propage vers l'espace non-déformé comme r et finalement II( ù 1)11 2. r1 approxime 2`avec une erreur dans l'espace déformé := llr2-q2ll. Cette erreur se propage vers la gauche, i.e. vers l'espace non-déformé, comme puis Ir q 3. Finalement s approxime r _ 11.'1 avec une erreur dans l'espace déformé e -r1. . Cette erreur se propage comme s ù r L'erreur finale peut ainsi être majorée par la somme des erreurs propagées : f _ Ils ù Pll Ils ù III+IIr ùqII+IIq ù II L'expression ci-dessus peut être immédiatement étendue à toute fonction de 15 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 E 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 Newton à: 20 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 25 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é 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 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 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 : F f` (:,.uf-T 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 25 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 30 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 10 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 15 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 estimer le déplacement du 20 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 25 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 de natures distinctes du fait de la différence entre les zones 30 de recouvrement des données ou de la présence de bruit, l'étape de recalage peut être mise en oeuvre sur des jeux de données S et D homogènes. 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 (10)

1. REVENDICATIONS1. Procédé de traitement d'image pour estimer une déformation d'un cerveau d'un patient, le procédé comprenant : - 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 du cerveau du patient . - 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é.
2. 2. Procédé selon la revendication 1, caractérisé en ce qu'il comprend en outre le calcul (600) des disparités entre des points de l'image du cerveau générée et une ou plusieurs images bidimensionnelles de contrôle acquise pendant l'opération.
3. 3. Procédé selon la revendication 1 ou la revendication 2, caractérisé en ce que l'image tridimensionnelle acquise avant l'opération est une angiographie par résonance magnétique. 30
4. 4. Procédé selon l'une des revendications 1 à 3, caractérisé en ce que les images bidimensionnelles acquises pendant l'opération sont des images ultrasonores en mode Doppler.
5. 5. Procédé selon l'une des revendications 1 à 4, caractérisé en ce que le traitement de 35 l'image tridimensionnelle pour obtenir l'arborescence artérielle cérébrale de référence25comprend une étape d'orientation de l'image tridimensionnelle par rapport à la tête du patient après l'installation de celui-ci sur une table d'opération.
6. 6. Procédé selon la revendication 5, caractérisé en ce que l'étape d'orientation de l'image tridimensionnelle comprend : - la localisation d'un jeu de points de calibrage sur la tête du patient après l'installation de celui-ci sur la table d'opération, - la détection de leurs homologues radiométriques dans l'image tridimensionnelle, et - la détermination de transformations rigides permettant de passer des positions des points du jeu de points de calibrage aux positions de leurs homologues dans l'image tridimensionnelle.
7. 7. Procédé selon la revendication 6, caractérisé en ce que les points de calibrage comprennent des points de la surface de la tête du patient après l'installation de celui-ci sur la table d'opération.
8. 8. Procédé selon l'une des revendications 1 à 7, caractérisé en ce que le traitement des images bidimensionnelles pour reconstituer au moins partiellement une arborescence artérielle cérébrale courante du patient comprend une segmentation des images bidimensionnelles pour obtenir la position médiane des vaisseaux et/ou artères du cerveau de sorte à reconstituer au moins partiellement l'arborescence artérielle cérébrale courante du patient.
9. 9. Procédé selon l'une des revendications 1 à 8, caractérisé en ce que l'étape de détermination d'un champ de déformation comprend la détermination de transformations non-linéaires, dites transformations élastiques, permettant le passage des positions des points de l'arborescence artérielle cérébrale de référence aux positions des points correspondants de l'arborescence artérielle courante.
10. 10. Système pour la détermination d'une déformation du cerveau d'un patient, caractérisé en ce qu'il comprend des moyens pour la mise en oeuvre du procédé selon l'une des revendications 1 à 9.