Brevet FR2900253A1 — publié le 20071026

Procédé de reconstruction de la distribution des propriétés optiques d'un milieu diffusant inhomogène Domaine technique de l'invention L'invention concerne un procédé de reconstruction de la distribution de propriétés optiques pour examiner un objet diffusant inhomogène par éclairement de l'objet: avec une lumière continue, le procédé comportant la formulation d'une équation de diffusion comportant un nombre d'onde initial. État de la technique La tomographie de fluorescence consiste à déterminer la répartition tridimensionnelle de fluorophores dans un objet comportant un milieu diffusant. L'objet est illuminé par une lumière dont la longueur d'onde correspond à l'excitation des fluorophores. Ensuite, on mesure la lumière de fluorescence, dont la longueur d'onde est décalée par rapport à la longueur d'onde d'excitation. Généralement, le milieu diffusant présente des inhomogénéités. Afin de prendre en compte les inhomogénéités, deux approches sont typiquement utilisées : La première approche consiste à normaliser les signaux de mesure de la fluorescence par rapport à la lumière d'excitation, comme décrit dans l'article Experimental three-dimensional fluorescence reconstruction of diffuse media by use of e normalized Born approximation de V. Ntziachristos et al. (Optics Letters Vol.26, No.12, 2001). Cependant, cette correction n'est pas suffisante, car les propriétés optiques du milieu à la longueur d'onde d'excitation et à la longueur d'onde de fluorescence sont différentes et les trajets optiques suivis par les photons de fluorescence et d'excitation sont différents, Cette première approche utilise des signaux continus et une géométrie plane. La seconde approche consiste à déterminer, de façon laborieuse, la cartographie du coefficient d'absorption et celle du coefficient de diffusion, comme décrit dans l'article Fluorescence-Enhanced Optical Tomography Using Referenced Measurements of Heterogeneous Media de R. Roy et al. (IEEE Transactions on Medical Imaging Vol.22, No.7, 2003). Cette seconde approche utilise des signaux fréquentiels et/ou temporels et est, ainsi, plus complexe. Par ailleurs, on cherche à développer des procédés qui permettent, en plus de la détermination de la répartition tridimensionnelle de fluorophores, de déterminer toutes les caractéristiques optiques d'un objet inhomogène par une approche générale. Objet de l'invention L'invention a pour but un procédé de reconstruction de la distribution de 20 propriétés optiques pour examiner un objet diffusant inhomogène par éclairement de l'objet avec une lumière continue. Selon l'invention, ce but est atteint par le fait que, l'équation de diffusion comportant au moins un terme d'inhomogénéité tenant compte des 25 variations spatiales des propriétés optiques, le terme d'inhomogénéité est incorporé dans un nombre d'onde effectif. Selon un mode de réalisation préférentiel de l'invention, la lumière continue étant une lumière d'excitation de fluorophores répartis dans l'objet, le 30 procédé comporte let reconstruction d'une image de tomographie de fluorescence.15 Selon un développement de l'invention, le carré du nombre d'onde effectif est défini par la somme du terme d'inhomogénéité et du carré du nombre d'onde initial. Selon un mode de réalisation particulier de l'invention, le terme d'inhomogénéité correspond à l'expression oD(r) V , D, étant le coefficient DÂ(r) de diffusion dépendant de la coordonnée spatiale r . La répartition des fluorophores peut notamment être déterminée par l'intermédiaire d'un formalisme utilisant des fonctions de Green. Selon un mode de réalisation préférentiel de l'invention, le procédé comporte un calcul de perturbations dans lequel les variations spatiales du nombre d'onde effectif sont traitées comme une perturbation, les variations spatiales des propriétés optiques correspondant aux variations spatiales du nombre d'onde effectif. Description sommaire des dessins D'autres avantages et caractéristiques ressortiront plus clairement de la description qui va suivre de modes particuliers de réalisation de l'invention donnés à titre d'exemples non limitatifs et représentés aux dessins annexés, dans lesquels : 25 La figure 1 illustre la propagation de la lumière dans un milieu diffusant. La figure 2 représente l'organigramme d'un mode de réalisation particulier du procédé selon l'invention. 30 Description de modes particuliers de réalisation Sur la figure 1, une source S de lumière d'excitation ponctuelle est placée à un point rs et émet une lumière ayant une première longueur d'onde TSF et une amplitude Q. La lumière émise se propage dans un volume V, ici schématisé par un cylindre mais pouvant être de forme quelconque, d'un milieu diffusant dans lequel sont répartis des éléments fluorescents F selon une distribution que l'on cherche à déterminer. Sur la figure 1, une première onde diffusive LSF émise par la source S excite un élément fluorescent FI qui émet ensuite un rayonnement ayant une seconde longueur d'onde ?FD l'intensité de cette seconde onde diffusive LFD est mesurée par le détecteur D. La distribution des éléments fluorescents doit alors être déterminée à partir des mesures effectuées par le détecteur D. Ces mesures sont effectuées à partir d'un ensemble de positions de source. Dans le cas de plusieurs détecteurs, les mesures effectuées par chacun des détecteurs seront prises en compte. Sur la figure 1, deux autres éléments fluorescents F2 et F3 sont représentés. Les éléments fluorescents FI, F2 et F3 étant répartis dans le volume V, le détecteur reçoit une densité de photons mesurée, composée par toutes les secondes ondes LFD émises par l'ensemble des éléments fluorescents FI, F2 et F3. Pour des raisons de clarté, les ondes correspondant aux éléments fluorescents F2 et F3 ne sont pas représentées. On définit une première fonction de transfert d'énergie G(sF,is,FF) représentative du transfert de l'énergie entre la source S de lumière d'excitation ponctuelle et les éléments fluorescents F, et une seconde fonction de transfert d'énergie G(4D,rF,r'D) représentative du transfert de l'énergie entre les éléments fluorescents F et le détecteur D.30 On peut montrer que la densité de photons 0 mesurée par le détecteur D peut être exprimée, de façon approximative, par une intégrale comportant le produit des première et seconde fonctions de transfert G ~(rSr1)) x Q(rÇ). f G( $FrSrF)IJ(rF)G(ÂFD,rF',rD)d F (1), V où le paramètre N(rF) = i • bu(rF) (2) dépend du rendement: quantique ri de l'élément fluorescent et de l'absorption locale 6 t(?) due aux éléments fluorescents. La densité de photons 0 est, bien entendu, également proportionnelle à l'amplitude Q(rs) de la lumière émise par la source S. Les première et seconde fonctions de transfert G peuvent ainsi être traitées comme des fonctions de Green déterminées en résolvant des équations de diffusion pour le volume V fini et pour une longueur d'onde Â donnée : V2G(r,ï-o) + k2G(Â,,ro) = -4arb( -ro) (3), étant la fonction de Dirac et k~ = le nombre d'onde. Dans ces expressions, Â est le coefficient d'absorption du milieu, D, est le coefficient de diffusion, et ro sont les variables spatiales de la fonction de Green et V est l'opérateur de dérivée spatiale nabla. Les paramètres Â et D, sont évalués aux longueurs d'onde correspondantes XSF et 2 FD. De plus, les solutions des équations de diffusion doivent satisfaire à des conditions aux limites sur la surface délimitant le volume, par exemple les conditions de Dirichlet ou les conditions de Neumann. Les expressions ainsi trouvées pour les fonctions de Green sont ensuite introduites dans l'équation (1). Pour un calcul discret, l'espace du volume V est divisé en N volumes élémentaires selon un maillage approprié. L'intégrale de l'équation (1) est alors remplacée par une somme sur les N volumes élémentaires dv . constituant le volume V : N 0(rs,rD) G(À.sF,rs,rj.)%(1 j)G(ÀFD,ri,YD)dv. (4). Ainsi, on associe à chaque volume élémentaire dv1 une première fonction de transfert d'énergie G(ÀS.F,i-s,i-j ) et une seconde fonction de transfert d'énergie G(2FD,Yf,YD) . Lorsqu'on considère une pluralité de Ns sources S de lumière d'excitation ponctuelle et une pluralité de ND détecteurs D, une représentation matricielle peut être choisie : [O]Nä x\'s = [G], o xN [ 1] NxN [GSF ]NxNs [Q]N, (5). ~o Dans cette équation, chaque colonne de [0] représente la mesure sur NDxNs ND détecteurs pour une source S donnée. Pour chaque source S et pour chaque détecteur D, une matrice de conversion J peut être calculée, ce qui permet de relier de manière linéaire 15 la densité de photons 0 mesurée au paramètre (3 : N PSD = JSD,I; (6). J=l L'ensemble des combinaisons source-détecteur permet de construire l'équation matricielle qui est ensuite résolue dans un algorithme de reconstruction, soit en calculant l'erreur entre les mesures expérimentales et 20 cette équation théorique matricielle (par exemple par l'intermédiaire de la méthode de rétroprojection d'erreurs ou d'algorithmes de type ART (Algebraic Reconstruction Technique)), soit en inversant directement la matrice J (par exernple par l'intermédiaire d'algorithmes de type SVD (Singular Value Decornposition)). 25 Dans le cas du calcul d'erreur entre la mesure expérimentale et la théorie, on utilise une distribution théorique des éléments fluorescents qui est ajustée au fur et à mesure de manière à réduire cette erreur. Dans le cas de l'inversion de la matrice, on obtient, dans une première étape, le paramètre (3 qui dépend de la distribution des éléments fluorescents par l'intermédiaire de l'absorption locale 6 (rF) due aux éléments fluorescents. La connaissance du paramètre 13 permet ainsi de déterminer la distribution et la concentration locale des éléments fluorescents. Afin de tenir compte de l'inhomogénéité de l'objet, on pourrait théoriquement reconstruire d'une part, la distribution volumique du coefficient de diffusion D, et d'autre part, la distribution volumique du coefficient d'absorption du milieu , ou du nombre d'onde k,, et ensuite remonter à la répartition de la fluorescence. Cependant, la technique tomographique utilisant des signaux continus ne permet pas de discriminer le coefficient de diffusion et le nombre d'onde. Ainsi, il convient d'introduire dans l'équation de diffusion au moins un terme d'inhomogénéité tenant compte des variations spatiales des propriétés optiques et ce terme d'inhomogénéité est incorporé dans un nombre d'onde effectif k'. Le formalisme analytique à base de fonctions de transfert décrit ci-dessus peut ensuite être utilisé, même pour un milieu inhomogène et dans le cas de signaux continus, et les approximations induites par ce formalisme sont ensuite corrigées. En effet, il peut être montré qu'il n'est pas nécessaire de déterminer le coefficient d'absorption du milieu et le coefficient de diffusion, pour remonter à la répartition tridimensionnelle de fluorophores dans l'objet inhomogène. Ainsi, afin de déterminer la densité de photons dans le milieu diffusant, on se propose, par exemple, de résoudre l'équation v'O(/ YS,Y) i ka~,Y) + vD''(Y) 0 ~(YS,Y) - -Q(YS) ~( ) (7), D~(Y) DÂ(Y) dans laquelle le nombre d'onde initial k, (r') peut dépendre de la position. Le terme d'inhomogénéité VD' (r) • o correspond à la variation spatiale du coefficient de diffusion 1))ä On introduit alors un nombre d'onde effectif k'2 = kÂ(r) + ODE (r) • 0 dépendant DÂ(r) du nombre d'oncle initial k, (rr) et de la variation spatiale du coefficient de diffusion. On prend donc en compte ce terme en tant que composante additionnelle au nombre d'onde initial. Dans le cas d'une variation spatiale nulle du coefficient de diffusion, le nombre d'onde effectif k' serait égal au io nombre d'onde initial. Le nombre d'onde effectif k' est modifié au cours d'itérations de façon à tenir compte de toutes les inhomogénéités du milieu et de déterminer les fonctions de Green représentatives du milieu. Le terme d'inhomogénéité concerne le coefficient de diffusion, qui prend en 15 compte les propriétés de diffusion et d'absorption du milieu. On peut montrer que la densité de photons 0 peut être décomposée en ondes sphériques, notamment si l'on se place loin des interfaces à l'intérieur d'un milieu diffusant. Dans ces conditions, on démontre que D ù '(r) , VO(rs,r) = f(r) 0(rs,r) (8), ~!)a(r) où f (i) est une fonction dépendant des coordonnées spatiales. On peut donc effectuer un changement de variable conformément à l'expression kit = k,(r) + .f (@) (9). Dans tous les cas, on peut considérer que l'approximation faite lors de la 25 décomposition du signal en ondes sphériques peut être intégrée dans les valeurs du paramètre k' qui représente, de manière globale, l'atténuations D?(@) 20 L'équation 7 ci-dessus peut, par exemple, être résolue par des techniques numériques ou par des techniques perturbatives. Les techniques numériques, par exemple de type éléments finis ou de type différences finies, permettent notamment de traiter aussi bien les fortes variations de k' que les faibles variations. Une approche perturbative détaillée ci-dessous permet notamment d'utiliser le formalisme analytique introduit ci-dessus. En utilisant un calcul de perturbations, le signal d'excitation mesuré cl) se décompose en un signal de référence (po connu caractérisant un milieu io homogène dont on connaît la solution et un résidu perturbatif b~ correspondant à la variation spatiale 8k du nombre d'onde effectif k' décomposé conformément à l'expression : k'2 (F) = ko + bk2(i-) (10). Les variations spatiales du nombre d'onde effectif sont ainsi traitées comme 15 une perturbation. Les variations spatiales des propriétés optiques correspondent aux variations spatiales du nombre d'onde effectif k'. La densité de photons peut ainsi être écrite de façon approximative : 0(rs,i')_(l)o+bil) _ Q(~ s) Go(rs'r + if (rs,r)bk2Go(r,r')dF = 4z D), 4n ~o(rs,r + 1 Q(rs) /'G(rS,Y)bk2Go(Y,Y')dF 4z 4zDÂ 20 (11). Dans cette approche perturbative, la fonction de transfert G est déterminée par l'équation de Dyson : G(Fs,F') _ )+ J_ f G(-is,r)bk2Go(F,F')dF (12). Pour résoudre l'équation (12), on procède de manière itérative en remettant à 25 jour progressivement la fonction de transfert G. Ainsi, les fonctions de transfert pour un milieu inhomogène peuvent être déterminées par l'intermédiaire de l'algorithme représenté à la figure 2. Selon cet algorithme, le nombre d'onde effectif k' est d'abord approximé par un nombre d'onde non-perturbé ko et la fonction de Green G est d'abord approximée par une fonction de Green non-perturbée Go, définis lors d'une première étape F1. La densité de photons (ptf, théorique est déterminée (étape F2) conformément à l'équation 7, numériquement ou analytiquement. D'autre part, des mesures expérimentales (F3), notamment du signal d'excitation, et un traitement des données expérimentales (F4) permettent d'obtenir la densité de photons (pexp expérimentale. Un test de convergence (étape F5) permet de vérifier, si la densité de photons théorique et la densité de photons (éeXp expérimentale se rapprochent dans la mesure d'un paramètre d'écarts prédéterminé. Si le module de leur différence n'est pas inférieur au paramètre d'écart s (sortie Non de F5), la densité de photons (pth théorique est à nouveau déterminée avec un nombre d'onde effectif k' corrigé dans une étape F6. Si le test de convergence est positif (sortie Oui de F5), la fonction de Green G est corrigée lors d'une étape F7. Ensuite, la répartition des fluorophores est reconstruite (étape F8). Le signal de fluorescence mesuré dépend des propriétés optiques à la longueur d'onde d'excitation et à la longueur d'onde de fluorescence. Les longueurs d'onde d'excitation et de fluorescence étant généralement relativement proches, c'est-à-dire l'écart étant de quelques dizaines de nanomètres, on peut éventuellement négliger la différence de l'atténuation aux deux longueurs d'onde. Parfois, cependant, la variation des propriétés optiques en fonction de la longueur d'onde est connue pour les principaux constituants du milieu (eau, sang, particules diffusantes) et l'on peut extrapoler le résultat, par exemple, par une régression linéaire. Il est également possible d'effectuer les mesures aux deux longueurs d'onde, excitation et fluorescence, et d'effectuer la reconstruction de deux distributions volumiques d'atténuation.30 Revendications 1. Procédé de reconstruction de la distribution de propriétés optiques pour examiner un objet (1) diffusant inhomogène par éclairement de l'objet (1) avec une lumière continue, le procédé comportant la formulation d'une équation de diffusion comportant un nombre d'onde initial (kÀ(r)), procédé caractérisé en ce que, l'équation de diffusion comportant au moins un terme d'inhomogénéité tenant compte des variations spatiales des propriétés optiques, le terme d'inhomogénéité est incorporé dans un nombre d'onde 1 o effectif (k' ). 2. Procédé selon la revendication 1, caractérisé en ce que, la lumière continue étant une lumière d'excitation (k ex) de fluorophores (4) répartis dans l'objet (1), le procédé comporte la reconstruction d'une image de 15 tomographie de fluorescence. 3. Procédé selon l'une des revendications 1 et 2, caractérisé en ce que le carré du nombre d'onde effectif (k') est défini par la somme du terme d'inhomogénéité et du carré du nombre d'onde initial (k),(r)). 4. Procédé selon la revendication 3, caractérisé en ce que le terme d'inhomogénéité correspond à l'expression oD'(r) • V , DÂ étant le coefficient Dx(r) de diffusion dépendant de la coordonnée spatiale i. 25 5. Procédé selon l'une quelconque des revendications 1 à 4, caractérisé en ce que la répartition des fluorophores (4) est déterminée par l'intermédiaire d'un formalisme utilisant des fonctions de transfert. 6. Procédé selon l'une quelconque des revendications 1 à 5, caractérisé en 30 ce qu'il comporte un calcul de perturbations dans lequel les variations 11 20 spatiales du nombre d'onde effectif (k') sont traitées comme une perturbation, les variations spatiales des propriétés optiques correspondant aux variations spatiales du nombre d'onde effectif (k').5