Brevet FR2979453A1 — publié le 20130301

Système et procédé pour estimer une quantité d'intérêt d'un système dynamique artère/tissu/veine L'invention concerne un système et un procédé pour estimer des paramètres hémodynamiques à partir d'images médicales. L'invention se distingue notamment des procédés connus par sa précision et selon le mode de réalisation par une rapidité d'exécution indispensable pour permettre un diagnostic thérapeutique en situation d'urgence. L'invention s'appuie notamment sur des techniques 10 d'imagerie de Perfusion par Résonance Magnétique (Perfusion Weighted Magnetic Resonance Imaging - PW-MRI selon une terminologie anglo-saxonne) ou par Tomodensitométrie (Computed Tomography - CT - selon une terminologie anglo-saxonne). Ces techniques permettent 15 d'obtenir des informations précieuses sur l'hémodynamique d'organes tels que le cerveau ou le coeur. Ces informations sont particulièrement cruciales pour un praticien cherchant à établir un diagnostic et à prendre une décision thérapeutique dans le traitement 20 en hyper-urgence de pathologies telles que les accidents vasculaires cérébraux. Pour mettre en oeuvre de telles techniques, on utilise un système d'imagerie médicale tel qu'illustré 25 à titre d'exemple par les figures 1 et 2 exploitant un appareil 1 d'imagerie par Résonance Magnétique Nucléaire ou par Tomodensitométrie. Celui-ci délivre une pluralité de séquences d'images numériques 12 d'une partie du corps, notamment du cerveau. Ledit appareil 30 applique pour cela une combinaison d'ondes électromagnétiques à haute fréquence sur la partie du corps considérée et mesure le signal réémis par certains atomes. L'appareil permet ainsi de déterminer la composition chimique et donc la nature des tissus biologiques en chaque point (ou voxel) du volume imagé. Des séquences d'images sont analysées au moyen d'une unité de traitement dédiée 4. Cette unité de traitement délivre in fine à un praticien 6, une estimation des paramètres hémodynamiques à partir des images de perfusion, au moyen d'une interface homme- machine 5 adaptée. Le praticien peut ainsi réaliser un diagnostic et décider de l'action thérapeutique qu'il jugera adéquate. Des images de perfusion par Résonance Magnétique Nucléaire ou par Tomodensitométrie sont obtenues en injectant un agent de contraste (par exemple un sel de gadolinium pour l'Imagerie par Résonance Magnétique) par voie intraveineuse et en enregistrant son bol au cours du temps au niveau de chaque voxel de l'image. Par souci de concision, nous omettrons les indices x,y,z pour identifier des voxels. Par exemple, au lieu de noter Sx,Y,z (t) le signal pour un voxel de coordonnées x,y,z, nous le noterons simplement S(t). Il est entendu que les opérations et les calculs décrits dans la suite 25 sont généralement effectués pour chaque voxel d'intérêt, de sorte à obtenir au final des images ou des cartes représentatives des paramètres hémodynamiques que l'on cherche à estimer. 30 Un modèle standard permet de relier l'intensité des signaux S(t) mesurée au cours du temps t à la concentration C(t) dudit agent de contraste. Par exemple, en Tomodensitométrie de Perfusion, le 35 signal pour chaque voxel est supposé directement proportionnel à la concentration : S(t)= k.C(t)+3. En Imagerie de Perfusion par Résonance Magnétique Nucléaire, on suppose par exemple qu'il existe une -k.TE.C(t) relation exponentielle S(t)=So-e . Dans les deux cas, SD représente l'intensité moyenne du signal avant l'arrivée de l'agent de contraste. Dans le cas de l'imagerie par Résonance magnétique nucléaire, k est une constante dépendant de la relation entre la susceptibilité paramagnétique et la concentration de l'agent de contraste dans le tissu et TE est le temps d'écho (echo time selon une terminologie anglo-saxonne). La valeur de la constante k pour chaque voxel étant inconnue, celle-ci est fixée à une valeur arbitraire pour tous les voxels d'intérêt. On obtient ainsi des estimations relatives et non pas absolues. Ces informations relatives restent toutefois pertinentes puisque l'on est intéressé principalement par la variation relative de ces valeurs dans l'espace, en particulier entre les tissus sains et les tissus pathologiques. Dans toute la suite nous supposerons les signaux d'intensité expérimentaux S(t) préalablement convertis en courbes de concentration CM. Par exemple, dans le 25 cas de l'imagerie de perfusion par résonance magnétique, nous avons C(t)=-1n , étant estimée en prenant, par exemple, la moyenne de SM avant l'arrivée de l'agent de contraste. 30 La conservation de la masse de l'agent de contraste dans le volume de tissu contenu dans chaque voxel à ( chaque instant s'écrit dC t) BF -[Ca(t)-Cv(t)1 . Ca(t) est dt la concentration de l'agent de contraste dans l'artère alimentant le volume de tissu (fonction d'entrée artérielle ou Arterial Input Fonction - AIF - selon une terminologie anglo-saxonne). BF est le flux sanguin dans le volume de tissu (Blood Flow selon une terminologie anglo-saxonne) et G(t) est la concentration de l'agent de contraste dans la veine drainant le volume de tissu (fonction de sortie veineuse ou Venous Output Function - VOF - selon une 10 terminologie anglo-saxonne). En supposant le système dynamique artère/tissu/veine linéaire et invariant dans le temps, on peut écrire Cv(t)=Ca(t)®h(t) où h(t) est la réponse 15 impulsionnelle du système - ou encore fonction densité de probabilité du temps de transit de l'agent de contraste dans le tissu - et ® désigne le produit de convolution. Une solution formelle de l'équation différentielle précédente avec condition initiale 20 C(t=0)=0 s'écrit alors C(t)=BF-Ca(t)OR(t) où R(t) est la fonction de répartition complémentaire du temps de transit dans le volume de tissu (residue fonction selon une terminologie anglo-saxonne) définie par R(t)=H(t)-fh(t)d-t- où H est la fonction généralisée 0 25 créneau de Heaviside. A partir de la réponse impulsionnelle et de la fonction de répartition complémentaire, on définit un nouveau paramètre hémodynamique, le temps de transit moyen dans le tissu (Mean Transit Time - MTT - selon une terminologie 30 anglo-saxonne) : +c +cc MTT, f r- h (t)dt = f R(t)dt (si limt - h(t) = 0 ) -> On peut définir également le volume sanguin dans le volume de tissu (Blood Volume ou BV selon une terminologie anglo-saxonne) par la relation BV=BF-MTT. Afin d'estimer les paramètres hémodynamiques tels que BF, mu ou BV ainsi que la fonction de répartition complémentaire R(t), il est donc nécessaire de déconvoluer les courbes de concentration C(t) par 10 des fonctions d'entrée artérielles Ca(t) que l'on supposera données par la suite. Pour effectuer cette opération de déconvolution de C(t) par Ca(t), le modèle standard de convolution 15 CM= BF - Ca(t)0 R(t) est tout d'abord discrétisé temporellement aux instants d'échantillonnage du signal en approchant numériquement l'intégrale de convolution CeORM. Sans perte de généralité, nous supposerons dans la suite l'échantillonnage périodique, 20 de période At=t1-ti_1. Par exemple, l'approximation de l'intégrale de convolution par la méthode des rectangles donne 25 Vi=1,N, C(ti)= BF - j Ca(T)- t -T)d-r BF - At - I -tk) 0 k=0 On se ramène ainsi à un système linéaire Ad =c de dimension N en posant ( Ça (t)o o Ça(t2) ca(ti) *.. . : .. o Ca(tN) ca(tN_i) ... ca(ti) ) A= 0t. , b= R(q) c(ti) c(t2) R(t2) c(tN) , c= R(tN) et d= BF. b. La matrice de Toeplitz triangulaire basse A est très mal conditionnée et presque singulière, de sorte qu'on ne peut pas inverser numériquement ce système linéaire sous peine d'obtenir des solutions dénuées de sens et des estimations des paramètres hémodynamiques aberrantes. Il faut donc recourir à diverses méthodes pour obtenir par exemple une pseudo-inverse À 1 de la ....... matrice A et par suite une estimation d de d par d= À-1.c . Parmi les nombreuses méthodes d'obtention d'une pseudo-inverse pour la matrice A, on trouve les méthodes classiques non paramétriques basées sur la troncature des valeurs singulières de A telles que les méthodes sSVD (Simple Singular Value Decomposition, - selon une terminologie anglo-saxonne), cSVD (Circular Singular Value Decomposition) et oSVD (Oscillation Index Singular Value Decomposition). La méthode sSVD a l'avantage d'être simple et 25 rapide. Néanmoins, elle souffre de deux inconvénients majeurs : - Elle est sensible au délai temporel T entre la fonction d'entrée artérielle Ca(t) et la courbe de concentration C(t), c'est-à-dire 30 qu'elle fournit des estimations des paramètres tels que BF et Aur qui dépendent de T alors qu'elles ne le devraient pas ; - En particulier, elle fournit des estimations des paramètres aberrantes lorsque ledit délai T est négatif, c'est-à-dire lorsque la fonction d'entrée artérielle Ca(t) est en retard par rapport à la courbe de concentration C(t). Ces inconvénients sont palliés par les méthodes cSVD et oSVD qui sont par construction insensibles au délai T et permettent de prendre en compte les délais négatifs. La méthode cSVD (ainsi que la méthode sSVD) souffre néanmoins d'un inconvénient important : elle n'est pas adaptative : la régularisation est déterminée une fois pour toute à l'aide d'un paramètre PSVD propre à l'algorithme qu'on peut interpréter comme la fréquence de coupure d'un filtre passe-bas. Or la régularisation devrait au contraire s'adapter à chaque courbe de concentration expérimentale C(t), en particulier à leur rapport signal-sur-bruit. En pratique, on devrait donc préalablement déterminer une valeur du paramètre PSVD adaptée à chaque jeu de données de perfusion d'intérêt en déterminant par exemple la valeur permettant d'optimiser un critère donné (e.g. l'erreur relative sur un paramètre, etc.). Or cela est stricto sensu impossible puisque les valeurs théoriques des paramètres ne sont pas connues. D'autre part, une telle calibration pour chaque modalité de perfusion (e.g. CT Perfusion ou MR PWI), chaque appareil de mesure, chaque jeu des paramètres d'acquisition et même chaque type de signaux de perfusion (e.g. perfusion dans la matière blanche, dans la matière grise, perfusion saine ou pathologique) n'est clairement pas souhaitable. En pratique, ces calibrations sont rarement effectuées et le paramètre PSVD est souvent fixé de manière assez arbitraire. La méthode oSVD permet de pallier en partie cet inconvénient en introduisant une régularisation semiadaptative la rendant moins sensible aux différentes conditions expérimentales et aux différents types de signaux de perfusion. Cette méthode est par construction invariante par délai et semi-adaptative. En conséquence, la méthode oSVD est sensée fournir des estimations des paramètres hémodynamiques de meilleure qualité et plus robustes vis-à-vis des conditions opératoires que ne le font les méthodes cSVD et sSVD. En revanche, la mise en oeuvre d'une telle méthode est complexe et non suffisamment documentée pour rendre celle-ci applicable efficacement en milieu hospitalier. Il en est de même pour la méthode cSVD. oSVD et cSVD restent à ce jour théoriques. En outre, comme nous le verrons plus loin, la complexité algorithmique de la oSVD peut être de l'ordre de N fois celle de la cSVD et de 4N fois celle de la sSVD. Si les temps de calcul pour des jeux de données typiques d'imagerie de perfusion avec des moyens de calcul matériels conventionnels (i.e. ordinateur personnel ou station de travail) sont de l'ordre de quelques secondes pour les méthodes sSVD et cSVD, ils peuvent donc atteindre plusieurs minutes avec la méthode oSVD. Cela est parfaitement rédhibitoire en situation d'hyper-urgence clinique, comme par exemple lors du traitement des accidents vasculaires cérébraux où l'on estime que 4 millions de neurones meurent chaque minute. La présente invention consiste selon un premier objet à implémenter la méthode oSVD pour pouvoir l'exploiter dans un système d'imagerie médicale. L'invention porte également sur un mode de réalisation particulièrement optimisé du point de vue du faible temps nécessaire à l'estimation d'une quantité d'intérêt. L'invention est illustrée préférentiellement - mais de manière non limitative en lien avec l'exploitation d'un système d'imagerie médicale pour estimer une ou plusieurs quantités d'intérêts (ou paramètres(s) hémodynamique(s) par imagerie de perfusion. Le mode d'implémentation optimisée procure une complexité algorithmique du même ordre que celle de la méthode cSVD. On obtient dès lors des temps de calcul pour la oSVD conforme à l'invention de l'ordre de quelques secondes avec des moyens de calcul matériels conventionnels, acceptables en situation d'hyper-urgence clinique. Selon un second objet, l'invention permet d'appliquer une approche équivalente pour favoriser la mise en en oeuvre de la méthode CSVD. A cette fin, il est prévu un premier mode de réalisation d'un procédé pour produire une estimation d d'une quantité d'intérêt d d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe. Un tel procédé est mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie médicale, ladite estimation consistant à calculer êI=À 1.C, À-1 étant une pseudo-inverse d'une matrice de convolution A et C décrivant une concentration d'un agent de contraste dans ledit voxel. Selon ce premier mode de réalisation, le procédé comporte : - une étape pour décomposer canoniquement A en valeurs singulières sous la forme A=U-S.VT où SC/iag(csi,...,GL) est la matrice diagonale des valeurs singulières classées par ordre croissant, VT est la transposée d'une matrice V, V=(Ifij et U=W sont deux matrices carrées de dimension LxL unitaires et réelles, LN, N étant un nombre d'échantillons déterminé, d, A et C étant respectivement de dimensions Lxl, LxL et Lxl ; - au moins une étape pour : 1 o élaborer une pseudo-inverse /V de A, 1 sous la forme /V = V- Vli - UT, UT désignant la transposée de U avec ( Wj =diag 0,...0,1/c1+1,...,1/c y L 0 Afin de répondre aux situations d'urgences cliniques, l'invention prévoit un second mode de réalisation optimisée du point de vue du temps de calcul pour produire une estimation d d'une quantité d'intérêt d d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe. Au même titre que le procédé précédent, celui-ci est mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie médicale, ladite estimation consistant à calculer (9= C , À-1 étant une pseudo- inverse d'une matrice de convolution A et C décrivant une concentration d'un agent de contraste dans ledit voxel. Selon l'invention, le procédé comporte toujours une première étape pour décomposer canoniquement A en valeurs singulières sous la forme A=U-S.VT où S= L) est la matrice diagonale des valeurs singulières positives classées par ordre croissant, VT et Ll=4) sont deux matrices carrées de dimension LxL unitaires et réelles, L>iV, N étant un nombre d'échantillons déterminé, d, A et C étant respectivement de dimensions Lxl, LxL et Lxl. En revanche, il comporte une étape itérative pour, à l'itération 0 Pour initialiser l'itération éventuelle de l'étape pour produire dj, un tel procédé peut comporter une étape initiale pour produire d0=A10.C. A titre d'exemple préféré, la condition de 5 régularisation peut consister en le calcul d'un index d'oscillation 04, ladite condition de régularisation étant satisfaite dès que 011F1301, Fa étant un seul prédéterminé. Un tel index d'oscillation peut être L calculé tel que 04 L max1 i = d! (i) -2di (i -1)+ di (i -2)1 . di ( . i=1 )L 1=3 10 Quel que soit le mode de réalisation conforme à l'invention choisi, la matrice de convolution A peut être bloc-circulante de dimension LxL définie par ( Ca(t0 0 - Ca(tN) - Ca(t2) - Ca (t0 . . A= At - Ca (t N) Ca(tN) At étant la 0 0 0 0 Ca(tN) --- Ca(ti) ) d'échantillonnage, Ca (t) étant une période 15 concentration de l'agent de contraste. Par ailleurs, pour pourvoir in fine délivrer au praticien un ou plusieurs paramètres hémodynamiques, un procédé conforme à l'invention peut comporter une étape pour produire une estimation d'un paramètre 20 hémodynamique à partir de l'estimation produite éi=diF . Pour mettre en oeuvre un procédé conforme à l'invention, celle-ci prévoit également une unité de traitement comportant des moyens de mémorisation, des moyens pour communiquer avec le monde extérieur et des 25 moyens de traitement. Les moyens pour communiquer sont aptes à recevoir du monde extérieur une donnée de expérimentale C et les moyens de traitement sont adaptés pour mettre en oeuvre un procédé pour estimer une quantité d'intérêt d selon l'invention. Préférentiellement, les moyens pour communiquer d'une telle unité de traitement peuvent délivrer une quantité d'intérêt estimée d selon un format approprié à une interface homme-machine apte à la restituer à un utilisateur. Si le procédé mis en oeuvre permet conformément à l'invention de produire l'estimation d'un paramètre hémodynamique, une telle unité de traitement peut être adaptée pour délivrer - selon un format approprié ladite estimation à une interface homme-machine apte à la restituer à un utilisateur. L'invention concerne également tout système d'analyse d'imagerie médicale qui comporte une unité de traitement adaptée conformément à l'invention et une interface homme-machine apte à restituer à un utilisateur une quantité estimée selon un procédé conforme à l'invention et mis en oeuvre par ladite unité de traitement. D'autres caractéristiques et avantages apparaîtront plus clairement à la lecture de la description qui suit et à l'examen des figures qui l'accompagnent parmi 25 lesquelles : - les figures 1 et 2, partiellement précédemment décrites, présentent deux variantes de réalisation d'un système d'analyse d'imagerie médicale ; 30 - Les figures 3 et 4 présentent respectivement une image de perfusion, obtenue par un appareil d'imagerie par Résonance Magnétique Nucléaire, d'une tranche d'un cerveau humain avant l'injection d'un agent de contraste et pendant la circulation de celui-ci dans les tissus dudit cerveau ; - les figures 5a et 5b présentent un signal de perfusion S(t) par Résonance Magnétique Nucléaire typique relatif à un voxel d'un cerveau humain ; - la figure 6 présente une courbe de concentration C(t) typique d'un agent de contraste circulant au sein d'un voxel d'un cerveau humain; - la figure 7 présente une fonction d'entrée artérielle Ca(t) typique ; - les figures 8 et 9 présentent des procédés conformes à l'invention pour mettre en oeuvre la méthode oSVD ; - la figure 8b présente un procédé conforme à l'invention pour implémenter la méthode cSVD ; - les figures 10 et 11 présentent respectivement un exemple de carte relative à des quantités d'intérêt estimées conformément à l'invention. De manière plus détaillée, la figure 1 permet de présenter un système d'analyse d'images médicales. Un appareil 1 d'imagerie par Résonance Magnétique Nucléaire ou par Tomodensitométrie est commandé à l'aide d'une console 2. Un utilisateur peut ainsi choisir des paramètres 11 pour piloter l'appareil 1. A partir d'informations 10 produites par l'appareil 1, on obtient une pluralité de séquences d'images numériques 12 d'une partie d'un corps d'un humain ou d'un animal. A titre d'exemple préféré, nous illustrerons les solutions issues de l'art antérieur ainsi que l'invention à l'aide d'images numériques issues de l'observation d'un cerveau humain. D'autres organes pourraient aussi être considérés. Les séquences d'images 12 peuvent optionnellement être stockées au sein d'un serveur 3 et constituer un dossier médical 13 d'un patient. Un tel dossier 13 peut comprendre des images de différents types, telles que des images de perfusion ou de diffusion. Les séquences d'images 12 sont analysées au moyen d'une unité de traitement 4 dédiée. Ladite unité de traitement comporte des moyens pour communiquer avec le monde extérieur pour recueillir les images. Lesdits moyens pour communiquer permettent en outre à l'unité de traitement de délivrer in fine à un praticien 6 ou à un chercheur, une estimation de paramètres hémodynamiques 14 à partir des images de perfusion 12, au moyen d'une interface homme-machine adaptée 5. L'utilisateur 6 du système d'analyse peut ainsi confirmer ou infirmer un diagnostic, décider une action thérapeutique qu'il jugera adéquate, approfondir des travaux de recherche... Optionnellement, cet utilisateur peut paramétrer le fonctionnement de l'unité de traitement 4 au moyen de paramètres 16. Par exemple, il peut ainsi définir des seuils d'affichage ou choisir les paramètres estimés qu'il souhaite visualiser. La figure 2 illustre une variante de réalisation d'un système d'analyse pour lequel une unité de prétraitement 7 analyse des séquences d'images 12 pour en déduire des données de perfusion 15 par voxel. L'unité de traitement 4 chargée d'estimer les paramètres hémodynamiques 14 est ainsi déchargée de cette action et met en oeuvre un procédé d'estimation à partir de données de perfusion 15 réceptionnés par ses moyens pour communiquer avec le monde extérieur. La figure 3 illustre un exemple d'image typique 12 35 d'une tranche de 5 millimètres d'épaisseur d'un cerveau humain. Cette image est obtenue par Résonance Magnétique Nucléaire. A l'aide de cette technique, on peut obtenir, pour chaque tranche, une matrice de 128 x 128 voxels dont les dimensions sont de 1,5 x 1,5 x 5 millimètres. A l'aide d'une interpolation bilinéaire on peut produire une image à plat de 458 x 458 pixels telle que l'image 20. La figure 4 illustre une image 20 similaire à celle présentée en liaison avec la figure 3. Toutefois cette image est obtenue après une injection d'un agent de contraste. Cette image est un exemple d'image de perfusion typique d'un cerveau. Les artères apparaissent ainsi clairement contrairement à la même image décrite en figure 3. Selon des techniques connues, il est possible de choisir une ou plusieurs fonctions d'entrée artérielle 21 dans l'hémisphère controlatérale à l'hémisphère pathologique pour estimer des paramètres hémodynamiques. La figure 5b permet d'illustrer un exemple de signal de perfusion S(t) par Résonance Magnétique Nucléaire tel que les données 15 délivrés par l'unité de prétraitement 7 décrite en liaison avec la figure 2. Le signal de perfusion est ainsi représentatif de l'évolution d'un voxel au cours du temps t à la suite d'une injection d'un agent de contraste. A titre d'exemple, la figure 5b décrit un tel signal sur une durée de 50 secondes. L'axe des ordonnées décrit l'intensité du signal dont l'unité est arbitraire. Pour obtenir un tel signal, l'unité de traitement 4 selon la figure 1 (ou en variante l'unité de prétraitement 7 selon la figure 2), analyse une séquence de n images de perfusion par Résonance Magnétique Nucléaire I1, 12, ..., Ii, ..., In à des instants ti,t2,...,ti,",tn comme le 35 décrit, à titre d'exemple, la figure 5a. Pour un voxel donné, par exemple pour le voxel V, on détermine un signal de perfusion S(t) représentatif de l'évolution du voxel au cours du temps t à la suite d'une injection d'un agent de contraste. La figure 6 présente une courbe de concentration déduite d'un signal de perfusion tel que celui-ci décrit en figure 5b. Comme déjà évoqué précédemment, il existe une relation entre un signal de perfusion et une courbe de concentration associée. Ainsi, en Imagerie de Perfusion par Résonance Magnétique Nucléaire, il existe -klEr(t) une relation exponentielle S(t)=S0.e où S0 est l'intensité moyenne du signal avant l'arrivée de l'agent de contraste, TE est le temps d'écho (echo time) et k est une constante dépendant de la relation entre la susceptibilité paramagnétique et la concentration de l'agent de contraste dans le tissu. La figure 6 permet ainsi de visualiser au cours du temps, l'évolution de la concentration d'un agent de contraste au sein d'un voxel. On note un pic de forte amplitude lors du premier passage (first pass, en langue anglaise) de l'agent de contraste dans le voxel suivi de pics d'amplitudes plus faibles lié à un phénomène de recirculation (second pass, en langue anglaise) dudit agent de contraste. La figure 7 illustre quant à elle, une fonction d'entrée artérielle typique Ca(t) représentative de la circulation d'un agent de contraste au sein d'un voxel artériel tel que le voxel 21 présenté en liaison avec la figure 4. La figure 7 permet de constater notamment que le phénomène de recirculation après un premier passage de l'agent de contraste est très faible. La figure 8 permet de décrire une implémentation conforme à l'invention de la méthode oSVD. Cette réalisation est appliquée préférentiellement à l'imagerie de perfusion. Un procédé pour estimer une quantité d'intérêt peut comporter une première étape initiale 100 pour sélectionner une fonction d'entrée artérielle Ca(t). Selon la méthode oSVD, Le modèle standard de la perfusion est discrétisé temporellement aux instants de 10 mesure sous la forme d'un système linéaire C=A.d avec: Ca(t0 0 ... Ca(tN) Ca(t2) R(ti) Ca(ti) A= At - Ca(tN) - - Ca(tN) R(t Ai) 0 - 0 0 . . . . 0 0 Ca(tN) - - - Ca(ti) C(ti) C(tN) et d= 8F- b. 0 0 A, b, C et d sont de dimensions respectives LxL, Lxl, Lxl et Lxl avec L>IV. Dans la suite, on note 15 14i), i=1,L, les composantes d'un vecteur V de dimension L. Une étape 101 consiste ainsi à construire la matrice de convolution circulaire bloc-circulante A. Selon l'invention, la matrice de convolution 20 circulaire bloc-circulante A est tout d'abord - en 102 - décomposée canoniquement en valeurs singulières sous la forme A= U - S. VT où Sdiag(csi,...,G L) est la matrice diagonale des valeurs singulières positives classées par exemple par ordre croissant (i.e. 0 C= U=S) et V = (vu) sont deux matrices carrées L x L unitaires réelles et VT désigne ici la transposée de V Pour tout voxel Vi d'intérêt, un procédé conforme à l'invention consiste à mettre en oeuvre une étape itérative 130 pour produire une estimation d d'une quantité d'intérêt d. Cette production est préférentiellement réalisée tant qu'une condition de régularisation 112 n'est pas vérifiée. Pour cela, on calcule au préalable en 111 et à chaque itération I, un critère tel qu'un index d'oscillation 04 (oscillation index OI, suivant une terminologie anglo-saxonne). A titre d'exemple préféré, ladite condition de régularisation est satisfaite dès que 04 est supérieur ou égal à une valeur donnée Fa prédéfinie strictement positive. L'itération de l'étape 130 s'arrête dès que ladite condition est satisfaite à l'itération 4=. Un tel procédé considère dès lors, en 114, que l'estimation d est OF. Selon l'invention, il devient possible de produire en 115 par exemple et respectivement des estimations des paramètres hémodynamiques BF, AfTT, BV, ou encore L L du vecteur b par ffexdiF(0, i=1,L At F (i) , At -IcliF (i) max(i) i=1,L F 1=1 i=1 diF et . Pour mettre en oeuvre l'itération ou les max(h i=1,L F itérations nécessaire(s) à la production de la quantité d'intérêt, l'invention prévoit une étape 110 pour initialiser tout paramètre nécessaire à ladite production. Selon un premier mode de réalisation conforme à l'invention et décrit en liaison avec les figures 8 et 9, la production 130 de la quantité d'intérêt d consiste en deux sous-étapes 131 et 132. Ainsi on -1 produit 131 à l'itération I une pseudo-inverse /V de -1 A sous la forme Aj =V-141-UT où la matrice diagonale 01diag(VVI,...,WL) est elle-même obtenue en procédant par itérations successives de la manière suivante. A la première itération (/=0) : = S-1= . Lors d'une itération suivante et éventuelle (1=1), 10 la matrice vy est obtenue en annulant le premier coefficient diagonal w1=1/61 de la matrice W0 tel que vg = diag( 0, 1 / 2 , 1 / et ainsi de suite en annulant à chaque fois le coefficient diagonal 14/1 de la matrice Mi_A jusqu'à ce que la condition de régularisation 15 OliF1301 soit remplie à l'itération /=/F. On a donc L'étape 131 permet de produire alors la pseudo-1 inverse Ai =1.1-1/11-UT puis on calcule en 132 d=BF.b par -1 dj= AI - c . 20 Avantageusement l'index d'oscillation 0/1 consiste en : L 01 = 1 1 (i -1) + (i-2) L max dj (i) i=1, L 1=3 Ce premier mode de réalisation, que nous 25 qualifierons de « direct » de la méthode oSVD, consiste donc à calculer le double produit matriciel -1 -1 AI =1/-Vli-UT puis le produit matrice-vecteur d=icv .c à chaque itération I de l'algorithme, jusqu'à ce que la condition de régularisation 01P01 soit remplie. Avantageusement, il peut être prévu en variante de pré- calculer une fois pour toute toutes les pseudo-inverses 1 AI , I=0,L-1 qui ne dépendent que de la fonction 5 d'entrée artérielle et non pas des courbes de concentration. Selon cette variante, l'étape 131 est réalisée préalablement à toutes celles mises en oeuvre pour un voxel Vi d'intérêt, par exemple à l'issue de l'étape 102 où A est décomposée canoniquement en 10 valeurs singulières. Nous pouvons constater, y compris en implémentant cette variante qui permet de négliger l'impact de 1 l'étape de pré-calcul des pseudo-inverses /V par 15 rapport au traitement de l'ensemble des voxels Vi d'intérêt, qu'en considérant que le nombre IF d'itérations nécessaires pour remplir la condition de régularisation OliF1301 est de l'ordre de N, la mise en oeuvre de la méthode oSVD est donc d'une complexité 20 algorithmique a priori de l'ordre de 0(L2xN)=0(4A3) pour chaque voxel d'intérêt, puisque les -1 multiplications matrice-vecteur d/= A/ .c requièrent O(L2) opérations chacune. 25 En particulier, la complexité d'un tel procédé conforme à l'invention est de l'ordre de N fois celle d'un algorithme pour mettre en oeuvre la cSVD et de 4N fois celle d'un algorithme pour mettre en oeuvre la sSVD. Le nombre de mesures N étant typiquement compris 30 entre 40 (e.g. en CT perfusion) et 70 (e.g. en MR PWI), l'algorithme oSVD est considérablement plus lent que les algorithmes cSVD et sSVD. Pourtant, en tant qu'elle est invariante par délai et semi-adaptative, la méthode oSVD fournit des estimations des paramètres hémodynamiques de meilleure qualité et plus robustes vis-à-vis des conditions opératoires que ne le font les méthodes cSVD et sSVD. Elle leur est donc préférable si le temps de calcul n'est pas un critère primordial. En situation d'urgence, un tel temps de mise en 10 oeuvre pour produire une estimation d'une quantité d'intérêt peut être prohibitif. On rappelle que le coût des moyens de calcul matériels conventionnels est typiquement proportionnel 15 à leur puissance: une machine 50 fois plus puissante coûte typiquement 50 fois plus cher, de sorte qu'il n'est pas économiquement viable de s'en remettre à des moyens matériels conventionnels plus puissants afin d'obtenir des temps de calcul acceptables pour la 20 méthode oSVD en situation d'hyper-urgence clinique. Certes, il est possible d'envisager des moyens de calcul matériels dédiés moins coûteux (e.g. cartes graphiques, cartes FPGA, Field Programmable Gate Array suivant une terminologie anglo-saxonne), mais par 25 définition de telles solutions matérielles dédiées peuvent être coûteuses à développer et d'une applicabilité limitée. La présente invention consiste en outre en une 30 implémentation rapide de la méthode oSVD, d'une complexité algorithmique du même ordre que celle de la méthode cSVD. L'invention permet ainsi d'obtenir des temps de calcul pour la oSVD de l'ordre de quelques secondes avec des moyens de calcul matériels 35 conventionnels, acceptables en situation d'hyper- urgence clinique. Un tel procédé conforme à un second mode de réalisation est illustré en liaison avec la figure 8. Un tel procédé reprend les étapes d'un procédé conforme au premier mode de réalisation. L'étape 130 est principalement différente. En 102, la matrice A est toujours décomposée canoniquement en valeurs singulières sous la forme A = U - S. V 7- où est la matrice diagonale des valeurs singulières positives. Celles-ci sont à présent 10 avantageusement classées par ordre croissant (i.e. 61 (TY2 ---cYL Le principe de ce second mode de réalisation consiste à éviter de devoir pré-calculer (étape(s) 131 -1 15 pour le premier mode de réalisation) les matrices A -1 et surtout les produits matrice-vecteur d=,41/ .c (en 132 selon le précédent mode de réalisation) lors de chaque itération de l'algorithme. Voyons comment cela devient possible selon ce deuxième mode de réalisation. 20 A l'itération 1>0 d'un procédé selon le premier mode de réalisation tel que décrit précédemment, nous avons donc Vil = di ag 00',0,14/14,..,WL de sorte que / ; 0 0 1 25 w1+1u1,1+1 - 0 1+1 - - - w1+1 uL,1+1 wLui,L - wLuL,L 2 L i En notant aii les coefficients de la matrice -1 Al1 = V- Vil -UT, nous avons donc par définition L du= 1 vikwkuik . k=I+1 -1 Le i-ème coefficient di(i) du vecteur di =A1 - c s'écrit alors N N L di (i) = 1 alu - ci =Ici - 1 vikwkujk j=1 j =1 k=I+1 puisque Ci = 0 si N On peut alors écrire : N L 0_1(0=1C./ - 1 vi kwkuj k j=1 k=1 N L N =1 Cj - 1 vikwkujk+ viiwiuji = di (i) +1 cj - viiwiuji j=1 U=I+1 ) j=1 de sorte à obtenir la formule de récurrence d'ordre un N di (i)= cli_i (i)- G i Il ci - uji j=1 On constate donc que l'on peut exprimer directement le vecteur di à partir du vecteur d1_1 obtenu à l'itération précédente de l'algorithme oSVD sans avoir à effectuer une quelconque opération matricielle mais N seulement en calculant une seule somme IC.i.t/j/- j=1 Grâce à cette formule de récurrence mise en oeuvre à l'étape 130 d'un procédé conforme à celui décrit en liaison avec la figure 8, on passe d'une complexité 5 algorithmique pour chaque voxel d'intérêt Vi de l'ordre de 0( L2 x /V) suivant l'implémentation directe à une complexité algorithmique de l'ordre de 0(L2+Lx/V) puisqu'il faut initialiser - étape 110 - l'algorithme en calculant do=A-1.0 de complexité 0(L2). 10 Le coût supplémentaire en temps de calcul de cette implémentation par rapport à l'algorithme cSVD n'est donc plus que de l'ordre de 0(L x N) = 0(2N2), à comparer à la complexité algorithmique de la cSVD elle-même en 15 0(L2) = 0(4N2) . En pratique, ce second mode d'implémentation rapide ou optimisée de la méthode oSVD est plus lente que l'algorithme cSVD seulement d'un facteur 2 ou 3 au lieu d'un facteur de l'ordre de N suivant l'implémentation directe conforme au premier 20 mode d'implémentation précédemment décrit. Les calculs peuvent donc être effectués à nouveau en quelques secondes sur matériel conventionnel. Ainsi, l'invention rend possible l'usage de la méthode oSVD 25 pour l'estimation de paramètres hémodynamiques par imagerie de perfusion (par exemple) sur moyens de calcul matériels conventionnels en situation d'hyperurgence clinique. 30 Un procédé conforme à l'invention permet de mettre en oeuvre la méthode cSVD selon un mode opératoire similaire au premier mode de réalisation d'un procédé conforme à l'invention pour mettre en oeuvre la méthode oSVD. La figure 8b illustre un tel procédé appliqué à la méthode cSVD. Ainsi, au même titre qu'un procédé conforme aux figures 8 et 9, un procédé conforme à l'invention pour mettre en oeuvre la méthode cSVD comporte une première étape 100 pour sélectionner une fonction d'entrée artérielle Ce. Selon la méthode cSVD, le modèle standard de la 10 perfusion est également discrétisé temporellement aux et d = BF. b. instants de mesure C= A. d avec: Ca(tO 0 Ca(tO Ca(tN) 0 0 C(ti) C(tN) 0 0 sous la forme d'un système linéaire Ca(tN) Ca(t2) Ca(tN) 0 0 Ca(tN) Ca(ti) ) A= At - C= R(ti) b R(tN) 0 0 15 A, b, C et d sont de dimensions respectives LxL, Lxl, Lxl et Lxl avec L1\1. Une étape 101 consiste ainsi à construire la matrice de convolution circulaire bloc-circulante A. Selon l'invention, la matrice de convolution 20 circulaire bloc-circulante A est tout d'abord - en 102 - décomposée canoniquement en valeurs singulières sous la forme A= U - S. VT où S. di ag(G L) est la matrice diagonale des valeurs singulières positives classées par ordre croissant (i.e. I.J(1-/ii) et V=(v1) sont deux matrices carrées LxL unitaires réelles et VT désigne ici la transposée de V. Pour tout voxel Vi d'intérêt, un procédé conforme à l'invention consiste à mettre en oeuvre une étape 130 pour produire dIF pour in fine produire l'estimation di=diF en 114 d'une quantité d'intérêt d. Dans le cas de la méthode cSVD, cette production de OF est préférentiellement réalisée en fixant une fois pour toute un entier IF tel que 0/F Enfin, l'invention ne s'applique pas seulement à l'imagerie de perfusion mais à tout type de données dont le traitement doit être effectué y compris en moindre temps.25 REVENDICATIONS1. Procédé mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie médicale pour produire une estimation d'un paramètre hémodynamique (14) à partir de signaux d'intensité expérimentaux S(t), ladite estimation étant réalisée à partir d'une estimation d d'une quantité d'intérêt d d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe, ladite estimation consistant à calculer â=À-1-c, À-1 étant une pseudo-inverse d'une matrice de convolution A et c décrivant une courbe de concentration d'un agent de contraste dans ledit voxel, ladite courbe de concentration résultant d'une conversion préalable desdits signaux expérimentaux, caractérisé en ce qu'il comporte : - une étape (102) pour décomposer canoniquement A en valeurs singulières sous la forme A = U - S - VT où S -=--- diag(c 7 1,...,c y L) est la matrice diagonale des valeurs singulières classées par ordre croissant, VT est la transposée d'une matrice V, V=(l'Y ..) Li et =(uii) sont deux matrices carrées de dimension LxL unitaires et réelles, L>N, N étant un nombre d'échantillons déterminé, d, A et c étant respectivement de dimensions Lxl, LxL et Lx1 ; - au moins une étape itérative (130) pour :o élaborer (130, 131) une pseudo-inverse -1 -1 Al de A, sous la forme Al =V-WI-UT, UT désignant la transposée de U, avec ^ Wi= diag ^-r 0 2. Procédé selon la revendication précédente, selon lequel il comporte l'étape (131) pour produire -1 Al , 0 3. Procédé mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie médicale pour produire une estimation d'un paramètre hémodynamique (14) à partir de signaux d'intensité expérimentaux S(t), ladite estimation étant réalisée à partir d'une estimation d d'une quantité d'intérêt d d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe, ledit procédé étant mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie médicale, ladite estimation consistant à calculer a=À-1 .c, étant une pseudo-inverse d'une matrice de convolution A etc décrivant une courbe de concentration d'un agent de contraste dans ledit voxel, caractérisé en ce qu'il comporte : une première étape (102) pour décomposer canoniquement A en valeurs singulières sous la forme A=U - S -VT où S=diag(ai,...05 L) est la matrice diagonale des valeurs singulières positives classées par ordre croissant, VT est la transposée d'une matrice V, V=(vii) et Huii) sont deux matrices carrées de dimension LxL unitaires et réelles, L>AT, étant un nombre d'échantillons déterminé, d, A et c étant respectivement de dimensions Lxl, LxL et Lxl ; au moins une étape itérative (130) pour, à l'itération 0 4. Procédé selon la revendication précédente, selon lequel il comporte une étape (110) pour produire d0 = AT4fc .30 5. Procédé selon l'une quelconque des revendications 1 à 4, pour lequel il comporte préalablement à chaque étape pour produire (130, 132) d1 (130, 132) la vérification (112) d'une condition de régularisation, ledit procédé s'interrompant à l'itération /F dès que ladite condition est satisfaite, l'estimation (114) de d étant a=c/iF. 6. Procédé selon la revendication précédente, pour lequel la condition de régularisation consiste en le calcul (111) d'un index d'oscillation 0//, ladite condition de régularisation étant satisfaite dès que (112) OliF_POI, PCU étant un seul prédéterminé. 7. Procédé selon la revendication précédente, pour lequel l'index d'oscillation est calculé tel que CU/= 1 1E (i) 2d/ + d/(i-2)1 . L max dl (i) . i=1,L 8. Procédé selon l'une quelconque des revendications précédentes, selon lequel la matrice de convolution A est bloc-circulante de dimension LxL définie (101) par Ca(ti) 0 -- Ca(tN) Ca(t2) Ca(ti) A- At - Ca(tN) Ca(tN) At 0 0 0 0 Ca(tN) - Ca (4) ) étant la période d'échantillonnage, Ca(t) étant une concentration de l'agent de contraste. 9. Unité de traitement (4) comportant des moyens de mémorisation, des moyens pour communiquer avec le monde extérieur et des moyens de traitement, caractérisée en ce que : - les moyens pour communiquer sont aptes à recevoir du monde extérieur une donnée expérimentale c décrivant une courbe de concentration d'un agent de contraste dans un volume élémentaire - dit voxel - d'un organe ; - les moyens de traitement sont adaptés pour mettre en oeuvre un procédé pour produire une estimation d'un paramètre hémodynamique selon l'une quelconque des revendications 1 à 8. 10. Unité de traitement selon la revendication précédente, pour lequel les moyens pour communiquer délivrent l'estimation d'un paramètre hémodynamique selon un format approprié à une interface homme-machine (5) apte à la restituer à un utilisateur (6). 11. Système d'analyse d'imagerie médicale comportant une unité de traitement (4) selon la revendication précédente et une interface homme- machine (5) apte à restituer à un utilisateur (6) un paramètre hémodynamique estimé par ladite unité de traitement selon un procédé conforme à l'une quelconque des revendications 1 à 8.30