Brevet FR2965951A1 — publié le 20120413

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 par application de méthodes probabilistes douces à l'imagerie de perfusion. Un tel procédé permet en outre d'estimer des fonctions de répartition complémentaire ou d'entrée artérielles et donc plus généralement toute quantité d'intérêt. L'invention se distingue notamment des procédés connus en ce quelle requière l'introduction d'informations a priori douces de nature physiologique ou hémodynamique sans contraindre ni forcer l'estimation recherchée par des hypothèses arbitraires et indésirables. L'invention s'appuie notamment sur des techniques 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 d'obtenir rapidement 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 en urgence de pathologies telles que les accidents vasculaires cérébraux. Pour mettre en oeuvre de telles techniques, on utilise un appareil 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 d'une partie du corps, notamment du cerveau. Ledit appareil 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. Cette unité de traitement délivre in fine à un praticien, une estimation des paramètres hémodynamiques à partir des images de perfusion, au moyen d'une interface homme-machine 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 S,,y,M 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 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. 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 signal pour chaque voxel est directement proportionnel à la concentration : S(t)=k-C(t)+50. En Imagerie de Perfusion par Résonance Magnétique Nucléaire, il existe k~T&C(t~ une relation exponentielle S(t)=S0-e- . Dans les deux cas, SD représente l'intensité moyenne du signal avant l'arrivée de l'agent de contraste. En ce qui concerne 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. De manière générale, nous noterons S(t)=W(C(t),Os) le modèle reliant le signal théorique S(t) à la concentration de l'agent de contraste théorique C(t), es étant le vecteur des paramètres libres de ce modèle. Par exemple, dans les cas de l'imagerie de 30 perfusion par résonance magnétique et de l'imagerie de perfusion par tomodensitométrie, Os=(Sp). La conservation de la masse de l'agent de contraste dans le volume de tissu contenu dans chaque voxel à chaque instant s'écrit '(t) -BF.[CQ (t)-C' (t)]. Ça (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) e s t la concentration de l'agent de contraste dans la veine drainant le volume de tissu (fonction de sortie veineuse ou Venous Output Fonction - VOF - selon une terminologie anglo-saxonne). En supposant le système dynamique artère/tissu/veine linéaire et invariant dans le temps, on peut écrire C'(t)=Ca(t)Oh(t) où h(t) est la réponse impulsionnelle du système - ou encore fonction densité de probabilité du temps de transit de l'agent de contraste dans le tissu - et OO désigne le produit de convolution. Une solution formelle de l'équation différentielle précédente avec condition initiale 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 t R(t)=H(t)- f h(z)d'c où H est la fonction généralisée 0 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 anglo-saxonne) . +co +co MTT= f t.h(t)dt= f R(t)dt (si tlimcot.h(t)=0). 0 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 . Les paramètres hémodynamiques tels que BF, MTT ou BV ainsi que la fonction de répartition complémentaire R(t) sont actuellement estimés comme suit, dans le cas de l'imagerie de perfusion par résonance magnétique nucléaire. Pour chaque voxel, le signal de perfusion expérimental SeXp(t) échantillonné aux instants de mesure tl,i=1,N, est converti en une courbe de concentration C(t) par la relation : Vi=1,N C(ti)=-k TEln Se~~t`~ La 0 constante k est fixée pour tous les voxels à une valeur arbitraire non nulle (e.g. k-TE =1). La constante SD est estimée en prenant, par exemple, sa moyenne avant que l'agent de contraste n'arrive. Notons que cela n'est possible que si l'acquisition des signaux de perfusion démarre suffisamment tôt par rapport au temps d'arrivée de l'agent de contraste (Bolus Arrivai Time - BAT - selon une terminologie anglo-saxonne). A partir de la courbe de concentration C(t), et en supposant connue la fonction d'entrée artérielle théorique associée Ça (t), le produit BF -R(t) est estimé par déconvolution numérique. Diverses approches ont été envisagées pour obtenir 30 les fonctions d'entrée artérielles théoriques Ca(t) afin de déconvoluer les courbes de concentration C(t). Dans une première approche, une fonction d'entrée a r t é r i e l l e globale expérimentale est choisie manuellement par le praticien. Elle peut être mesurée, par exemple, au niveau de l'artère sylvienne controlatérale ou de l'artère carotide interne pour l'imagerie de perfusion du cerveau, ou obtenue par des mesures supplémentaires, par exemple optiques. Si elle permet d'obtenir des signaux avec hauts rapports signal sur bruit, cette approche a néanmoins de nombreux inconvénients. Tout d'abord, elle nécessite une intervention humaine et/ou des mesures supplémentaires. Ceci n'est pas souhaitable en situation d'urgence clinique et rend les procédures et les résultats finaux plus difficilement reproductibles. Ensuite et surtout, cette fonction d'entrée artérielle globale ne correspond pas aux fonctions d'entrée artérielles locales de chaque voxel. Elle en diffère en termes de délai (car les fonctions d'entrée artérielles locales sont en général en retard par rapport à la fonction d'entrée artérielle globale prise en amont du système vasculaire) et de dispersion (car la propagation de l'agent de contraste est plus lente en aval qu'en amont du système vasculaire). Or on sait que ces phénomènes ont au final une incidence considérable sur les estimations des paramètres hémodynamiques puisque, par symétrie du produit de convolution, ces défauts se répercutent directement dans l'estimation de la fonction de répartition complémentaire. Ainsi, par exemple, on n'obtient pas au final une estimation du véritable temps de transit moyen (MTT ) entre la fonction d'entrée artérielle locale et la fonction de sortie veineuse locale mais seulement un temps de transit moyen entre la fonction d'entrée artérielle globale et la fonction de sortie veineuse. Afin de pallier ces discordances, certains auteurs ont introduit de nouveaux paramètres descriptifs tels que le paramètre TMAX=argmaxR(t) quantifiant le délai t entre la fonction d'entrée artérielle globale et les fonctions d'entrée artérielles locales, quand bien même ils n'interviennent pas dans le modèle de perfusion standard original (dans lequel la fonction d'entrée artérielle est la véritable fonction d'entrée artérielle locale pour chaque voxel). D'autres méthodes tendent à minimiser l'influence de ces discrépances des fonctions d'entrée artérielles locales sur l'estimation des paramètres hémodynamiques. Elles introduisent toutefois de nouvelles inconnues dans le problème global et ne font que l'éluder. Selon une deuxième approche, une fonction d'entrée artérielle globale est obtenue automatiquement à partir d'images de perfusion via des techniques du signal telles que le partitionnement l'analyse en composantes indépendantes Component Analysis - ICA - selon une anglo-saxonne). Si cette approche de traitement de données ou (Independent terminologie permet de s'affranchir de l'intervention humaine, elle ne résout pas les problèmes de délai et de dispersion inhérents aux fonctions d'entrée artérielles globales et introduit de nouvelles inconnues (e.g. il se peut qu'on obtienne des fonctions de sortie veineuses au lieu des fonctions d'entrée artérielles). Selon une troisième approche, des fonctions d'entrée artérielle locales sont obtenues automatiquement à partir des images de perfusion à l'aide de techniques de traitement du signal et de critères de sélection. On recherche par exemple la « meilleure » fonction dans le voisinage immédiat du voxel tissulaire courant où l'on souhaite estimer les paramètres hémodynamiques ou les fonctions de répartition complémentaires. L'objet de cette troisième approche est d'obtenir au final des estimations moins biaisées et plus précises en s'affranchissant, du moins en partie, des problèmes de délai et de dispersion. Cependant, rien ne garantit a priori et a posteriori que les fonctions d'entrée artérielle locales ainsi obtenues sont des approximations pertinentes de la « véritable fonction locale » pour le voxel d'intérêt. Par exemple, cette « véritable » fonction pourrait ne pas se situer dans le voisinage considéré (s'il est trop petit) ou, au contraire, pourrait être confondu avec une autre fonction d'entrée artérielle (s'il est trop grand). De plus, cette « meilleure » fonction d'entrée artérielle locale est recherchée parmi des fonctions d'entrée artérielle « normales » (i.e. avec un temps d'arrivée de l'agent de c o n t r a s t e de courte/grande précocité, de grande amplitude, etc.). Or il s'agit précisément de distinguer fonctions d'entrée artérielles normales et fonctions d'entrée artérielles pathologiques, par exemple ischémiques. Par conséquent, même si les résultats finaux peuvent être meilleurs qu'avec une approche globale, les incertitudes sur ces fonctions d'entrée artérielles locales et, a fortiori, sur les paramètres hémodynamiques ou les fonctions de répartition complémentaires demeurent en grande partie. Pour effectuer l'opération de déconvolution de la courbe de concentration expérimentale C(t) par la fonction d'entrée artérielle Ça (t) théorique donnée par les méthodes décrites précédemment, le modèle standard de convolution C(t)=BF-Ca(t)OR(t) est tout d'abord discrétisé temporellement, par exemple, suivant l'approximation de la méthode des rectangles : ti i Vi=l, N, C(ti)=BF. f ça(z) R(t-z)d'r~BF Ot ~Ca(ti) R(ti tk) 0 k=0 où At est la période d'échantillonnage. On se ramène 5 alors à un système linéaire Ad =c en posant Ça (0 0 0 \ R(tl) A=Ot. Ca(t2) Ca(0 b= R(t2) 0 ~Ca(tN) Ca(tN-1) Ca(0) R(tN) d = BF.b En pratique, la matrice A est très mal 10 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 aberrantes. Il faut donc recourir à diverses méthodes pour obtenir par exemple une pseudo- 15 inverse A-1 de la matrice A et par suite une estimation d de d par d=A-1.c. Parmi ces méthodes d'obtention d'une pseudo-inverse, on peut citer la décomposition en valeurs singulières tronquée (Truncated Singular Value Decomposition, - (T) SVD - ou 20 Simple Singular Value Decomposition, -sSVD - selon des terminologies anglo-saxonnes) et la déconvolution de Hunt dans le domaine fréquentiel. De manière plus générale, on peut minimiser un 25 critère du type Ad -c 2 + rd 2 où rd 2 est un terme de régularisation privilégiant certaines solutions et obtenir une estimation de d par d= argmin(Ad-c 2+ rd 2) . Parmi ces méthodes, nous pouvons citer la c= régularisation de Tikhonov, l'utilisation de la transformée en ondelettes, etc. Une fois d obtenue, on peut obtenir une estimation BF de BF par BF = d(t1) = d(0) puisque, par définition, lim R (t) =1 . Cependant, on constate que l'on estime N souvent BF par BF =maxd(ti) , par exemple dans le cadre Z=1 des méthodes à base de Décomposition en Valeurs Singulières, afin de compenser les sous-estimations systématiques en d(0) inhérents à ces méthodes. Par suite, on obtient une estimation b de b par b= d BF N puis une estimation de MTT = I R(t)dt par MTT =ot-1b(tz) o 1=1 en suivant, par exemple, l'approximation par la méthode des rectangles. Enfin, on obtient typiquement une estimation de BV par B V = BF.MTT même si l'estimation d'un produit n'est pas le produit des estimations. On trouve de nombreuses variantes de ces méthodes à base de fonction (s) d'entrée artérielle (s) . par exemple, les fonctions d'entrée artérielles expérimentales peuvent être préalablement ajustées à un modèle théorique paramétrique ou semi-paramétrique Ca(t,0a) où Oa est un vecteur de paramètres, afin d'accroître artificiellement le rapport signal sur bruit. Le signal peut être également sur-échantillonné artificiellement afin de rendre la déconvolution numérique plus stable ou encore de s'affranchir des éventuels problèmes dus au phénomène de recirculation, lorsqu'il y a chevauchement temporel du signal de circulation de l'agent de contraste (first pass) et du signal de recirculation (second pass). Mais ces variantes restent basées sur des méthodes de déconvolution numérique telles que la décomposition en valeurs singulières tronquée. D'après certaines études comparées, parmi les différentes méthodes testées sur données synthétiques censées simuler des données réelles typiques, la déconvolution par décomposition en valeurs singulières et ses variantes (décomposition en valeurs singulières linéaire tronquée (sSVD), circulaire tronquée (cSVD) ou circulaire tronquée lissée (oSVD)) donnent au final les meilleures estimations des paramètres BF et MUT en termes de biais (erreur systématique par rapport à la valeur véritable), précision (écart-type de l'estimation par rapport à la valeur véritable en fonction du rapport signal-sur-bruit des signaux d'entrée) et robustesse par rapport aux diverses fonctions de répartition complémentaires R(t) et aux diverses fonctions d'entrée artérielles Ca(t) pouvant se présenter en pratique suivant les patients, les 20 types de tissus, les pathologies, etc. Cependant, outre les problèmes relatifs à la sélection des fonctions d'entrée artérielles décrits précédemment, cette famille de méthodes numériques 25 souffre de graves problèmes inhérents. Tout d'abord, les estimations d de BFb ne sont pas décroissantes au cours du temps mais oscillantes, à tel point qu'elles peuvent parfois prendre des valeurs 30 négatives. Or, R(t) qui est la quantité d'agent de contraste restant dans le voxel à un instant t, est une fonction nécessairement décroissante et positive. Des méthodes ad hoc telles que la oSVD permettent de réduire ces oscillations aberrantes mais elles subsistent car elles sont inhérentes à la décomposition en valeurs singulières. C'est la raison pour laquelle avec cette famille de méthodes, on estime les flux N sanguins BF par BF =maxd(ti) alors que dans le cadre du Z=i modèle standard de perfusion, on devrait au contraire les estimer par BF=d(0) car limR(t)=1. En prenant le t_>o, maximum au lieu de la valeur à l'origine, on espère seulement gommer l'effet de ces oscillations sur les estimations de BF. Celles-ci ne peuvent donc être parfaitement satisfaisantes. En particulier, la précision des estimations des BF pouvant être atteinte par des méthodes plus rigoureuses d'estimation des fonctions de répartition complémentaires et donc des paramètres hémodynamiques reste inconnue. Ainsi, ces méthodes de déconvolution numérique entrent en contradiction avec le modèle standard de la perfusion puisqu'elles fournissent des solutions ne respectant pas les propriétés dudit modèle. Afin de s'affranchir de ce problème et obtenir des estimations de fonctions de répartition complémentaires physiologiquement admissibles, des modèles paramétriques R(t,OR) pour ces fonctions de répartition complémentaires ont été introduits, OR étant le vecteur des paramètres dudit modèle. Ces modèles sont ajustés aux signaux expérimentaux, par exemple par la méthode de Bayes. Cependant, cette approche peut apparaître prématurée à ce jour. En effet, il faudrait préalablement disposer d'estimations non-paramétriques de ces fonctions de répartition complémentaires pour déterminer des modèles paramétriques ou semiparamétriques adéquats à les décrire. Or, des simulations de Monte-Carlo ont montré que si ces modèles ne sont pas parfaitement adaptés pour décrire correctement toutes les formes de fonctions de répartition complémentaires pouvant se présenter en pratique, alors les estimations résultantes des paramètres hémodynamiques tels que MTT ou BF deviennent aberrantes. Le choix des modèles théoriques pour les fonctions de répartition complémentaires est donc critique et il ne peut se faire convenablement qu'en confrontant les modèles à des données expérimentales. D'où la nécessité de disposer de méthodes non-paramétriques pour estimer ces fonctions de répartition complémentaires, afin éventuellement de les remplacer dans un second temps par des méthodes paramétriques ou semi-paramétriques classiques permettant d'obtenir des estimations physiologiquement admissibles des fonctions de répartition complémentaires. Comme mentionné précédemment, les estimations obtenues par des méthodes de déconvolution telles que la décomposition en valeurs singulières sont en contradiction avec la définition même de R(t)=H(t)- f h(z)dz dans le cadre du modèle standard de o perfusion. Elles ne sont pas physiologiquement et physiquement admissibles. Il n'est donc pas possible d'ajuster des modèles théoriques paramétriques ou semi- paramétriques à ces estimations. A fortiori, il n'est pas possible de comparer plusieurs modèles et de sélectionner les plus adaptés à décrire les fonctions de répartition expérimentales. Il n'est donc plus possible de progresser dans la modélisation et la compréhension des phénomènes de perfusion à cause des défauts inhérents à de méthodes numériques telles que la décomposition en valeurs singulières. Par ailleurs, il s'avère que le problème sous- jacent à l'estimation non-paramétrique des fonctions de répartition complémentaires et des paramètres hémodynamiques n'est pas un « simple » problème de déconvolution de courbes de concentration empiriques par des fonctions d'entrée artérielles empiriques. En effet, ce pourrait être le cas si les fonctions d'entrée artérielles étaient réellement des données du problème, connues avec une certitude absolue et une précision infinie. Or on ne dispose au mieux que des signaux artériels expérimentaux mesurés et donc connus seulement aux bruits de mesure près, qu'il faut de plus convertir préalablement en courbes de concentration. Autrement dit, en supposant les fonctions d'entrée artérielles empiriques mesurées égales aux fonctions d'entrée artérielles théoriques, on néglige les bruits de mesure sur les signaux expérimentaux et les incertitudes provenant de l'estimation ou de la conversion des courbes de concentration à partir de ces signaux. Le modèle standard de convolution C(t)=BF.Ca(t)OR(t) porte uniquement sur des signaux théoriques qui ne peuvent pas, en général, être mesurés directement. En réalité, les signaux mesurés sont la somme des signaux théoriques et du bruit de mesure. Ainsi, les signaux de perfusion tissulaire et artériel s'écrivent respectivement Sem, (t) = Sth (t)+fit et SaeXp (t) = Sath (t)+ fit où t et e sont des processus stochastiques d'espérances mathématiques nulles modélisant les bruits de mesure. On a de plus Sth(t)=SDe Cth(t) et Sath(t)=SDe Cach(t) dans le cas de l'imagerie de perfusion par résonance magnétique ou de manière plus générale Sth(t)=Y(Cth(t),Os) comme décrit précédemment. Le modèle standard de perfusion théorique appliqué aux signaux expérimentaux mesurés par imagerie de perfusion doit donc s'écrire Cth(t)=BF-Cath(t)xOR(t), soit, dans le cas de l'imagerie de perfusion par Résonance magnétique nucléaire : ln Sexp Vit) - t = BF - ln s w -e OR(t) 0 Oa ou de manière équivalente : Sexp (t)+t Saexp (t)+e ln =BF-ln OR(t) SO SOa non pas . CeXp(t)=CBF-LCaeXp (t)j0R(t)+t qui est le modèle mathématique de convolution implicite erroné sur lequel se basent la plupart des méthodes de déconvolution. On voit d'ailleurs qu'il est préférable d'écrire le modèle standard de perfusion sous la forme : iSexp(t)=s0e Cth(t)+t Saexp (t)= SOae-Ceh(t) +st Cth ( t ) = BF.Cath (t) O R(t) afin d'éviter de devoir prendre le logarithme de variables aléatoires qui ne sont pas toujours positives. Il est connu que les incertitudes de mesure sur le signal à déconvoluer ont une influence considérable sur le résultat final de l'opération de déconvolution : une infime variation sur le signal d'entrée due à ces incertitudes peut engendrer une variation considérable sur le résultat final. C'est d'ailleurs pour pallier ces problèmes et réduire ces instabilités que des méthodes de déconvolution telles que la régularisation de Tikhonov ou la Décomposition en Valeurs Singulières ont été introduites. A fortiori, l'influence des bruits de mesure et des incertitudes sur les fonctions d'entrée artérielles, aujourd'hui entièrement négligée, est encore plus considérable : le bruit de mesure artériel e et les incertitudes sur SDa interviennent désormais dans la matrice de convolution A et sont donc propagées et amplifiées. Cette négligence des erreurs et des incertitudes sur les fonctions d'entrée artérielle provoque des erreurs importantes sur l'estimation des paramètres hémodynamiques, en même temps qu'une illusion de précision sur ces estimations. Certaines méthodes visent à gommer les bruits de mesure dans les fonctions d'entrée artérielles afin de minimiser ce problème qui reste éludé à ce jour. Il serait préférable de disposer de méthodes permettant de propager les incertitudes sur les fonctions d'entrée a r té r i e l l e s su r l'estimation des paramètres hémodynamiques et des fonctions de répartition complémentaires afin de maîtriser et de quantifier les erreurs d'estimation. D'autre part, il conviendrait également d'éviter d'écrire le modèle standard de la perfusion comme Cth(t)=BF.Cath(t)OR(t). En effet, le modèle standard de la perfusion définit en premier lieu la réponse 30 impulsionnelle hW du système dynamique artère/tissu/veine à partir de laquelle la fonction de répartition complémentaire R(t) est calculée : R(t)=H(t)- f h(z)dz . Il convient donc d'écrire le modèle 0 standard de la perfusion théorique en fonction de la réponse impulsionnelle h(t) comme étant t Cth(t)=BF.Cath(t)XO H(t)- f h(z)dz , d'ajuster ce modèle afin 0 d'estimer la réponse impulsionnelle h(t) aux instants 5 de mesure tj =1,N par h pour obtenir ensuite une estimation R de la fonction de répartition H(ti)-Ot.lh(ti) i=2 (méthode d'approximation des rectangles). En particulier, mieux vaut ne pas estimer la réponse 10 impulsionnelle h(t) à partir de l'estimation de la fonction de répartition complémentaire R. Du point de vue numerique, on comprend aisément qu'il est préférable d'estimer la dérivée (h(t)) pour estimer ensuite la primitive (R(t)) plutôt que l'inverse. On 15 constate néanmoins que la fonction de répartition complémentaire R(t) est estimée directement au lieu de la réponse impulsionnelle h(t). D'autre part, le problème de la déconvolution est 20 mal posé admettant a priori une infinité de solutions possibles. A fortiori, le problème mal posé auquel l'on fait face en imagerie de perfusion pour estimer des paramètres hémodynamiques, des réponses impulsionnelles, des fonctions de répartition 25 complémentaires ou des fonctions d'entrée artérielles, est un problème de déconvolution triplé d'un problème de propagation des bruits de mesure et d'un problème de conversion des signaux expérimentaux en courbes de concentration. 30 Pour tenter de se ramener à un problème bien posé admettant éventuellement une unique solution, on doit complémentaire par, par exemple, R M = rajouter a priori de l'information, des contraintes sur la solution recherchée parmi l'infinité de solutions possibles a priori. C'est la raison pour laquelle, on trouve de nombreuse méthodes de déconvolution et d'estimation dans l'état de la technique, chaque méthode injectant, de manière plus ou moins explicite, ou plus ou moins directe, un type particulier d'information a priori. Ainsi et à titre d'exemple, nous pouvons citer la régularisation de Tikhonov classique telle que décrite précédemment. Dans sa version la plus courante, la matrice F vaut r = aIN. Ceci pénalise les solutions de grande norme euclidienne, le scalaire a quantifiant le poids de cette pénalisation, de cette contrainte. Pour sa part, la Décomposition en Valeurs Singulières Tronquée, classique en imagerie de perfusion, consiste à annuler les valeurs singulières de la matrice de convolution plus petites qu'un certain seuil qui joue donc le rôle de paramètre de régularisation. Ces petites valeurs singulières sont reliées aux composantes du signal de plus hautes fréquences, de sorte que la Décomposition en Valeurs Singulières Tronquée agit comme un filtre passe-bas. Si les fonctions de répartition complémentaires sont certainement des signaux de basse fréquence, la troncature des petites valeurs singulières ne correspond cependant pas directement et seulement à un ajout d'information de nature physiologique mais surtout à un ajout d'information de nature purement algébrique, en l'occurrence l'appartenance à un sous-espace de dimension donnée. On peut constater à ce titre que des oscillations de trop hautes fréquences demeurent. Les signaux estimés ont en effet tendance à suivre les bruits de mesure. De plus, le seuil de troncature ne correspondant pas au poids d'une contrainte explicite (mais seulement à forcer les solutions à appartenir à un certain sous-espace), il est difficile de déterminer des critères pour déterminer sa valeur. Réciproquement, il n'est pas garanti que l'on puisse convenablement optimiser un tel critère, par exemple la rugosité de la fonction de répartition complémentaire dans le cas de la Décomposition en Valeurs Singulières Tronquée Lissée (oSVD), par Décomposition en Valeurs Singulières Tronquée car l'action de la méthode sur le critère est très indirecte. Les méthodes actuelles d'estimation de paramètres hémodynam iques ou de fonction de répartition complémentaires ne permettent pas en outre d'obtenir une estimation de la précision sur ces estimations (i.e. écart-type de l'estimateur) et encore moins de quantifier la confiance qu'on peut accorder aux estimations et aux estimations de la précision sur ces estimations. En particulier, il est difficile de quantifier la qualité de l'ajustement (Goodness of fit, selon une terminologie anglo-saxonne) du modèle standard de la perfusion par les méthodes de déconvolution et d'estimation. On peut constater ainsi que les méthodes d'estimation à base de Déconvolution en Valeurs Singulières ont par exemple tendance à sur-ajuster (overfitting, selon une terminologie anglo- saxonne) les signaux expérimentaux. Les signaux estimés ne sont pas lisses mais, au contraire, ont tendance à suivre le bruit de mesure. On peut obtenir des valeurs faibles pour mesurer l'ajustement des modèles (statistique du x2 - Sum of Squared Errors, selon une terminologie anglo-saxonne), mais ces faibles valeurs sont trompeuses et la statistique du x2 devrait être écartée si on veut détecter et prendre en compte les phénomènes de sur-ajustement. Il n'est donc pas évident d'introduire de « bonnes » mesures d'ajustement pour les méthodes de déconvolution et d'estimation des paramètres hémodynamiques suivant l'état de la technique. Il est donc souhaitable de disposer de méthodes pour lesquelles une telle mesure d'ajustement est définie de manière univoque, prend en compte le sur-ajustement et permet de quantifier seulement l'adéquation des données au modèle standard de la perfusion. En résumé, les méthodes actuelles d'estimation des paramètres hémodynamiques, des réponses impulsionnelles ou des fonctions de répartition complémentaires sont empreintes de nombreuses erreurs méthodologiques et de nombreuses approximations non maîtrisées. L'interprétation des résultats est ainsi souvent rendue difficile et l'imagerie de perfusion n'est pas pleinement exploitée. L'objet de la présente invention est de répondre à l'ensemble des inconvénients engendrés par l'utilisation des méthodes connues. L'invention a pour principal objet de fournir de nouvelles méthodes permettant de rechercher l'estimation h de la réponse impulsionnelle h parmi des solutions vérifiant certaines contraintes de nature physiologique et hémodynamique, sans pour autant introduire des contraintes ad hoc non forcement vérifiées et surtout invérifiables par l'expérience ou de nature non physiologique ou hémodynamique. De plus, ces méthodes permettent de déterminer l e s poids de telles contraintes de manière automatique et univoque, sans avoir à recourir à des méthodes ad hoc. L'invention permet ainsi de prétendre que le problème de l'estimation des paramètres hémodynamiques, des fonctions de répartition complémentaires ou des fonctions d'entrée artérielles, est à présent bien posé. L'invention permet d'y apporter la seule et unique solution (éventuellement multiple). Parmi les principaux avantages conférés par l'invention, nous pouvons citer, de manière non limitative, le fait de pouvoir : - traduire quantitativement des informations physiologiques qualitatives ou semi-qualitatives douces sur les paramètres hémodynamiques, les réponses impulsionnelles, les fonctions de répartition complémentaires ou les fonctions d'entrée artérielles ; - prendre en compte explicitement les incertitudes et les erreurs sur les fonctions d'entrée artérielles et sur les signaux expérimentaux et, les propager sur les estimations des quantités d'intérêt ; - améliorer l'estimation desdits paramètres, des réponses impulsionnelles, des fonctions de répartition complémentaires ou des fonctions d'entrée artérielles en termes de biais (erreur systématique), de précision (erreur statistique) et de robustesse vis-à-vis des différentes configurations pouvant se présenter dans la pratique ; - obtenir des intervalles de confiance - et même des paris sur lesdits intervalles de confiance - sur les estimations des paramètres hémodynamiques, des réponses impulsionnelles, des fonctions de répartition complémentaires ou des fonctions d'entrée artérielles, pour préciser et augmenter la c o n f i a n c e que l'on peut accorder à ces estimations ; obtenir des estimations non-paramétriques des réponses impulsionnelles, des fonctions de répartition complémentaires ou des fonctions d'entrée artérielles admissibles du point de vue physiologique et hémodynamique et compatible avec le modèle standard de la perfusion, pour permettre dans une étape ultérieure l'ajustement, la comparaison et enfin la sélection de modèles théoriques paramétriques ou semi-paramétriques pour lesdites réponses impulsionnelles, lesdites fonctions de répartition complémentaires ou lesdites fonctions d'entrée artérielles ; quantifier dans quelle mesure les fonctions d'entrée artérielles en entrée peuvent être effectivement des fonctions d'entrée artérielles admissibles pour le tissu dans chaque voxel d'intérêt ; fournir des mesures objectives et quantitatives de l'adéquation d'un modèle de perfusion global, permettant au final la comparaison et la sélection des modèles de perfusion globaux les plus appropriés. A cette fin, il est prévu un procédé pour estimer une quantité d'intérêt parmi une pluralité d'un système dynamique artère/tissu/veine d'un volume élémentaire dit voxel - d'un organe. Un tel procédé est destiné à être mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie de perfusion, et comporte une étape pour estimer ladite quantité d'intérêt à partir d'une donnée de perfusion expérimentale. Selon l'invention, ladite étape consiste à évaluer, suivant la méthode de Bayes, une distribution marginale a posteriori pour ladite quantité d'intérêt par : - l'assignation d'une distribution de probabilité directe de la donnée de perfusion sachant les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ; - l'assignation d'une distribution de probabilité conjointe a priori desdites quantités, par l'introduction d'une information douce relative à la réponse impulsionnelle dudit système dynamique. L'invention prévoit également un procédé pour estimer une quantité d'intérêt parmi une pluralité d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe, ledit système dynamique étant linéaire, invariant dans le temps et déterminé formellement par la relation C(t)=BF.Ca(t)OR(t) où 0) est la concentration d'un agent de contraste circulant dans un voxel, Ca(t) est la concentration dudit agent de contraste dans l'artère alimentant ledit voxel, BF est le flux sanguin dans ledit voxel, O désigne le produit de convolution et R(t) est la fonction de répartition complémentaire du temps de transit dans ledit voxel. Un tel procédé est mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie de perfusion, et comporte une étape pour estimer ladite quantité d'intérêt à partir d'une d o n n é e de perfusion expérimentale. Selon l'invention, ladite étape consiste à évaluer, suivant la méthode de Bayes, une distribution marginale a posteriori pour ladite quantité d'intérêt par : - l'assignation d'une distribution de probabilité directe de la donnée de perfusion sachant les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ; - l'assignation d'une distribution de probabilité conjointe a priori desdites quantités, par l'introduction d'une information douce relative à la fonction de répartition complémentaire du temps de transit R(t) dans ledit voxel. L'invention prévoit que l'assignation d'une distribution de probabilité conjointe a priori desdites 20 quantités puisse être réalisée en outre par l'introduction d'une information douce relative à la concentration dudit agent de contraste dans l'artère alimentant le voxel. 25 L'invention prévoit que de tels procédés puissent être mis en oeuvre par itérations successives pour une pluralité de voxels considérés. Selon certaines variantes de réalisation, un 30 procédé conforme à l'invention peut comporter : - une étape pour calculer une information complémentaire sous la forme d'un intervalle de confiance associé à une quantité d'intérêt estimée ; 10 15 - une étape pour calculer une information complémentaire sous la forme d'un taux de pari sur un intervalle de confiance associé à une quantité d'intérêt estimée ; - une étape pour calculer une information complémentaire sous la forme d'une mesure de l'adéquation du produit de : - l'assignation de la distribution de probabilité directe de la donnée de perfusion sachant les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ; - l'assignation de la distribution de probabilité conjointe a priori desdites quantités. L'invention prévoit que de tels procédés puissent comporter une étape pour délivrer une quantité d'intérêt estimée - voire toute information complémentaire associée à ladite quantité d'intérêt estimée - à une interface homme-machine apte à la restituer à un utilisateur. Selon un second objet, l'invention prévoit une unité de traitement comportant des moyens de mémorisation, des moyens pour communiquer avec le monde extérieur et des moyens de traitement. Les moyens pour communiquer sont aptes à recevoir du monde extérieur une donnée de perfusion expérimentale et les moyens de traitement sont adaptés pour mettre en oeuvre un procédé pour estimer une quantité d'intérêt conforme à l'invention. Les moyens pour communiquer d'une telle unité de traitement peuvent - selon un mode de réalisation préféré - délivrer une quantité d'intérêt estimée selon un format approprié à une interface homme-machine apte à la restituer à un utilisateur. L'invention prévoit en outre un système d'analyse d'imagerie de perfusion comportant une unité de traitement conforme à l'invention et une interface homme-machine apte à restituer à un utilisateur une quantité estimée selon un procédé également conforme à l'invention. 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 lesquelles : - les figures 1 et 2 présentent deux variantes de 20 réalisation d'un système d'analyse d'imagerie de perfusion ; - 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, 25 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 30 perfusion S@) par Résonance Magnétique Nucléaire typique relatif à un voxel d'un cerveau humain ; - la figure 6 présente une courbe de concentration Qt) 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 ; - la figure 8 présente un procédé conforme à l'invention ; - la figure 9 présente une carte relative à des volumes sanguins cérébraux estimés conformément à l'invention ; - la figure 10 présente une carte relative à des flux sanguins cérébraux - en cas d'ischémie cérébrale - estimés conformément à l'invention ; - la figure 11 présente une carte relative à des temps de transit moyen estimés conformément à l'invention ; - la figure 12 concerne une carte relative à une probabilité qu'un flux sanguin cérébral soit dans un intervalle de confiance. La figure 1 permet de présenter un système d'analyse d'images de perfusion. 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 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, I2, , In à des instants ti, t2,..., ti,..., tn comme le 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 -k.TE.C(t) une relation exponentielle S(t)= Swe où SD 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 décrit un exemple de procédé - conforme à l'invention - pour estimer une quantité d'intérêt p a r m i u n e pluralité d'un système dynamique artère/tissu/veine d'un voxel d'un organe. Un tel procédé peut être mis en oeuvre par une unité de traitement d'un système d'analyse d'imagerie de perfusion tel que le système décrit en liaison avec les figures 1 ou 2 et adapté en conséquence. Un procédé conforme à l'invention comporte principalement une étape 56 pour assigner une ou plusieurs distributions marginales a posteriori de différentes quantités d'intérêt que l'on cherche à estimer, telles que des paramètres hémodynamiques, les valeurs de la réponse impulsionnelle théorique aux instants de mesure ou celles de la fonction de répartition complémentaire. Il comporte en outre une étape 57 pour calculer ladite estimation. Pour assigner une telle distribution marginale à posteriori, il est nécessaire de configurer 50 l'unité de traitement. Cette configuration peut être réalisée de préférence par l'unité de traitement elle-même, à partir d'un ou plusieurs paramètres de configuration. La configuration peut également se traduire par la constitution d'une bibliothèque d'une ou plusieurs distributions marginales a posteriori, bibliothèque préétablie et mémorisée dans une mémoire de programme de ladite unité. L'invention prévoit que ladite bibliothèque peut être enrichie au fur et à mesure de son utilisation voire délivrée par une unité externe de calcul apte à réaliser ladite configuration à partir du ou des paramètres de configuration et apte à coopérer avec l'unité de traitement pour délivrer ladite bibliothèque. Un procédé conforme à l'invention peut comporter ainsi des étapes de configuration mises en oeuvre préalablement à l'assignation 56, parmi lesquelles sont nécessaires et suffisantes . - l'assignation 54 de la distribution directe de probabilité des signaux expérimentaux sachant tous les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ; - l'assignation 53 de la distribution de probabilité conjointe a priori de tous ces paramètres. L'invention permet d'estimer une ou plusieurs quantités d'intérêt dans des cas d'application divers : - les fonctions d'entrée artérielles théoriques sont supposées connues avec une certitude absolue et une précision infinie ; - les fonctions d'entrée artérielles locales diffèrent de la fonction d'entrée artérielle globale donnée notamment par un délai temporel inconnu qu'il s'agira d'estimer également ; - le signal d'entrée artériel est seulement mesuré, éventuellement à un délai temporel près ; - les fonctions d'entrée artérielles ne sont pas données et les signaux d'entrée artériels ne sont pas mesurés - cas d'espèce a priori le plus proche du réel. Les étapes de configuration peuvent dépendre du cas d'application considéré. La figure 8 permet d'illustrer un procédé conforme à l'invention selon un premier exemple d'application pour lequel un signal d'entrée artériel est mesuré, éventuellement à un délai temporel z près. Un modèle de perfusion expérimental M peut s'écrire comme : =11(c,os)+ =11(a,Osa)+a = BF - Ab(t-'c) = BF -B(t-'c)a T b=[R(t1),...,R(tN)] est un vecteur des valeurs de la fonction de répartition complémentaire théorique inconnue. Ce vecteur peut être éventuellement exprimé à partir du vecteur des valeurs de la réponse impulsionnelle h=[h(0=0,...,h(tN)IT P a r exemple, N-1 -T b= 1,1,1-At-h(t2),...,1-At- h(tN) en suivant la méthode i=1 d'approximation des rectangles à gauche ou N -T b= 1,1-At-h(t2),...,1-At-lh(tN) en suivant la méthode i=2 d'approximation des rectangles à droite ; T c=[C(t1),...,C(tN)] est un vecteur des valeurs de la concentration théorique inconnue de l'agent de contraste dans le voxel ; T s=[S(t1),...,S(tN)] est un vecteur de mesures expérimentales du signal d'intensité par imagerie de perfusion, réel ou complexe ; OS est un vecteur des paramètres reliant S(t) à CO ; T =[(t1),...,(tN)] est un vecteur des bruits de mesures sur ledit signal mesuré ; sa = [sa (t, ),..., sa (t,)]T est un vecteur des mesures du signal d'entrée artériel expérimental Sa(t) ; Osa est un vecteur des paramètres reliant S,(t) à Ça (t) ; 5a est un vecteur des bruits de mesures sur ledit signal artériel ; a = LCa (q),...,Ca T est un vecteur des valeurs de la fonction d'entrée artérielle théorique Ca(t) inconnue ; A est une matrice de convolution associée au vecteur 10 des valeurs de la fonction d'entrée artérielle théorique inconnue a obtenue en approchant numériquement l'intégrale de convolution par la méthode d'approximation des rectangles, des trapèzes, ou par des méthodes d'ordre supérieur telles que les méthodes 15 de Simpson, de Boole, de Gauss-Legendre, etc. A peut être également une matrice de convolution circulante, telles que celles utilisées dans le cadre de la méthode de dé convolution par décomposition en valeurs singulières circulaire (circular Singular Value 20 Decomposition, cSVD suivant une terminologie anglo-saxonne). Nous noterons O le vecteur des paramètres hémodynamiques du modèle Al, ici par exemple O =(BF,z). 25 En liaison avec les figures 1, 2 et 8, un procédé mis en oeuvre par une unité de traitement 4 peut comporter deux premières étapes de configuration 51d et 51e consistant à introduire respectivement une information Is sur le vecteur des mesures du signal 30 expérimental s et une information 4, sur le vecteur des mesures du signal d'entrée artériel expérimental sa. Au sens de l'invention, on distingue une information dure d'une information douce. Une information dure correspond à toute proposition Booléenne considérée comme certaine - c'est-à-dire dont la probabilité vaut 1. Par exemple les propositions Booléennes telles que « cette courbe est lisse » et « ce signal suit ce modèle » constituent des informations dures. Par opposition, une information douce concerne toute proposition Booléenne qui consiste à indiquer, et à seulement indiquer avec une certaine probabilité de telles propositions Booléennes. Au final, cela revient à introduire des propositions booléennes telles que « cette courbe est plus ou moins lisse » ou « ce signal suit plus ou moins ce modèle ». Si l'on s'en tient, par exemple, seulement aux deux premiers moments EG, a)=(0,0) et EW,e)=(6s,6sa) du couple de vecteurs des bruits de mesure réels correspondant, alors le Principe du Maximum d'Entropie (entropie de Shannon différentielle sous la mesure de référence de Lebesgue) requiert que les vecteurs et a doivent être regardés comme des processus stochastiques mutuellement indépendants, blancs, stationnaires Gaussiens et d'écarts-types respectifs 6s et Usa . De manière plus générale, nous noterons (Es,Esa) les paramètres caractérisant le couple des vecteurs des bruits de mesure G, a ) . Par exemple (Es,Esa)=(o-s,o-sa). Un procédé conforme à l'invention comporte une étape de configuration 54 pour assigner une distribution de probabilité conjointe directe du couple des vecteurs des mesures (s,s,) sachant les vecteurs a et b, le vecteur O, les vecteurs de paramètres Os et Os, et enfin le couple des vecteurs (Es,Esa) . Ladite distribution de probabilité conjointe directe s'écrit alors p(s,sa a,b,O,Os,Osa,Es,Esa,Is,Isa,M) . Par exemple, nous avons : p(s,sa a,b,O,Os,Osa,6s,6sa,Is,Isa,M)oc N [s(ti)-W(BF.Ab(ti -z),es)j2 [sa(ti)- (a(ti),esa )j2 (66a ) N exp 2 + 2 i=1 2o-26a L'invention permet en variante, de pouvoir exprimer différemment ladite distribution de probabilité conjointe directe dans le cas où l'on s'intéresse, non pas directement au vecteur b des valeurs de la fonction de répartition complémentaire théorique inconnue, mais au vecteur des valeurs h de la réponse impulsionnelle. Pour cela, il suffit d'exprimer b en fonction de h. La distribution de probabilité conjointe directe 54 s'écrit alors par exemple p(s,sQ a,h,O,Os,Os ,6s,6s,,Is,Is,M) cc (66a) -Nex N [s(ti)-kf(BF.Ab(ti-2),os)T [sa(ti)-W(a(ti),Osa)T p 26-2 262 a Nous considèrerons par la suite, sans perte de généralité, que l'on s'intéresse préférentiellement au 20 vecteur h. On pourrait en outre assigner de la même manière la distribution de probabilité conjointe directe du couple des vecteurs des mesures de signaux complexes, en multipliant les distributions de probabilité 25 directes de leurs parties réelles et imaginaires. L'invention prévoit également une variante lorsque les vecteurs des valeurs des signaux de perfusion s et sa peuvent être convertis précisément en des vecteurs de valeurs de courbes de concentration ceXp , p a r ceXp(ti) =- 1 lns(ti) k - TE so perfusion par Résonance Magnétique Nucléaire. Il est possible de procéder ainsi, par exemple, lorsque les intensités moyennes SD des signaux de perfusion par résonance magnétique nucléaire avant l'arrivée de 10 l'agent de contraste peuvent être estimées avec précision et indépendamment des autres paramètres. Selon cette variante, on peut montrer qu'il est alors légitime d'assigner une distribution de probabilité gaussienne pour le couple des vecteurs (cexp,cexpa) 15 sachant tous les paramètres telle que : h(cexp,cexpa a,h'e'6'6a,Is,ISa,M)x ~~~ -NeX~ ~cexp(ti)-BF.Ab(012 +[cexpa(ti)-a(012 a p 2 2 i=1 2o-26a 6 et 6a sont désormais les écarts-types des bruits de mesure sur ceXp et cexpa respectivement. Plus généralement, nous utiliserons le terme 20 « donnée de perfusion expérimentale » pour désigner indifféremment un vecteur de valeurs d'un signal de perfusion expérimental s ou sa conversion en un vecteur des valeurs d'une courbe de concentration cm. Dans la suite nous notons, sans perte de généralité, s 25 et sa les données de perfusion expérimentales. Par ailleurs, le procédé comporte trois étapes de configuration 51a, 51b et 51c consistant à introduire respectivement une information Io sur les paramètres hémodynamiques 0 du modèle A!, une information Ih sur exemple bi =1,N en imagerie de la réponse impulsionnelle h et une information la sur la fonction d'entrée artérielle a. Les informations introduites aux étapes 51a à 51e constituent, au sens de l'invention, des paramètres de configuration pour configurer c'est-à-dire - rendre apte une unité de traitement à assigner 56 une distribution marginale à posteriori et ainsi estimer 57 une quantité d'intérêt. A partir de ces informations - ou paramètres de configuration, un procédé conforme à l'invention peut comporter une étape 53 pour assigner une distribution de probabilité conjointe a priori qui s'exprime sous la forme p(a,h,O,es,es,,ES,ES IO,Ih,la,IS,IS,,M) Appliquant le Principe du Maximum d'Entropie, cette distribution peut typiquement se factoriser en : p(a'h,e,es,esa'ES'ESa o,Ih,Ia,IS,ISa,M) = p(O lo,M) p(h,ES Ih,IS,M) p(a'ESa Ia'ISa'M) p(eS'OSa IS'ISa'M) Un tel procédé comporte ainsi une étape 52a consistant à assigner la distribution de probabilité a priori p(O Io,M). A titre d'exemples, on peut assigner une 20 distribution a priori non-informative. Par exemple, si l'information I0 consiste seulement à savoir que BF et z appartiennent aux intervalles respectifs [BFmin,BFinax] et [zmin,zmaxi, alors la distribution de probabilité a priori p(010,M) est une distribution a priori uniforme 25 de Bayes-Laplace : p(Olo,M)=p(BF,z ID,M)= x[BFmin,BFinax ] (BF)x[imin ,imax ] (z )/(BFinax - BFmin )/(zmax -zmin A l'autre extrémité, on peut également assigner des distributions de probabilité informatives telles que les distributions des fréquences relatives empiriques 30 ou les distributions de probabilité marginales a15 posteriori obtenues à partir d'expériences antérieures, par exemple à partir de techniques d'imagerie de perfusion quantitatives telles que la Tomographie par Emission de Positrons (Positron Emission Tomography, PET suivant une terminologie anglo-saxonne) ou l'ASL (Arterial Spin Labelling). Un procédé conforme à l'invention comporte en outre une étape de configuration 52b consistant à assigner la distribution de probabilité a priori p(h,Es Ih,Is,M)= p(h Es,Ih,M)- p(ES Is,M) à partir de l'information Ih (ou Ig si on s'intéresse seulement à la fonction de répartition complémentaire). Cette information est constituée d'éléments d'information 15 dure et douce. Par exemple, en ce qui concerne hW, les éléments d'informations quantitatives dures N h(ti)=h(0)=0 et h(t)dt=ot-1h(ti)=1 suivant la méthode p i=2 des rectangles à droite peuvent être considérés. Ces deux informations permettent de réduire à 11T-2 le 20 nombre des valeurs/paramètres à estimer. Comme nous allons le voir, un choix judicieux peut consister à conserver les valeurs h'=[h(t2),...,h(tN_1)IT et à exprimer N-1 h(tN) par h(tN)=1/Ot- h(ti) . D'autre part, nous pouvons i=2 également considérer l'information quantitative dure 25 Vt>-0, hW 0. Il reste donc à assigner la distribution de probabilité conjointe a priori du vecteur h en combinant ces informations quantitatives dures avec une information purement qualitative et douce. Considérons par exemple l'information qualitative douce Ih telle 30 que « hW est plus ou moins lisse ». Cette information qualitative Ih peut se traduire d'abord en l'information quantitative douce « le carré de la dérivée seconde de h(t) est plus ou moins grand en moyenne ». Après discrétisation temporelle aux instants de mesure tl=1,N, une approximation numérique d'ordre 2 de la dérivée seconde de h est donnée par exemple par : d2h(tt) _ h(tt-1)+h(tt+l)-2h(tt)+O(Ot2 bi = 2,N-1, - dt2 Ot2 On peut également utiliser des approximations 10 numériques d'ordre supérieur, par exemple la formule d'ordre 4 : bi =3,N-2, d2h(ti) -12h(tt-2)+3h(tt-1)-2h(tz)+3h(tz+i)-12h(tt+2) +o(ot4) dt2 Ot2 Ces approximations numériques peuvent s'écrire sous dZh' une forme matricielle = Dh où D est par exemple la dt2 /0 0 0 0 0~ 1 -2 1 0 . 0 1 -2 1 0 15 matrice carrée D= 1 0 de At~ . 0 0 1 -2 1 0 0 1 -2 1 ~0 0 0 0 0) dimension N dans le cas de l'approximation d'ordre 2 et h' est tel que défini précédemment. Le carré de la norme euclidienne de la dérivée seconde de h' peut s'écrire alors 20 d2h = (Dh)T (Dh) = hT (DTD)h = hTWlh . 2 dt2 Supposer qualitativement la réponse impulsionnelle hW plus ou moins lisse peut donc se traduire de d'h' _ ,jhTWh est plus ou moins grande. On cherche donc la distribution 5 de probabilité a priori p(h' h(t1),h(tN),I) - les points extrêmes /40 et h(G) devant être traités séparément - parmi toutes les distributions de probabilité continues de même norme euclidienne d 2h' Principe du Maximum d'Entropie - qui consiste à choisir 10 parmi toutes ces distributions à support l'hyper- quadrant [0,+o1N2 celle ayant la plus grande entropie de Shannon difJférentielle (sous la mesure de référence de Lebesgue) car c'est la plus incertaine et donc la plus honnête on obtient la distribution 15 conditionnelle gaussienne multivariée tronquée sur [0,+00r 2 d'espérance mathématique vectorielle constante X12 W1-1 62 (ou la gaussienne tronquée sur [0,1]N z pour le vecteur des valeurs de la fonction de répartition 20 complémentaire) s12(h-M)T W (h-M) 26Sz où Cl(pi,el,6s) est la constante de normalisation -1 f exp -E12(h-M)TW1(h-M) dh' 2o-S2 h'E[O,+OEJN-2 Il apparaît donc deux nouveaux hyperparamètres dans 25 notre modèle global de perfusion, l'espérance manière quantitative par supposer que dtz dt2 . En suivant alors le M=(III ,...,)T et de matrice variance-covariance p(h'h(ti),h(tN),~l,~l,Es,Ih,M)=Cl(µ1,,6s)exp mathématique scalaire µl et la variance fractionnaire inverse el qui joue le rôle de paramètre de régularisation pour notre modèle de perfusion. El quantifie le poids de l'information a priori qualitative douce Ih par rapport au poids de l'information quantitatives procurée par les données expérimentales s et sa. D'autre part, par définition : p(h(tl),h(tN) Es,Ih,M) =p(h(tl) Es4,M) p(h( )I ,Ih1,M) = N-1 .5E1401.5 h(tN)-1/At~ h(ti) i=2 Cependant, afin de tenir compte du fait que le véritable délai z n'est presque jamais égal à un instant de mesure tl,i=1,N, on pourra poser p h(t1) µ1,1~E1,1,Es,I1h,M = C11(µ11,E11,6s)ex E1,12[h(tl)-µ1,1] 2 26s2 h(tl)E[0,11 où Ci,i (111,i' Ei,i' 6s est la constante de normalisation C1 1 01.11 1, E1 1, 6s ~ _ f exp h(t1)E[0,1] E1,12 [h(t1)-Nl,lj 215s2 dh (t1) de sorte à indiquer et seulement indiquer que h(z)=0+. Au final la distribution de probabilité a priori de h à support l'hyper-quadrant [0,+or pour l'information Ih peut s'écrire p(h µ1,E1,Es,Ih,M) = N-1 h(tN)-1/At+ E h(ti) i=2 Cl(µ1,s1,6s)exp -E12(h-M)TW1(h-M) 26s2 Ô [h(t1)1.5 De manière générale, nous noterons Eh (ou Eh) le vecteur des paramètres de la distribution a priori pour h, par exemple Eh -(µ1,e1) ou Eh -(µ1,E1,E1,1) Il reste encore à assigner une distribution de 5 probabilité a priori pour Eh. On peut typiquement assigner la distribution non-informative impropre de Bayes-Laplace-Lhoste-Jeffreys ( ll a) -1 p(µ1,£1,£11 Ih,M)0CE1-IE11-1 sur I". I". I". De la même manière, on peut introduire 10 l'information quantitative douce Ih : « le carré de la dérivée première de hW est plus ou moins grand en moyenne ». On peut ainsi obtenir au final une nouvelle distribution de probabilité a priori gaussienne 15 multivariée tronquée pour h par application du Principe du Maximum d'Entropie. Il convient de remarquer que les informations douces Ih et Ih telles que décrites précédemment sont 20 les seules informations, les seule contraintes qu'il est pleinement légitime d'introduire dans notre problème : ce ne sont pas des hypothèses arbitraires qui peuvent être vérifiées ou pas par l'expérience mais, au contraire, il s'agit juste de l'expression de 25 propriétés physiologiques fondamentales sans lesquelles le problème de l'estimation des paramètres hémodynamiques, des réponses impulsionnelles ou des fonctions de répartition complémentaires n'aurait en fait aucun sens : elles sont logiquement nécessaires et 30 suffisantes à notre problème. Toute autre information serait à l'inverse une simple hypothèse pouvant être éventuellement vérifiée ou falsifiée par l'expérience. Conformément à introduire d'autres on peut indiquer et impulsionnelle suit fonctionnelle sans l'aide d'une l'invention, on peut toutefois informations douces. Par exemple, seulement indiquer que la réponse plus ou moins une certaine forme pour autant la contraindre - à dure - à suivre exactement information cette forme. L'ajout de ce type d'informations semi- quantitatives douces permet de déterminer dans quelle mesure les réponses impulsionnelles peuvent être décrites par les formes fonctionnelles proposées. Supposons donc par exemple qu'on veuille introduire l'information douce Ih « h(t) suit plus ou moins un modèle paramétrique ou semi-paramétrique Mh:f(t,Oh) », 15 Oh étant le vecteur des paramètres dudit modèle. Comme dit précédemment, de tels modèles paramétriques ont été proposés, par exemple le modèle F à deux paramètres donné par : MTT f (t,Oh) = t P e-t/P =(MTT,) ( ) : fonction Gamma d'Euler MTT /P /r MTT>0,P>0 MTT 20 Notons que dans estimé directement, le premier moment l'intégrale dece cas, le paramètre MTT peut être sans avoir à estimer numériquement de la réponse impulsionnelle (ou la fonction de répartition complémentaire) comme décrit précédemment. 25 Indiquer que h suit plus ou moins une certaine forme fonctionnelle f(t,Oh) revient alors à indiquer quantitativement que la norme euclidienne des résidus 2 N h-f(t,Oh) =~[h(tl)-f(tl,Oh)j2 est plus ou moins grande. i=1 En appliquant à nouveau le Principe du Maximum d'Entropie comme précédemment, on trouve de la même manière que la distribution de probabilité a priori de h est la gaussienne multivariée tronquée sur l'hyper- quadrant [0,+4 3 E32 N p(h Eh,ES,Oh,Ih,Mh,M)=C3(s3,6S,Oh)exp - 21E11(0-f (ti,eh)12 26S i=1 OÙ C3(s3,6S,Oh) est la constante de normalisation f exp - s3 2 h - f (t, eh hE[0,+oo]N 26S C3(s3,6S,Oh)= De manière équivalente, on pourrait introduire une information douce Ih telle que « hW suit plus ou moins un vecteur de valeurs donné h=[h(tl),...,h(tN)iT » Par application du Principe du Maximum d'Entropie, on obtient alors la distribution de probabilité a 15 priori : _ 2 N p(h Eh,Es,Ih,M) oc C4 (s4,6s,h)exp - &4 2 1[11(0- h(t.)12 26s i=1 où C4(s4,6,h) est la constante de normalisation: --1 hE[0,+oo]N L'invention permet également de pouvoir combiner 20 plusieurs informations douces sur la réponse impulsionnelle hW (ou la fonction de répartition complémentaire) et leurs distributions de probabilité a Priori correspondantes. Ainsi, si pl h Eh,Ih,M I,..., pl h Eh,Ih,M I sont n distributions de 25 probabilité a priori traduisant des informations C4(s4,6s,h)= f 2 _ exp - s4 h - h 2 d 26s2 Ih,...,Ih, de paramètres de régularisation respectifs EL..,Eh (avec par exemple Eh =(el,µl), Eh =(E3,Oh), etc.), alors une distribution a priori pour h tenant compte de c e s n informations p e u t s' é c r i r e p(h Eh,Ih,M)=n p(h Eh,Ih,M) en notant Eh =(Eh,...,En) e t k=1 Ih = Ih A ... A Ih Pour encoder une information la sur la fonction d'entrée artérielle locale, un procédé conforme à l'invention comporte en outre une étape 52c consistant à assigner une distribution de probabilité a priori p(a,Esa Ia,Isa,M)=p(aEsa,Ia,M).p(Esa Isa,M). Une telle distribution est assignée de la même manière que celle relative à la réponse impulsionnelle h, e n introduisant et en combinant des informations dures et/ou douces sur la fonction d'entrée artérielle. Par exemple, on peut spécifier que la fonction d'entrée artérielle est plus ou moins lisse, positive, unimodale, bimodale, nulle à l'origine, asymptotiquement nulle, d'aire donnée ou encore indiquer et seulement indiquer qu'elle suit plus ou moins un certain modèles paramétrique ou semiparamétrique Ça (t,Oa où Oa e s t un v e c t e u r de paramètres, comme par exemple le modèle « tri-Gamma » à 25 onze paramètres: Ca(t,0a)= a0-1 -(t-to)lf3o b(t tl ) ai-1 -(t tl)/~3l 1-a -b)(t- a2-1 -(t-t2)/~32 -t0) e - e (t2) e POa°T (ao) + PlaIF(al) + P2a2E(a2) ea =(a,b,a0,P0,t0,a1,P1,t1,a2,P2,t2) T ( ) : fonction Gamma d'Euler Comme décrit précédemment, on obtient une distribution de probabilité a priori gaussienne 2 (e â) N p(aEâ,Esa,Oa,I1a,Ma,M)=C(Eâ,6sa,0a)exp~- 21[a(ti)-Ca(ti,ea) 2o-sa i=1 où C(Ea,6sa,0a) est la constante de normalisation : -1 De même, il est possible d'indiquer et de seulement indiquer (sans quoi on se ramène au cas où la fonction d'entrée artérielle est donnée) que le vecteur a suit 10 plus ou moins un vecteur de valeurs donné T cc =~cc(tl),...,cc(tN)~ tel une fonction d'entrée artérielle locale on obtient la distribution de probabilité a priori (Ela)2 a-Ca(t,ea) II 2 2o-sa 2 C(Ea,6sa ,0a)= f exp aE[O,+oeiN p(a Eâ,Esa ,Ia,M)=Ca (E2 a,o-sa,a)ex(N 2 2 ~[a(ti)-a(ti)1 26Sa i=1 15 où ( ,o-a,-a) est la constante de normalisation 2 exp -(En a-u 2 f 2o-sa 2 aE[O,+coN e --1 Ça (Ea2, o-sa , . = Un procédé conforme à l'invention comporte en outre une étape de configuration 52d consistant à assigner une distribution de probabilité a priori p(Es,Esa,Os,Osa Is,Isa,M) qui peut typiquement se factoriser en : / p(Es,Esa,es,esa IS'ISa,M) -p(ES IS,M)'p(ESa ISa,M)'p(es IS,M)'p(esa ISa,M) On a par exemple la distribution de probabilité a priori non-informative . p(ES,ESa,OS,OSa IS,ISa ,M) -P(os IS,M)'P(6Sa ISa,M)'P(SO IS,M)'P(SOa ISa,M) (us.us ) -1 (S0max S0min) (SOmax SOmin ) a a a En tenant compte des différents hyperparamètres introduits dans les étapes 52b à 52d, la distribution de probabilité conjointe a priori 53 peut se réécrire comme p(a,Ea,h,Eh,O,Os,Osa,Es,Esa I0,Ih,Ia,Is,Isa,M) . Pour simplifier les expressions suivantes, nous notons I=(IO,Ih,Ia,IS,ISa,M) l'ensemble des informations introduites comme autant de paramètres de configuration de l'unité de traitement. Etant données une distribution de probabilité directe 54 et une distribution de probabilité conjointe a priori 53 telles que décrites précédemment, on obtient la distribution de probabilité conjointe a posteriori 55 de tous les paramètres en appliquant la règle de Bayes : p(a,Ea,h,Eh,O,Os,Osa,Es,Esa s,sa,I)= p(a,Ea,h,Eh,0,0s,esa,Es,Esa I).p(s,sa a,h,O,Os,Osa,Es,Esa,I) p(s,sa I) p(a,Ea,h,Eh,0,0s,esa,Es,Esa I).p(s,sa a,h,O,Os,Osa,Es,Esa,I) La configuration préalable étant réalisée, l'invention permet à présent d'estimer une quantité d'intérêt que nous noterons 6 parmi tous les éléments du vecteur E=(a,Ea,h,Eh,O,Os,esa,Es,Esa). Par exemple, 6 = BF ou 0=1, ou bien encore 6 = h(O , etc. Un procédé conforme à l'invention comporte donc une 10 étape 56 consistant à évaluer la distribution marginale a posteriori pour 6, s o i t p(9 s,sa,I)= J p(E s,sa,I)d(E\9) . E\0 A partir de cette distribution marginale a posteriori, on peut obtenir 57 des estimations 9 de O. 15 Par exemple, on obtient l'estimateur de Bayes sous 2 fonction de coût quadratique L 0-0 =0-0 Q où est la norme euclidienne en prenant l'espérance mathématique de cette distribution 9Q =J 9 - p(9 s,sa,I)d9 . 0 De même, on peut obtenir l'estimateur du maximum a -P 20 posteriori 9 (Maximum a posteriori Estimator, MAP suivant une terminologie anglo-saxonne) par 0 =argmaxp(9 s,sa,I) . e A titre d'exemple, on peut obtenir la distribution 25 de probabilité marginale a posteriori de la valeur oc h(ti),i=1,N de la réponse impulsionnelle en chaque instant de mesure ti, en marginalisant tous les autres instants . bi=1,N, p(h(ti)s,sa,I)= f p(hs,sa,I)dh(tl)...dh(ti#i)...dh(tN) h(tj ) ji et par suite obtenir des estimations de ces valeurs telles que hQ(ti) ou hP(ti) . De la même manière, il est également possible d'obtenir 57 une estimation h(x) (ou R(x)) de la valeur de la réponse impulsionnelle h(x) à un instant x quelconque, non forcément confondu avec un instant de mesure y Il suffit pour cela d'introduire h(x) dans l'expression des valeurs de la fonction de répartition complémentaire R(t) comme décrit précédemment et de calculer sa distribution marginale de probabilité. Il est même souhaitable d'introduire de tels instants x1,...,xL dans le problème d'estimation puisque les approximations numériques des intégrales telles que t +oc) t f «1-)dz f «1-)dz =1 f Ça (z)- o o o t-i H(t-r)- f h(v)dv o d'c ainsi que l'approximation numérique de dtzh ou dh seront d'autant meilleures que le nombre d'instants pris en compte sera grand. Les estimations résultantes seront donc également d'autant meilleures. Ensuite, puisque MTT=Eh= I t-h(t)dt, on peut obtenir des estimations de ce paramètre telles que l'estimation N « en moyenne » MTTQ =Ot-ltl -hQ(tz) ou l'estimation la i=2 P N hP plus probable MTT=Ot (ti) en appliquant une i=2 méthode d'intégration numérique, par exemple ici, la méthode des rectangles à droite. En outre, un procédé selon l'invention peut comporter une étape 58 pour obtenir une estimation de la précision sur l'estimation du paramètre e, voire des intervalles de confiance pour ces estimations. L'invention prévoit qu'un tel procédé puisse comporter en outre une étape 59 pour obtenir des paris sur lesdits intervalles de confiance. Par exemple, la précision sur l'estimation eQ peut être quantifiée par la matrice variance-covariance de la distribution marginale a posteriori pour e : 2 T = f /e-êQyep(e s,sa,I)de 0 On a alors par exemple un (hyper-)intervalle de confiance à « un sigma » p o u r e J= e-Q - dag (o--0-Q), eQ + dag (o--0-Q) et la probabilité que e appartienne à cet intervalle pJ = f p(e s,sa,I)de ou, de 0EJ manière équivalente, le taux de pari p' (odds 1-pJ suivant une terminologie anglo-saxonne). Il est également possible d'obtenir des estimations, des intervalles de confiance et des paris sur ces intervalles de confiance pour le paramètre BV =BF-MTT , le vecteur des valeurs de la fonction de répartition complémentaire b, le vecteur des valeurs de la fonction de sortie veineuse v= Ah ou encore le vecteur c=BF.Ab des valeurs de la courbe de concentration théorique car, connaissant la distribution de probabilité conjointe de plusieurs variables aléatoires, on peut calculer la fonction densité de probabilité d'une fonction quelconque en dépendant. Par exemple, par linéarité de l'espérance mathématique, on obtient immédiatement les estimations des valeurs de la fonction de répartition complémentaire R(t) à partir de celles de la réponse impulsionnelle i bQ = 1,1-Ot.h..,1-Ot. i=2 N bP= 1,1-Ot.hP(t2),...,1-Ot.EhP(tN) i=2 en suivant par exemple la méthode des rectangles à droite. De la même manière, de la distribution de probabilité conjointe p(c,Os s,sa,I) de c et Os, on peut déduire également celle du vecteur sth-LSth(t1 ,"Sth(tNg ' ~ des valeurs du signal de perfusion théorique Sth(t) puisque sth =Y(c,Os) . A partir de cette distribution, on peut enfin obtenir 57 des estimations stil ainsi que des hyper-intervalles de confiance 58 et des paris 59 sur ces hyper-intervalles suivant une méthode similaire à celle décrite précédemment. L'invention prévoit que les intervalles de confiances ou paris sur lesdits intervalles de confiance, puisse permettre d'ajuster 62 la configuration de l'unité de traitement. Ainsi, il est possible, de modifier les paramètres de configuration et proposer des estimations de meilleure qualité. L'invention prévoit en outre, qu'un procédé puisse comporter une étape 60 pour calculer des résidus r(ti)=S(tl)-sth(ti),i=1,N entre les signaux par imagerie de perfusion théoriques et expérimentaux. L'invention 30 permet alors de calculer diverses statistiques ou distances D(s,sth) entre ces vecteurs, la plus classique étant la somme des carrés des résidus ou khi carré N z2 ,D-(02 (Sum of Squared Errors, S SE selon une i=1 terminologie anglo-saxonne). Ces diverses statistiques permettent de quantifier l'adéquation du modèle de perfusion aux données expérimentales. On obtient ainsi des « cartes d'erreur » pour le modèle pour chaque voxel d'intérêt. Pour les raisons évoquées précédemment (i.e. problème de sur-ajustement), la quantification de l'adéquation du modèle de perfusion aux données expérimentales s et sa pourra se faire de manière particulièrement avantageuse en calculant la probabilité des données expérimentales (s,sa) sachant le modèle de perfusion en chaque voxel d'intérêt : p(s,sa I)=J p(s,sa E,I)p(8I)dE Dans ce cas, la carte d'erreur serait basée sur 1-p(s,sa I) . L'invention prévoit également que l'on puisse appliquer 61 divers tests statistiques ou diverses techniques de diagnostic graphiques tels que les courbes de Henri (Q-Q plot suivant une terminologie anglo-saxonne) ou les cartes de Poincaré (return maps suivant une terminologie anglo-saxonne) afin de vérifier si les résidus «0 sont bien indépendants, identiquement distribués et gaussiens, etc. L'invention prévoit ainsi, par un processus itératif et essais successifs, de pouvoir corriger et affiner 62 le processus 50 de configuration, notamment de corriger et affiner des modèles de perfusion théoriques afin de progresser dans la modélisation, la compréhension et le traitement des phénomènes de perfusion et d'obtenir in fine de meilleures estimations des paramètres hémodynamiques, des réponses impulsionnelles, des fonctions de répartition complémentaires, des fonctions d'entrée artérielles ou encore des fonctions d'entrée veineuses. En liaison avec la figure 8, nous nous proposons de décrire à présent un procédé conforme à l'invention selon un deuxième exemple d'application pour lequel les fonctions d'entrée artérielles théoriques sont supposées être données avec une certitude absolue et une précision infinie à un délai temporel z près. Un modèle de perfusion expérimental M peut s'écrire alors comme : =11(c,os)+ = BF - Ab(t-r) = BF -B(t-i)a où toutes les quantités sont telles que définies précédemment sauf a =LCa (q)," ,Ca (tN)jT qui e s t le vecteur des valeurs de la fonction d'entrée artérielle théorique désormais supposé connu. 20 La distribution de probabilité conjointe directe s'écrit alors p(s a,h,O,Os,Es,I) et la distribution de probabilité conjointe a priori p(h,Eh,O,Os,Es I) avec I =(Io,Ih,IS,M) . Ces distributions s'assignent comme décrit précédemment. La règle de Bayes devient p(h,Eh,Q,QS,ES s,a,I) x 25 p(h,Eh,O,OS,ES a,I)' p(s a,h,O,Os,Es,I) On obtient ensuite des estimations, des intervalles de confiance et des paris sur ces intervalles de confiance pour tout paramètre 6EE=(h,Eh,O,Os,Es) de la même manière que celle décrite précédemment. 30 En liaison avec la figure 8, nous nous proposons de décrire à présent un procédé conforme à l'invention selon un troisième exemple d'application pour lequel les fonctions d'entrée artérielles théoriques ne sont pas données et où les signaux d'entrée artériels ne sont pas mesurés. Un modèle de perfusion expérimental M peut s'écrire alors comme : =11(c,os)+ =BF-Ab(t)=BF-B(t)a où toutes les quantités sont telles que définies précédemment. La distribution de probabilité conjointe directe 15 s'écrit toujours p(s a,h,O,Os,Es,I) et la distribution de probabilité conjointe a priori s'écrit désormais p(a,Ea,h,Eh,O,OS,ES I) avec I =(Ip,Ih,Ia,IS,M) . La règle de Bayes devient p(a,Ea,h,Eh,O,OS,ES s,I) oc p(a,Ea,h,Eh,O,Os,Es I)-p(s a,h,O,Os,Es,I) On obtient ensuite des estimations, des intervalles de confiance et des paris sur ces intervalles de confiance pour tout paramètre 6EE=(a,Ea,h,Eh,O,OS,ES) de la même manière que celle décrite précédemment. 25 Il convient de remarquer, quelque soit l'exemple d'application retenu, qu'à l'inverse des méthodes connues, un procédé d'estimation conforme à l'invention est une méthode exacte, dans le sens où elle consiste 30 et consiste seulement à traduire des informations qualitatives, quantitatives ou semi-quantitatives dont 20 on dispose a priori sur les quantités d'intérêt dans la Théorie des Probabilités Bayésienne afin de déterminer univoquement l'information a posteriori sur ces mêmes quantités, procurée par le biais des mesures expérimentales. Aucune hypothèse arbitraire qui pourrait ne pas être vérifiée - en particulier sur les réponses impulsionnelles ou les fonctions de répartition complémentaires - n'est nécessaire car l'invention introduit seulement les distributions de probabilités les plus incertaines codifiant les différentes contraintes qualitatives douces logiquement nécessaires et suffisantes pour résoudre le problème. Il s'en suit que, dans les cas où les fonctions d'entrée artérielles ne sont pas supposées données mais tout au plus mesurées, ces méthodes permettent de tester si les signaux d'entrée artériels proposés peuvent effectivement correspondre aux « véritables » fonctions d'entrée artérielles locales pour chaque voxel considéré : en effet, si les signaux d'entrée artérielles ne sont pas adéquats et ne correspondent pas aux véritables fonctions d'entrée artérielles locales, il se peut qu'il n'y ait pas de jeu de paramètres (paramètres hémodynamiques, fonction de répartition complémentaire, etc.) qui soit solution du problème tout en respectant les différentes contraintes posées a priori (ou du moins dont la probabilité n'est pas négligeable a priori). Dans ce cas, la Théorie des Probabilités interprétera les signaux d'entrée artériels proposés comme du bruit : les écarts-types 6a/E, seront beaucoup plus grands par rapport à ceux typiquement obtenues à partir de signaux d'entrée artériels plus adéquats. L' invention trouve ainsi une application particulièrement intéressante pour tester les différentes méthodes de sélection ou d'estimation de fonctions d'entrée artérielles globale ou locales suivant l'état de l'art. S'il s'avère que les estimations, en particulier celles sur leurs paramètres dé régularisation Ea, obtenues à partir de ces fonctions d'entrée artérielles globale ou locales sont trop souvent aberrantes d'un voxel à un autre, on pourra conclure qu'il faut soit introduire de nouvelles méthodes (locales) de sélection plus appropriées pour ces fonctions soit recourir à des méthodes ne nécessitant pas la donnée de fonctions d'entrée artérielles ou la mesure préalables de signaux d'entrée artériels, ce qui est précisément l'objet de la troisième méthode décrite précédemment. A titre d'exemple d'application, nous pouvons citer les principales étapes de mise en oeuvre de l'invention au moyen d'un système d'analyse d'imagerie de perfusion adapté, tel que celui décrit en figures 1 ou 2 : - ouverture d'un dossier patient ou prise en compte de séquences d'images par l'unité de traitement 4 (ou de prétraitement 7) pour sélectionner des séquences d'images d'intérêt - en particulier, sélection des images I1 à In de perfusion au cours du temps à partir desquelles sont obtenus les signaux de perfusion S(t) pour chaque voxel, tel qu'illustré en figure 5a ; - prévisualisation au moyen d'une interface homme- machine 5 des images pour permettre à un utilisateur 6 d'identifier des tranches ou des zones d'intérêt ; - configuration de l'unité de traitement 4 à partir des paramètres de configuration (informations introduite) pour permettre la mise en oeuvre du procédé d'estimation conforme à l'invention ; - choix de la (ou les) quantité(s) d'intérêt à estimer ; estimation par l'unité de traitement 4 de paramètres hémodynamiques 14, tels que BF, z ou MUT, pour un organe tel que le cerveau humain ; - estimation optionnelle d'autres paramètres tels que Eh, Es, Os ou encore Ea, Esa ou Osa lorsque la fonction d'entrée artérielle n'est pas donnée ; - estimation optionnelle des vecteurs des valeurs des réponses impulsionnelles h, des fonctions de répartition complémentaires b des vecteurs des valeurs des fonctions d'entrée artérielles a, des vecteurs des valeurs des fonctions de sortie veineuses v, des vecteurs des valeurs des concentrations théoriques c, des signaux théoriques stil ou encore des résidus r; - délivrance desdites quantités d'intérêt estimées 14 à l'interface homme-machine 5 pour que celle-ci les présente in fine par exemple sous la forme de cartes où l'intensité ou la couleur de chaque pixel dépend de la valeur calculée pour en restituer la teneur au praticien ; - affichage optionnel des estimations des réponses impulsionnelles, des fonctions de répartition complémentaires, des fonctions d'entrée artérielles, des fonctions de sortie veineuses, des concentrations théoriques, des signaux théoriques ou encore des résidus pour certains voxels d'intérêt choisis par l'utilisateur ; - affichage optionnel de cartes de confiance ou de cartes de paris pour un ou plusieurs paramètres d'intérêt tels que les paramètres hémodynamiques ; - affichage optionnel d'intervalles de confiance ou de paris sur ces intervalles de confiance pour certaines réponses impulsionnelles, fonctions de répartition complémentaires, fonctions d'entrée artérielles, etc. pour certains voxels d'intérêt choisis par l'utilisateur ; - affichage optionnel de cartes d'erreur pour une ou plusieurs distances entre données expérimentales et modèle de perfusion non-paramétrique global, en particulier affichage de la probabilité des données expérimentales sachant le modèle de perfusion global ; - sélection assistée de ladite zone pathologique d'intérêt par le praticien, caractérisée par une anomalie de la distribution d'un ou de plusieurs paramètres hémodynamiques, des réponses impulsionnelles (ou des fonctions de répartition complémentaires) ou des fonctions d'entrée artérielles locales; - estimation, par l'unité de traitement, du volume de la zone de tissus anormalement perfusée éventuellement connexe à une z o n e l é s é e et p o u r laquelle le praticien pourra décider une action thérapeutique (thrombolyse intraveineuse afin de résorber le caillot sanguin par exemple) ; - estimation, par l'unité de traitement, de certaines quantités, telles que le ratio des volumes des zones lésée et anormalement perfusée sur lesquelles, un praticien pourra peaufiner son diagnostic et sa prise de décision thérapeutique (thrombolyse intraveineuse afin de résorber un caillot sanguin par exemple). Comme décrit précédemment, le processus de configuration 50 de l'unité de traitement 4 peut être réalisé par l'unité elle-même (mise en oeuvre du processus 50). En variante, ladite configuration peut consister à mémoriser et sélectionner une bibliothèque de distributions de probabilité conjointes a posteriori selon les quantités d'intérêt que l'on souhaite estimer. L'élaboration de cette bibliothèque peut être réalisée à l'aide d'une unité dédiée et apte à coopérer avec l'unité de traitement 4. En variante, il peut exister des itérations à la suite de l'estimation d'intervalles de confiance, de paris sur ces intervalles de confiance pour certaines quantités d'intérêts afin d'affiner ladite configuration. L'élaboration de distances entre données expérimentales et modèle de perfusion non-paramétrique global, en particulier affichage de la probabilité des données expérimentales sachant le modèle de perfusion global peut également induire une mise à jour de la configuration. L'invention prévoit donc d'afficher des estimations des paramètres sous forme de « cartes de paramètres » où l'intensité ou la couleur de chaque voxel dépend de la valeur calculée, par exemple de manière linéaire. De même, l'invention prévoit éventuellement d'afficher les écarts-types de ces estimations sous la forme de « cartes de confiance » ainsi que les paris sur les intervalles de confiance correspondant sous la forme de « cartes de paris ». Pour ce qui est des estimations des vecteurs de valeurs des réponses impulsionnelles, des fonctions de répartition complémentaires, des fonctions d'entrée artérielles, des fonctions de sortie veineuses, des courbes de concentration, des signaux de perfusion ou encore des résidus, l'invention prévoit leur affichage sous forme de séries temporelles pour chaque voxel où l'utilisateur le demande. Enfin, l'invention prévoit l'affichage de distances entre signaux expérimentaux et modèle de perfusion non-paramétrique ou de la probabilité des données expérimentales sachant ce modèle sous la forme de « cartes d'erreur ». Les figures 9 à 12 permettent d'illustrer un mode d'affichage sous la forme de cartes, de certaines quantités d'intérêts tels que des paramètres hémodynamiques 14 estimés conformément à l'invention voire des écarts-types ou probabilités qui leur sont associés. Ainsi pour un cerveau humain analysé à l'aide d'imagerie par Résonance Magnétique Nucléaire, la figure 9 permet de visualiser une estimation des volumes sanguins cérébraux (cerebral blood volumes, CBV suivant une terminologie anglo-saxonne). Une telle carte (458 x 458 pixels) permet de mettre en évidence une zone ischémique probable 80. En effet, il est possible de constater à l'aide d'une interface 6 adaptée, une nette augmentation du paramètre CBV dans le territoire de l'artère cérébrale postérieure droite par rapport à l'hémisphère controlatérale. Une vasodilatation consécutive à l'ischémie peut être révélée par une lecture de la carte telle qu'illustrée en figure 9. La figure 10 permet d'illustrer une carte (458 x 458 pixels) relative à l'estimation des flux sanguins cérébraux en cas d'ischémie cérébrale. On peut constater en analysant la carte, une diminution du paramètre (cerebral blood volumes, CBF suivant une terminologie anglo-saxonne) dans le territoire de l'artère cérébrale postérieure droite par rapport à l'hémisphère controlatérale consécutive à l'ischémie. Une telle carte permet de mettre en évidence une zone ischémique probable 80. La figure 11 permet d'illustrer une carte (458 x 458 pixels) relative à l'estimation des temps de transit moyen MET. On peut constater en analysant la carte, une nette augmentation des MU dans le territoire 81 de l'artère cérébrale postérieure droite par rapport à l'hémisphère controlatérale consécutive à l'ischémie. La figure 12 permet de décrire une carte (458 x 458 pixels) relative à l'estimation de la probabilité que le flux sanguin cérébral soit dans l'intervalle de confiance [CBF-o-cBF,CBF+o-cBF]. On peut constater en analysant la carte que, mis à part des ajustements aberrants, les probabilités sont centrées autour de 0,68. C'est une valeur à laquelle on est en droit de s'attendre si la distribution de probabilité a posteriori de CBF suit une loi normale. Grâce aux cartes présentées précédemment, l'invention permet de mettre à la disposition d'un utilisateur tout un ensemble d'informations utiles, informations qui ne pouvaient être disponibles à l'aide des techniques connues de l'état de l'art. Cette mise à disposition est rendue possible par une adaptation de l'unité de traitement 4 selon les figures 1 ou 2 en ce que ses moyens pour communiquer avec le monde extérieur de l'unité 4 sont aptes à délivrer les paramètres estimés 14 selon un format approprié à une interface homme-machine 5 apte à restituer à un utilisateur 6 lesdits paramètres estimés sous la forme par exemple de cartes telles qu'illustrées par les figures 9 à 13. Grâce à l'invention, les informations délivrées sont ainsi plus nombreuses et justes. Les informations dont dispose le praticien sont ainsi de nature à accroitre la confiance du praticien dans sa détermination d'un diagnostic et sa prise de décision thérapeutique. Pour améliorer les performances du système selon l'invention, celle-ci prévoit que l'unité de traitement puisse être dotée de moyens pour paralléliser des calculs sur les voxels de l'image pour lesquels l'estimation des paramètres hémodynamiques, des fonctions de répartition complémentaires ou des fonctions d'entrée artérielles est requise. Cela peut être effectué en utilisant des technologies matérielles telles que les microprocesseurs graphiques (Graphical Processor Unit ou GPU, selon une terminologie anglo-saxonne) ou les grappes de calcul (clusters) ou logicielles telles que les Méthode de Monte-Carlo parallèles, etc. En variante, l'unité de traitement conforme à l'invention peut s'appuyer sur des moyens de calcul distants. Les temps de calculs peuvent ainsi être encore considérablement réduits. REVENDICATIONS1. Procédé pour estimer (57) une quantité d'intérêt (14) parmi une pluralité 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 (4) d'un système d'analyse d'imagerie de perfusion, et comportant une étape pour estimer ladite quantité d'intérêt à partir d'une donnée (15) de perfusion expérimentale, caractérisé en ce que ladite étape consiste à évaluer - suivant la méthode de Bayes - une distribution marginale a posteriori (56) pour ladite quantité d'intérêt par : - l'assignation (54) d'une distribution de probabilité directe de la donnée de perfusion sachant les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ; - l'assignation (53) d'une distribution de probabilité conjointe a priori desdites quantités, par l'introduction (51b) d'une information douce relative à la réponse impulsionnelle dudit système dynamique. 2. Procédé pour estimer (57) une quantité d'intérêt (14) parmi une pluralité d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe, ledit système dynamique étant linéaire, invariant dans le temps et déterminé formellement par la relation C(t)=BF-CQ(t)®R(t) où C(t) est la concentrationd'un agent de contraste circulant dans un voxel, Ça (t) est la concentration dudit agent de contraste dans l'artère alimentant ledit voxel, BF est le flux sanguin dans ledit voxel, désigne le produit de convolution et R(t) est la fonction de répartition complémentaire du temps de transit dans ledit voxel, ledit procédé étant mis en oeuvre par une unité de traitement (4) d'un système d'analyse d'imagerie de perfusion, et comportant une étape pour estimer ladite quantité d'intérêt à partir d'une donnée (15) de perfusion expérimentale, caractérisé en ce que ladite étape consiste à évaluer (56), suivant la méthode de Bayes, une distribution marginale a posteriori pour ladite quantité d'intérêt par : l'assignation (54) d'une distribution de probabilité directe de la donnée de perfusion sachant les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ; l'assignation (53) d'une distribution de probabilité conjointe a priori desdites quantités, par l'introduction (51b) d'une information douce relative à la fonction de répartition complémentaire du temps de transit R M dans ledit voxel. 3. Procédé selon la revendication précédente, caractérisé en ce que l'assignation (53) d'une distribution de probabilité conjointe a priori desdites quantités est réalisée en outre par l'introduction (51c) d'une information doucerelative à la concentration dudit agent de contraste dans l'artère alimentant le voxel. 4. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il e s t mis en oeuvre par itérations successives pour une pluralité de voxels considérés. 5. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il comporte une étape pour calculer (58) une information complémentaire sous la forme d'un intervalle de confiance associé à une quantité d'intérêt estimée. 6. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il comporte une étape pour calculer (59) une information complémentaire sous la forme d'un taux de pari sur un intervalle de confiance associé à une quantité d'intérêt estimée. 7 Procédé selon la revendication 4, caractérisé en ce qu'il comporte une étape pour calculer (60) une information complémentaire sous la forme d'une mesure de l'adéquation du produit de : - l'assignation (54) de la distribution de probabilité directe de la donnée de perfusion sachant les paramètres intervenant dans le problème de l'estimation des quantités d'intérêt du système dynamique artère/tissu/veine du voxel considéré ;- l'assignation (53) de la distribution de probabilité conjointe a priori desdites quantités. 8. Procédé selon l'une quelconque des revendications 5 à 7, caractérisé en ce qu'il comporte une étape pour délivrer une quantité d'intérêt estimée (14) à une interface homme-machine (5) apte à la restituer à un utilisateur (6). 9. Procédé selon l'une quelconque des revendications 5 à 8, caractérisé en ce qu'il comporte une étape pour délivrer toute information complémentaire associée à ladite quantité d'intérêt estimée, à une interface homme-machine (5) apte à la restituer à un utilisateur (6). 10. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce q u e la donnée de perfusion expérimentale consiste en un vecteur de valeurs d'un signal de perfusion expérimental ou en la conversion de ce dernier en un vecteur de valeurs d'une courbe de concentration. 11. 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 de perfusion expérimentale (15);les moyens de traitement sont adaptés pour mettre en oeuvre un procédé pour estimer une quantité d'intérêt (14) parmi une pluralité d'un système dynamique artère/tissu/veine d'un volume élémentaire - dit voxel - d'un organe, selon l'une quelconques des revendications 1 à 10. 12. Unité de traitement selon la revendication précédente, caractérisé en ce que les moyens pour communiquer délivrent (57) une quantité d'intérêt estimée (14) selon un format approprié à une interface homme-machine (5) apte à la restituer à un utilisateur (6). 13. Unité de traitement selon la revendication précédente, caractérisé en que les moyens pour communiquer délivrent (58, 59) une information complémentaire associée à une quantité d'intérêt estimée (14) selon un format approprié à une interface homme-machine (5) apte à la restituer à un utilisateur (6). 14. Système d'analyse d'imagerie de perfusion comportant une unité de traitement (4) selon l'une quelconques des revendications 11 à 13 et une interface homme-machine (5) apte à restituer à un utilisateur (6) une quantité estimée (14) selon un procédé conforme à l'une quelconque des revendications 1 à 10 et mis en oeuvre par ladite unité de traitement (4).