La théorie des sondages - Images des mathématiques

ou juste avant la balise de fermeture -->

10% de réduction sur vos envois d'emailing --> CLIQUEZ ICI

ou juste avant la balise de fermeture -->

82 Introduction La théorie des sondages ne fait pas partie, en France, des connaissances usuelles des statisticiens même si la pratique des sondages, quant à elle, est très répandue. Elle s’est développée à partir des années 1930 dans le monde anglo-saxon ainsi qu’en Inde. Bien que reposant sur les mêmes principes que la statistique mathématique classique elle en diffère sensiblement dans son esprit en raison d’objectifs spécifiques. Cette théorie se consacre essentiellement au problème de la sélection de l’échantillon (voir ci-après la notion de plan de sondage qui fait le pendant des plans d’expériences en théorie classique) et à la recherche d’estimateurs. Du fait qu’elle porte sur des populations finies, d’existence bien concrète, elle ne peut ignorer les contraintes du monde réel, ce qui n’est peut-être pas sans lien avec le faible intérêt qu’elle suscite chez nous. Plan de sondage et probabilités d’inclusion On considère une population comprenant N individus parfaitement identifiés par un numéro d’ordre. Pour ce qui suit il nous suffira de ne retenir que ces numéros d’ordre et nous définissons ainsi la population U = {1,...,k,...,N}. Notons que les vocables de population et individus sont purement conventionnels. Les parties de cette population sont appelées échantillons. Dans cette présentation nous n’envisagerons que la situation où l’échantillon à sélectionner est de taille (cardinal) fixée, notée n, et désignerons simplement par S l’ensemble des échantillons de taille n. Commençons par donner quelques définitions. Définition 1. On appelle plan de sondage une loi de probabilité définie sur S. Concrètement le plan de sondage définit, pour chaque échantillon, la probabilité qu’il soit sélectionné via le mécanisme aléatoire utilisé. Définition 2. Soit un plan de sondage p et Sk l’ensemble des échantillons contenant l’individu k. On appelle probabilité d’inclusion de l’individu k ∈ U : πk = s∈Sk p(s) La théorie des sondages Michel LEJEUNE* On présente ici les principes et enjeux simples de la théorie des sondages pouvant susciter des désirs d'approfondissement. * Université Pierre Mendès France, LABSAD (BSHM), 38040 Grenoble Cedex 09. michel.lejeune@upmf-grenoble.fr Il s’agit donc de la probabilité que cet individu appartienne à l’échantillon sélectionné. Le fait d’avoir un échantillon de taille n se traduit par k∈U πk = n. En effet, soit Ik la variable indicatrice de la sélection de l’individu k, on a πk = E(Ik ) et : k∈U πk = k∈U E(Ik ) = E( k∈U Ik ) = n. On définit de même la probabilité qu’un couple d’individus {k,l} soit dans cet échantillon, en considérant l’ensemble Skl des échantillons contenant ce couple. Définition 3. On appelle probabilité d’inclusion d’ordre 2 des individus k et l : πkl = s∈Skl p(s). Pour illustrer ces notions considérons le plan simple sans remise (PSSR), correspondant au cas où p est uniforme sur S : ∀s ∈ S,p(s) = N n −1 . Comme il y a N − 1 n − 1 échantillons comprenant un individu donné, on a : ∀k ∈ U, πk = N − 1 n − 1 N n −1 = n N . Ce résultat, à savoir que la probabilité d’être sélectionné, pour un individu donné, est égale au taux de sondage n N , est tout à fait intuitif. On montre aisément que l’on a πkl = n(n − 1) N(N − 1) . Remarque 4. Notre définition d’échantillon exclut le cas du plan simple avec remise (PSAR), lequel correspond à la théorie classique de l’échantillonnage aléatoire. Les notions précédentes et la théorie générale peuvent aussi être développées dans cette situation. En pratique on n’effectue pas de plans avec remise car, on le sent bien intuitivement, le risque d’observer plusieurs fois le même individu constitue une perte d’information. La théorie montre, dans diverses situations, qu’un estimateur sans remise est meilleur qu’un estimateur avec remise et, même, qu’il est préférable, dans le cas d’un échantillon avec remise, de ne prendre en compte chaque individu qu’une seule fois. Les estimateurs On s’intéresse maintenant à une certaine « variable » réelle observable sur chaque individu. Notons yk la valeur prise par l’individu k. On souhaite estimer une caractéristique de la variable, le plus souvent son total (ou sa moyenne). Pour le total ty, Horvitz et Thompson ont proposé un estimateur pour un plan quelconque. Soit Y1,..., Yn les variables aléatoires correspondant à la sélection de n individus selon le plan de sondage, cet estimateur est : t H T y = n i=1 1 πi Yi . La théorie des sondages 83 84 Les 1 πi sont appelés poids de sondage. Cet estimateur est sans biais. En effet : t H T y = k∈U 1 πk Ik yk , d’où E(t H T y ) = k∈U 1 πk E(Ik )yk = k∈U yk = ty . Il est généralement retenu car on montre qu’il est le seul à avoir cette propriété parmi les estimateurs fonctions linéaires des Yi . De plus on sait exprimer sa variance et un estimateur sans biais de sa variance, via les probabilités d’inclusion d’ordres 1 et 2. On peut évidemment étudier cet estimateur par une voie directe. Pour les plans PSAR et PSSR il est clair que N n n i=1 Yi est sans biais pour ty (chaque Yi ayant une loi uniforme sur les yk , son espérance est my = 1 N ty , la moyenne de la population). Notons que l’estimateur de Horvitz-Thompson pour la moyenne de la population est la moyenne de l’échantillon Y = 1 n n i=1 Yi. Le calcul de la variance montre les difficultés inhérentes à la théorie des sondages. Pour le plan PSAR les Yi sont indépendantes et le calcul est immédiat. Pour le PSSR interviennent les covariances et le calcul demande une certaine adresse (le lecteur pourra tenter une démonstration en utilisant n i=1 Yi = k∈U Ik yk , Var(Ik ) = πk (1 − πk ) et Cov(Ik ,Il) = πkl − πkπl). On trouve : Var(Y) = N − n N − 1 v2 n où v2 est la variance de la population : v2 = 1 N N k=1(yk − my )2. Un apport essentiel de la théorie des sondages est d’établir la précision d’une estimation (vulgairement : la fourchette). Celle-ci est définie par la demi-longueur d’un intervalle de confiance à 95 %. En vertu du théorème central limite qui donne une bonne approximation pour les tailles usuelles d’échantillons, cet intervalle repose sur la loi de Gauss (moyennant une deuxième approximation due à l’estimation de la variance v2). Dans le cas PASR on adopte ainsi pour précision sur l’estimation de la moyenne : 1, 96 1 − n N s2 n où s2 est la variance de l’échantillon qui estime sans biais N N − 1 v2. Le principe de stratification Dans cette brève présentation il n’est pas possible de présenter les principaux plans de sondage. Nous nous concentrons sur le plan stratifié en raison de ses liens avec les pratiques très répandues d’utilisation de quotas et de redressement d’échantillons. Ces procédures reposent sur la même idée : toute partition de la population en strates fortement homogénes (variances internes aux strates faibles vis-à-vis de la variance globale) doit permettre d’accroître la précision. Dans le plan stratifié (PSTRAT) on effectue des sondages, par exemple de type PSSR, indépendamment dans chaque strate et, de facto, se pose le problème du choix des tailles des « sous-échantillons ». L’estimateur naturel de ty qui est aussi celui d’Horvitz-Thompson est obtenu, dans le cas de sous-échantillons de type PSSR, en recomposant les estimateurs des moyennes par strate, soit H h=1 NhY h pour les H strates. Par l’indépendance mutuelle sa variance est immédiate. Un problème intéressant, mais purement théorique, est celui dit de « l’allocation optimale », à savoir de trouver les tailles nh minimisant cette variance sous la contrainte H h=1 nh = n (pour le lecteur intéressé : les taux de sondage nh/Nh sont alors proportionnels aux écarts-type des strates). Généralement on stratifie à taux de sondage constants dans les strates ce qui garantit un gain de précision par rapport à un sondage PSSR. En effet la « fourchette » est multipliée par un facteur (1 − η2)1/2 où η2 est un coefficient bien connu des statisticiens, compris entre 0 et 1, qui mesure en quelque sorte le lien existant entre le critère de stratification et la variable d’intérêt (rapport de la variance intra-strates à la variance totale pour cette variable). Cette stratification, dite aussi « à la proportionnelle », est qualifiée par les praticiens comme produisant un échantillon représentatif vis-à-vis du critère de stratification. Notons que dans ce cas particulier les probabilités d’inclusion sont identiques pour tous les individus. Ceci nous amène à la « méthode des quotas » dont la pratique est systématique dans les instituts. Elle consiste à sélectionner un échantillon qui soit un modèle réduit de la population sur certains critères de partitionnement dont on peut penser qu’ils ont un lien avec la thématique de l’enquête. Comme les effectifs des strates dans la population doivent être connus pour chaque critère d’une part et qu’il y a de fortes contraintes de mise en œuvre d’autre part, on se limite généralement à un critère géographique (région), à la classe d’âge, au sexe et, parfois, à la catégorie socio-professionnelle. La différence avec un plan PSTRAT tient au fait que, pour des raisons de faisabilité, la proportionnalité des effectifs de l’échantillon à ceux de la population ne concerne que les effectifs des critères « à la marge » et non des critères croisés entre eux. Il n’est pas inutile, ici, pour les citoyens que nous sommes tous, d’ouvrir une parenthèse pratique. Dans la presse on déclare communément que les résultats d’un sondage proviennent d’un échantillon obtenu par la méthode des quotas et/ou d’un échantillon représentatif de la population selon la région, l’âge, etc. Ainsi il y a, chez les sondeurs, un véritable culte des quotas, la plupart d’entre eux étant convaincus qu’il suffit d’effectuer ce modèle réduit pour garantir de bonnes estimations. Or l’essentiel n’est pas là et, en fait, la mention consacrée dans la presse n’est d’aucun intérêt quant à juger de la qualité réelle de l’échantillon. Ce qui nous importe avant tout est de savoir dans quelle mesure un plan aléatoire, incluant ou non des contraintes de quotas, a-t-il été respecté. L’approche théorique permet de déterminer le gain apporté par l’utilisation de quotas par rapport à un PSSR. Il est de même nature que celui présenté en stratification à ceci près que le coefficient η2 repose ici sur la variance expliquée par un modèle à effets additifs des critères de quotas (pour les initiés : modèle d’analyse de variance sans interactions). Le constat empirique est bien décevant pour les critères usuels, le facteur de réduction de la fourchette (1 − η2)1/2 ne passant que très rarement sous la valeur 0,90. Ce constat indique au passage que ces critères de quotas sont des déterminants extrêmement ténus du comportement ou des opinions des individus (il n’en reste pas moins que les quotas sont utiles, voire même nécessaires, pour limiter les biais importants découlant des obstacles de terrain). La pratique des redressements relève des mêmes principes que les quotas mais elle intervient en aval du recueil de l’échantillon. Elle consiste, dans le but d’améliorer les estimations, à « caler » l’échantillon observé sur la population pour divers critères disponibles. Comme pour les quotas le calage se fait à la marge, critère par critère. Il s’effectue en attribuant des poids aux individus en perturbant le moins possible leurs poids de sondage. Mathématiquement il s’agit de déterminer les poids wi tels que les contraintes de marges soient respectées pour les différents critères retenus et tels que : n i=1 d wi, 1 πi soit minimal (où d est la distance choisie). La résolution d’un tel problème nécessite une procédure itérative. En fait l’algorithme utilisé en pratique a été introduit tout à fait empiriquement : la recherche de poids correcteurs étant immédiate pour un seul critère, on calera successivement chacun des critères de façon cyclique jusqu’à quasi convergence vers les bonnes marges. Ce n’est que postérieurement que les travaux initiés par J.-C. Deville à l’INSEE ont donné un cadre théorique aux procédures de calage, notamment en introduisant la notion de distance mentionnée ci-dessus (l’algorithme usuel relève d’une distance implicite peu commune mais proche de la distance quadratique). Les estimateurs par redressement n’ayant pas d’expression explicite leurs propriétés sont difficiles à établir. Ainsi, pour déterminer biais et variance, on doit se contenter d’approximations pour n grand. Brièvement disons que le La théorie des sondages 85 86 biais, pour des distances raisonnables, reste négligeable et que leur précision est très proche de celle découlant des estimateurs avec quotas (rapport égal à 1 + O(1/n1/2)). Ainsi, comme pour ces derniers, le gain par rapport à un PSSR, pour une variable d’intérêt donnée, sera d’autant plus élevé que les critères retenus auront un lien fort avec elle. Ceci permet d’orienter le choix de ces critères sur la base des liens observés dans l’échantillon. La modélisation des non-réponses Il existe diverses sources d’erreur dans les sondages : erreur de mesure, biais de sélection, biais de couverture de la population et, souvent la plus importante, la non-réponse, à savoir qu’un individu dûment sélectionné par le plan de sondage n’a pu être observé pour une raison ou une autre. Les non-réponses, systématiquement présentes et même dans de fortes proportions, sont une cause potentielle de biais. Pour y remédier différents schémas du mécanisme de non-réponse ont été proposés pour lesquels s’appliquent des modèles appropriés. Par exemple, si le fait de répondre ou non à une enquête est indépendant de la variable d’intérêt conditionnellement à un ensemble de variables auxiliaires (par exemple âge, sexe, profession, statut matrimonial, etc.) le recours à un redressement sur ces variables est efficace. Il existe aussi des non-réponses partielles (seules certaines variables n’ont pu être observées ici ou là) débouchant sur des modèles d’imputation. Le traitement des non-réponses totales ou partielles donne lieu actuellement à de nombreux développements théoriques. Pour en savoir plus ARDILLY (P.), Les techniques de sondage, Technip, (1994). DEVILLE (J.-C.) et SARNDAL (C.-E.), Calibration estimators in survey sampling, Journal of the American Statistical Association, vol. 87, p. 376-382, (1992). TILLÉ (Y.), Théorie des sondages, Dunod, (2001). 64 comme marginal par rapport aux « grandes mathématiques ». Il suivait alors le séminaire d’algèbres d’opérateurs de Jacques Dixmier à l’Ecole Normale Supérieure, où l’on parlait notamment d’algèbres de von Neumann. Ces algèbres sont des géné- ralisations non commutatives, ou si l’on préfère, quantiques, de la théorie de la mesure. Elles avaient été introduites par von Neumann dès les années 20 ou 30 pour donner un fondement mathématique à la mécanique quantique qui venait d’être découverte. Etant donné un espace de Hilbert H on considère l’algèbre L(H) des opérateurs bornés sur H. On considère des algèbres d’opérateurs, c’est-à-dire des sous-algèbres A de L(H) telles que si un opérateur T appartient à A, il en est de même de son adjoint T ∗. On impose de plus à la sous-algèbre A une condition topologique : on dit que A est une algèbre de von Neumann si elle est stable par convergence faible ou forte (c’est la convergence simple sur tout vecteur de H) ; on dit que A est une C∗ -algèbre si elle est stable pour la convergence normique (c’est la convergence uniforme sur la boule unité de H). Parmi les algèbres de Von Neumann, un rôle fondamental est joué par les facteurs, c’est-à-dire les algèbres de von Neumann dont le centre est réduit aux scalaires. Murray et von Neumann, dans les années 1930, avaient tenté une classification des facteurs, qu’ils répartirent en trois classes. Les facteurs de classe I et II sont les plus proches du cas commutatif, avec une notion de trace et de dimension : une trace est une forme linéaire positive τ qui vérifie la propriété τ(xy) = τ(yx). Le type I est celui des algèbres de matrices Mn(C), ou en dimension infinie, de l’algèbre L(H) des opérateurs bornés sur un espace de Hilbert. On a alors la trace usuelle des Alain Connes : une autre vision de l’espace Il est difficile de présenter Alain Connes sans commencer par évoquer brièvement son impressionnant curriculum vitae. Né en 1947, élève de l’Ecole Normale Supé- rieure (1966-1970), après un passage au CNRS, il est successivement professeur à Paris 6, directeur de recherche au CNRS, et depuis 1984 professeur au Collège de France, où il occupe la chaire d’Analyse et Géométrie. Parallèlement, il est depuis 1979 professeur à l’Institut des Hautes Etudes Scientifiques à Bures-surYvette. Il partage sa vie mathématique entre ces deux lieux : l’IHES lui offre le calme pour s’adonner à ses recherches et la possibilité de rencontrer des mathématiciens et des physiciens théoriciens de tous pays et de toutes spécialités ; le Collège de France lui donne l’occasion, chaque année, de présenter dans un cours ses résultats les plus récents. Depuis 2003 il est aussi professeur à l’Université de Vanderbilt aux Etats-Unis. Il a reçu les plus hautes distinctions internationales : la médaille Fields (1982) et le prix Crafoord (2001), deux prix dont le prestige est comparable au prix Nobel, qui, on le sait, n’existe pas pour les mathé- matiques. En France, où il est membre de l’Académie des Sciences depuis 1983, il vient de recevoir, en décembre 2004, la médaille d’or du CNRS. La thèse : classification des facteurs de type III Pour comprendre qui est Alain Connes, il faut discerner quelles furent les différentes étapes de sa vie mathématique. Dans un premier temps, il s’est imposé comme un jeune mathématicien au talent exceptionnel, en résolvant un problème reconnu difficile par les spécialistes, mais considéré par beaucoup à l’époque opérateurs, qui a la propriété d’intégralité : elle prend des valeurs entières sur les projecteurs (éléments p tels que p2 = p = p∗). Mais il y a aussi des facteurs admettant une trace qui prend des valeurs réelles quelconques sur les projecteurs : c’est le type II. Murray et von Neumann découvrirent aussi l’existence d’objets hautement non commutatifs, dit facteurs de type III, qui n’admettent aucune trace. Jusqu’aux années 1970, ces facteurs étaient restés mystérieux et résistaient à toute tentative de classification. Lorsque Alain Connes arrive au séminaire Dixmier, on en est encore là, mais les travaux de Powers, Araki et Woods ont produit de nouveaux exemples de facteurs de type III, et même une infinité de tels facteurs deux à deux non isomorphes. Le génie d’Alain Connes a été d’appliquer à ces objets une théorie encore nouvelle et peu exploitée, due au mathématicien japonais Minoru Tomita. Un facteur de type III n’ayant pas de trace, on remplace la notion de trace par celle de poids. Un poids est une forme linéaire positive ϕ, mais a priori ϕ(yx) =/ ϕ(xy). La surprise, c’est que la non commutativité engendre une dynamique, une évolution au cours du temps donnée par un groupe d’automorphismes à un paramètre du facteur. Plus précisément, on a pour tous les t ∈ R des automorphismes σt du facteur M, avec la loi de groupe σt ◦ σs = σt+s, qui permettent de corriger la non commutativité via la formule dite KMSβ : ϕ(yx) = ϕ(xσiβ(y)) où σiβ s’obtient par prolongement analytique, pour un certain β > 0. Connes a montré que ce groupe à un paramètre est en fait indépendant du poids, modulo les automorphismes inté- rieurs, et donne lieu à des invariants spectraux qui permettent de classifier les facteurs de type III. Connes a ainsi introduit dans sa thèse (1973) les facteurs dits IIIλ où λ est un nombre réel entre 0 et 1. Cette solution, par un jeune thésard encore inconnu, d’un problème ouvert depuis des décennies a profondément impressionné Jacques Dixmier et l’ensemble des spécialistes des algèbres d’opérateurs. Cependant, la majorité des mathématiciens français ignoraient alors la théorie de von Neumann, et le problème de la classification des facteurs n’était pas considéré comme un des principaux défis des mathématiques. Ce sont surtout les physiciens théoriciens qui ont reconnu le génie d’Alain Connes. La mécanique statistique quantique, ainsi que certains modèles de théorie des champs, utilisent en effet les algèbres de von Neumann de façon essentielle. La statistique de Boltzmann associée à un hamiltonien H est donné (sur une observable A) par le poids ϕ(A) = Trace(Ae−βH )/Trace(e−βH ) et l’évolution dans le temps par les automorphismes σt(A) = eitH Ae−itH . Le lecteur vérifiera sans peine la formule KMSβ ci-dessus. Le groupe à un paramètre associé à un poids modélise l’évolution dans le temps d’un système statistique quantique associé à un état de température donné (comme d’habitude, β = 1/kT où k est la constante de Boltzmann et T la température absolue). Aussi, lorsque Connes devient visiteur à l’IHES, c’est en tant que physicien théoricien. Des facteurs aux feuilletages : vers la géométrie non commutative Cette arrivée de Connes à l’IHES est un tournant de sa carrière scientifique. En effet, il entre en contact avec des mathématiciens qui, s’ils ignorent tout de la théorie des facteurs de type III, jonglent quotidiennement avec des faisceaux et des feuilletages, avec des groupes d’homologie et des tenseurs de courbure. La géométrie différentielle et la géométrie algébriques sont considérées par la communauté mathématique comme plus centrales, pour ne pas dire plus nobles, que la théorie des algèbres d’opérateurs. Mais Connes ne se laisse pas impressionner. Il raconte, avec sa modestie habituelle, et quelqu’exagération, qu’il ne comprenait rien aux conversations de ses collègues avec qui il déjeûnait à la cafétéria de l’Institut. Mais qu’on ne s’y trompe pas, Alain Connes a une faculté étonnante d’assi- 65 miler de nouvelles notions. Et pour apprendre une théorie mathématique, au lieu de se plonger dans les livres, il pré- fère discuter avec d’autres mathématiciens, se faire expliquer puis retrouver par lui-même toute la théorie. Et il fait cela très vite. Aussi a-t-il très rapidement compris la théorie des feuilletages, et vu le lien avec la théorie de von Neumann. Un feuilletage est un objet géométrique qui est localement trivial comme un mille-feuille, mais qui globalement a une structure dynamique non triviale. Connes définit alors l’algèbre de von Neumann associée à un feuilletage, qui est en général un facteur, souvent de type III, le groupe à un paramètre de la théorie de Tomita ayant une interprétation géométrique très simple. Il définit aussi la C∗ -algèbre associée à un feuilletage, et dans l’esprit d’Atiyah et Singer montre un théorème d’indice pour les opérateurs elliptiques le long des feuilles d’un feuilletage transversalement mesuré. D’où l’idée de géométrie non commutative. L’espace des feuilles d’un feuilletage est un ensemble non dénombrable sur lequel on serait bien en peine de définir, au sens classique, une théorie de la mesure, ou une topologie, voire une structure différentiable. Expliquons cela par un exemple, celui du feuilletage du tore bidimensionel par une droite de pente irrationnelle. On se donne un carré dans le plan ; après recollage des côtés opposés, on obtient un tore, qui est une variété compacte de dimension 2. On se donne une direction D dans le plan. Partant d’un point x1 du bord inférieur I, on construit une suite de segments de droites dans le carré de la façon suivante : le premier segment est issu de x1 ; chaque segment est parallèle à D et a son origine et son extrêmité sur le bord du carré ; chaque fois qu’un point du bord est extrêmité d’un segment, l’origine du segment suivant est le point opposé du bord. On obtient ainsi une famille de segments tracés à l’intérieur du carré, et après recollage des bords, on obtient une trajectoire sur le tore, appellé feuille (du point de vue de la géométrie riemannienne, c’est une géodésique de la métrique plate sur le tore). Le feuilletage est donné par la partition ainsi obtenue de la variété (ici le tore) en feuilles. Voir les figures ci-contre où l’on a noté x2,x3,··· les points obtenus sur le segment I. Il faut distinguer deux cas. Supposons que la pente de D soit rationnelle. Le processus est périodique (sur la figure, la pente est égale à 2). Il n’y a qu’un nombre fini de segments, la trajectoire sur le tore est périodique. L’espace des trajectoires est bien décrit par l’algèbre commutative des fonctions continues sur un intervalle (ici la moitié de l’intervalle I). C’est un espace usuel, dit commutatif. Prenons au contraire pour la pente de D un nombre θ irrationnel. On a tracé les quatre premiers points situés sur I de la trajectoire de x1 et les deux premiers points de celle de y1. Chaque trajectoire est infinie (on ne revient jamais au point de départ) et dense. Si on veut décrire l’espace des trajectoires par des fonctions continues sur I, celles-ci doivent avoir la même valeur sur les points d’une même trajectoire, donc être constantes. L’algèbre de ces fonctions est l’algèbre des nombres complexes. Elle ne donne aucun renseignement sur l’espace Tθ des feuilles. Alain Connes propose de décrire Tθ par un algèbre non commutative. En ce sens, l’espace des feuilles est un espace non commutatif. Un autre exemple d’espace non commutatif est l’espace des orbites de l’action d’un groupe discret sur une variété compacte. Le cas le plus simple, intimement lié au feuilletage ci-dessus, est celui où la variété est le cercle muni de l’action du groupe Z engendrée par une rotation d’angle 2πθ. Là encore les orbites sont denses si θ est irrationnel et l’espace des orbites doit être décrit par une algèbre non commutative. Pour définir l’algèbre du feuilletage du tore ci-dessus, on considère des noyaux k(x,y) définis sur les couples de points x et y situés sur une même feuille, continus et à support compact sur chaque feuille. On multiplie ces noyaux par le produit usuel de convolution des noyaux : k1 ∗ k2(x,y) = k1(x,z)k2(z,y)dz où l’intégrale est prise sur la feuille, pour la mesure de Lebesgue. De tels noyaux peuvent être interprêtés comme des familles d’opérateurs indexés par les feuilles. En complétant l’algèbre ainsi obtenue, on définit alors d’une part l’algèbre de von Neumann du feuilletage, et d’autre part la C∗-algèbre du même feuilletage. De même, dans le cas de l’action de Z on considère des matrices à support fini ai,j indexées par des couples (i,j) d’éléments d’une même orbite, et qu’on multiplie par le produit usuel des matrices, et on définit ainsi l’algèbre de von Neumann et la C∗- algèbres associées. 66 67 Notons que dans les deux exemples ci-dessus les facteurs obtenus sont de type II. Cela est dû au fait que l’action de la rotation sur le cercle préserve la mesure de Lebesgue (ou pour le feuilletage, le flot des trajectoires préserve une mesure transverse), d’où une trace sur l’algèbre. Mais on construit facilement des exemples sans mesure invariante, et on obtient alors le type III. Figure 1 – Pente irrationnelle. Figure 2 – Pente rationnelle. 68 Ces géométries « non commutatives » sont donc décrites au moyen d’algèbres non commutatives qui jouent le rôle d’espaces de fonctions. Ainsi un espace mesuré X est décrit par l’algèbre L∞(X) des fonctions mesurables bornées, qui est une algèbre de von Neumann commutative. Pour l’espace des feuilles d’un feuilletage, c’est l’algèbre de von Neumann (non commutative) associée qui joue le rôle d’algèbre de fonctions mesurables bornées sur l’espace des feuilles. De même pour la topologie : un espace topologique compact X est décrit par l’algèbre des fonctions continues sur X. Pour un feuilletage, la C∗ -algèbre (non commutative) associée est considérée comme l’algèbre des fonctions continues sur l’espace « non commutatif » des feuilles. Ce qui est fait pour un feuilletage peut aussi se faire pour l’action d’un groupe discret sur un espace compact. L’espace des orbites est lui aussi un espace « non commutatif » décrit par une C∗ -algèbre (du point de vue topologique) ou une algèbre de von Neumann. D’autres espaces non commutatifs sont des espaces d’orbites de relations d’équivalences, ou bien sont définis par des groupoïdes. Un autre cas intéressant est celui des groupes : le dual d’un groupe, c’est-à-dire l’espace des classes d’équivalences de représentations unitaires irréductibles, est aussi un espace non commutatif : dans le cas des groupes de Lie ou des groupes p-adiques, cet espace désingularise le dual décrit par la théorie des représentations, mais dans le cas des groupes discrets, on obtient un espace non commutatif hautement non trivial et encore mal connu. K-théorie et homologie cyclique Un invariant de topologie algébrique qui se généralise sans difficulté au cas non commutatif est la K-théorie. C’est un groupe abélien K(X) attaché à un espace topologique X (disons, compact) et qui classifie (à isomorphisme stable près) les fibrés vectoriels sur X. Deux fibrés E1 et E2 sur X sont stablement isomorphes s’il existe un fibré F tel que les fibrés sommes directe E1 ⊕ F et E2 ⊕ F soient isomorphes. Sur l’ensemble des classes d’isomorphisme stable de fibrés, on définit l’addition par la somme directe des fibrés, est K(X) n’est autre que le groupe (dit de Grothendieck) des différences formelles de classes de fibrés, exactement comme dans la construction de l’anneau Z des entiers relatifs à partir des entiers naturels. On sait depuis Atiyah et Singer que cet invariant joue un rôle crucial dans la théorie de l’indice. Par un théorème dû à Serre, les fibrés sur X correspondent à certains modules sur l’algèbre C(X). On peut alors définir la K-théorie d’une C∗ -algèbre et donc d’un espace topologique non commutatif. Dans le début des années 1980, Connes s’est intéressé à cet invariant, et a compris, avec Georges Skandalis, l’intérêt de la K-théorie bivariante développée à ce Figure 3 – Feuilletage du tore. moment là par le mathématicien russe Gennadi Kasparov. Ils ont montré, dans l’esprit de Grothendieck et en réinterprétant une des preuves d’Atiyah et Singer, comment le théorème de l’indice pour les opérateurs elliptiques sur des variétés compactes se ramenait à une propriété de fonctorialité en théorie de Kasparov. De là la généralisation aux feuilletages était naturelle. Il devenait donc naturel de calculer la K-théorie d’espaces non commutatifs comme les feuilletages, les groupes ou les actions de groupes. Dès 1980 lors d’un congrès à Kingston (Ontario), Alain Connes a eu l’intuition d’une interprétation géométrique de la K-théorie de la C∗ -algèbre d’un feuilletage ou d’un groupe. De ses discussions avec le topologue Paul Baum, qu’il rencontre alors, jaillira l’idée d’une conjecture, désormais célèbre sous le nom de conjecture de Baum-Connes : le groupe de K-théorie analytique (c’est-à-dire la K-théorie de la C∗ -algèbre du groupe ou du feuilletage) est isomorphe, via une flèche d’indice, à un groupe de Kthéorie dit géométrique, construit à partir du classifiant des actions propres du groupe (ou du groupoïde d’holonomie du feuilletage). Cette conjecture est liée à la théorie des représentations des groupes de Lie (séries discrètes), à la topologie (conjecture de Novikov sur l’invariance par homotopie des hautes signatures), à la géométrie riemannienne (conjecture de Gromov-Lawson sur les obstructions à l’existence de métriques riemanniennes à courbure scalaire positive), à l’algèbre (conjecture des idempotents). Précisée par la suite avec N. Higson, la conjecture de Baum-Connes sera le point de départ des travaux de nombreux mathématiciens pendant au moins 20 ans. Parmi les résultats les plus spectaculaires, citons : – le travail de N. Higson et G. Kasparov en 1996 qui établit la conjecture pour les groupes ayant la propriété dite de Haagerup ou a-T-menabilité (Gromov) ; c’est une classe de groupes relativement vaste, elle contient les groupes moyennables, mais exclut par exemple les groupes discrets possédant une propriété de rigidité dite propriété T de Kazhdan. – la thèse de Vincent Lafforgue en 1998, qui pour la première fois prouve la conjecture de Baum-Connes pour certains groupes discrets ayant la propriété T. Elle a également permis de démontrer la conjecture pour tous les groupes localement compacts connexes. Malgré le nombre de résultats établissant la conjecture pour des groupes particuliers, un cas aussi simple que celui de SL(3,Z), le groupe des matrices 3 sur 3 à coefficients entiers et de déterminant égal à 1 (ce groupe a la propriété T), est ouvert, sans qu’on ait aucune idée d’approche. Inversement, on peut tenter de construire des groupes suffisemment sauvages pour être des contrexemples à la conjecture. L’idée puissante de M. Gromov de construire des groupes aléatoires en sorte qu’ils aient des propriétés très éloignées des groupes pour lesquels la conjecture est actuellement prouvée, est prometteuse. Mais pour l’instant , elle n’a réussi à fournir que des contrexemples à une conjecture, dite conjecture de Baum-Connes à coefficients, qui est plus forte que la conjecture de Baum-Connes classique. Mais la K-théorie n’est pas le seul invariant intéressant en topologie algébrique. Une question naturelle est de définir l’homologie, ou la cohomologie d’un espace topologique non commutatif. Pour cela, il faut faire un détour par la géométrie différentielle. On sait que dans le cas classique, l’homologie (ou la cohomologie) peut être définie comme (co)homologie singulière (au moyen de triangulations par des simplexes ou des cycles singuliers) ou bien de façon équivalente comme cohomologie de C˘ ech, décrite au moyen de cocycles associés à un recouvrement par une famille d’ouverts. Il s’agit là de définitions topologiques de la (co)homologie à coefficients entiers, naturellement invariantes par homéomorphisme, mais difficilement généralisables dans le cas non commutatif. Si l’espace topologique est de plus muni d’une structure de variété différentiable, alors il y a une autre définition de la cohomologie (à coefficients réels ou complexes), celle de de Rham, obtenue à partir du complexe des formes différentielles. C’est cette dernière définition qui est retenue par Connes pour le cas non commutatif. Mais le prix à payer, c’est qu’il faut choisir une certaine sous-algèbre dense de la C∗ -algèbre, jouant le rôle de l’algèbre des fonctions lisses (de classe C∞, par exemple). Le point de départ est dans des calculs de géométrie non commutative : le calcul de caractères de Chern de modules de Fredholm fait apparaître une généralisation de la notion de trace. Ainsi un 2-cocycle cyclique est une forme trilinéaire sur une algèbre vérifiant les formules : τ(a0, a1, a2) = τ(a1, a2, a0) τ(a0a1,a2,a3) − τ(a0,a1a2,a3) + τ(a0,a1,a2a3) − τ(a3a0,a1,a2) = 0 A partir de la notion de n-cocycles, Alain Connes définit en 1981 la cohomologie cyclique d’une algèbre. Il développe cette théorie purement algébrique, découvre la longue suite exacte qui permet de la calculer. L’homologie 69 cyclique a depuis été abondemment utilisée par les algébristes, indépendemment des motivations d’Alain Connes qui, lui, revient toujours à son idée : la géométrie non commutative. L’un des problèmes les plus difficiles est de bien choisir la sous-algèbre dense de la C∗ -algèbre : il faut qu’elle soit suffisamment petite pour qu’on puisse calculer sa cohomologie cyclique, mais suffisamment grosse pour avoir la même K-théorie que la C∗ -algèbre. Une bonne partie des articles d’Alain Connes dans les années 1980 tournent autour de cette problématique, appliquée au cas des feuilletages ou au cas essentiellement équivalent des actions de groupes discrets sur des variétés. Il définit dans ce cadre la classe fondamentale transverse, analogue de la classe fondamentale d’une variété (qui appartient à l’homologie de la variété). En particulier, il donne une interprétation de la classe de Godbillon-Vey d’un feuilletage de codimension 1 en terme de cohomologie cyclique et donc d’accouplement avec la K-théorie. Un feuilletage de codimension 1 est donné par une équation ϑ = 0, où ϑ est une 1-forme satisfaisant la condition d’intégrabilité dϑ ∧ ϑ = 0. On a donc dϑ = α ∧ ϑ où α est une 1-forme. On considère alors la 3-forme α ∧ dα, et un calcul simple montre que sa classe de cohomologie ne dépend que du feuilletage. C’est l’invariant de Godbillon-Vey du feuilletage. Alain Connes tire de son interprétation de cette classe un corollaire frappant : si la classe de Godbillon-Vey est non nulle, alors le flot des poids de l’algèbre de von Neumann préserve une mesure de masse finie (et elle est de type III). Ce superbe théorème de Connes peut aussi se démontrer de manière élémentaire sans aucune homologie cyclique, mais il montre toute la force et la beauté du point de vue géométrie non commutative. On voit là une des caracté- ristiques de la pensée d’Alain Connes qui est sa profonde unité. En développant la K-théorie et l’homologie cyclique des feuilletages, il n’oublie pas son point de départ, la classification des facteurs de type III. Un autre succès important de cette méthode de géométrie différentielle non commutative est d’avoir donné la première preuve de la conjecture de Novikov pour les groupes hyperboliques au sens de Gromov. Dans ce cas l’algèbre de « fonctions lisses » sur le dual du groupe (vu comme espace non commutatif) est l’algèbre de fonctions à décroissance rapide définie par Paul Jolissaint. Triplets spectraux et indice transversal Jusque là, les théorèmes d’indices rencontrés dans le cadre des feuilletages sont des théorèmes d’indice longitudinaux : on considère des opérateurs elliptiques le long des feuilles du feuilletage, et l’indice d’un tel opérateur est un élément de la K-théorie de la C∗ -algèbre du feuilletage, c’est-à-dire la K-théorie de l’espace non commutatif des feuilles. Pour obtenir un indice qui soit un nombre, on l’évalue sur une classe d’homologie de cet espace, c’est-à- dire de la cohomologie cyclique de l’algèbre du feuilletage. Mais il y a un problème plus difficile auquel Alain Connes va s’attaquer dans les années 1990, en collaboration avec Henri Moscovici. C’est la question du théorème d’indice transverse : cette fois on veut vraiment considérer un opérateur elliptique sur l’espace non commutatif des feuilles, c’est donc un opérateur transversalement elliptique, dont Connes et Moscovici montrent qu’il définit un élé- ment de la K-homologie du feuilletage (la K-homologie est la théorie duale de la K-théorie) , donc une application de la K-théorie vers Z, l’anneau des entiers. Le but du théorème d’indice transversal est de donner cette application de la K-théorie vers Z de façon concrète, c’est-à-dire au moyen d’un cocycle cyclique, pour lequel une formule explicite est donnée. Il y a d’abord une première difficulté : fabriquer de tels opérateurs transversalement elliptiques, sans aucune hypothèse sur le feuilletage : ainsi on ne veut pas se restreindre au cas où le feuilletage aurait une métrique riemannienne transversale invariante par holonomie. Pour cela, Connes et Moscovici procèdent en deux étapes ; d’abord on grossit l’espace pour prendre celui de toutes les métriques transversales, ensuite on admet, au lieu d’opérateurs elliptiques, des opérateurs hypoelliptiques. Moyennant quoi on obtient un triplet spectral sur l’espace non-commutatif des feuilles. Etant donné un espace non commutatif dont la topologie est décrite par une C∗ -algèbre A, on se donne un espace de Hilbert H dans lequel est représentée A, et un opérateur (non borné) autoadjoint D à résolvante compacte (l’opérateur 1 + D2 est d’inverse compact) et tel que les commutateurs [D, a] soient bornés pour a dans une certaine sous-algèbre dense de A. Dans le cas de l’opérateur de Dirac sur une variété riemannienne, ce commutateur redonne la métrique, c’est-à-dire l’élément de longueur infinitésimal ds2. La notion de triplet spectral permet de définir l’analogue non commutatif de la notion de variété riemannienne. 70 71 Il reste alors à calculer le caractère de Chern de ce triplet spectral. Ce qu’ils font en utilisant la notion de trace de Dixmier et de résidu de Wodzicki. La formule obtenue est à la fois très simple par son élégance, et très compliquée par le nombre de termes impliqués dès qu’on veut l’expliciter. Ainsi, même dans le cas d’un feuilletage de codimension 1, il faut une bonne centaine de pages pour mener le calcul ... Dans le cas général, cela est quasiment impossible, sauf si l’on peut comprendre un principe permettant d’organiser et de simplifier ces calculs. Comme c’est en général le cas en mathématiques, un tel principe simplificateur est fourni par la notion de symétrie. Classiquement, la symétrie est décrite par un groupe ; ici, dans une situation non commutative, c’est une algèbre de Hopf ou groupe quantique : il s’agit d’un objet qui se comporte comme un groupe, sauf que l’ensemble des éléments du groupe n’est pas vraiment un ensemble ou un espace, mais un espace non-commutatif. L’idée est que ce groupe quantique agit sur l’espace non commutatif, d’où l’on déduit une application caractéristique qui permet de pousser la cohomologie cyclique du groupe quantique (qui se calcule, c’est la cohomologie de Gelfand-Fuchs) dans la cohomologie cyclique de notre espace non commutatif. Il se fait alors que le caractère du triplet spectral est dans l’image de cette application caractéristique, et alors tout se calcule explicitement grâce à la cohomologie de Gelfand-Fuchs, au moyen de polynômes universels analogues à ceux qui apparaissent dans les théorèmes d’indice classiques. Retour à la physique C’est là qu’a lieu, de façon inattendue, un retour à la physique, via les algèbres de Hopf. Les calculs des physiciens en théorie quantique des champs reposent sur des méthodes de développement perturbatifs où les termes sont des intégrales divergentes, qui nécessitent une renormalisation. Ces techniques de renormalisation font apparaître la combinatoire des diagrammes de Feynmann. Des formules empiriques dites de Bogoliubov-Parashiuk ramènent le calculs de diagrammes compliqués à des diagrammes plus simples. Or en 1998, le physicien Dirk Kreimer découvre que ces formules qui n’étaient a priori que de simples recettes, traduisent l’existence d’un objet mathé- matique, qui n’est autre qu’une algèbre de Hopf ou groupe quantique. C’est là que Connes rencontre Kreimer, et ils découvrent ensemble que l’algèbre de Hopf de Kreimer et celle de Connes-Moscovici sont essentiellement les mêmes. Autrement dit, ce sont les mêmes règles de symétries quantiques qui régissent d’une part les calculs de théorie quantique des champs (permettant de calculer, par des méthodes perturbatives, des quantités physiquement observables), d’autre part les calculs de géométrie non commutative (donnant explicitement l’indice d’opérateurs transversalement elliptiques sur des feuilletages). Un pas de plus a ensuite été franchi par Connes et Kreimer pour comprendre l’origine mathématique de cette algèbre de Hopf et son rôle dans le processus de renormalisation : ce processus n’est autre que la décomposition de Birkhoff, et ceci établit un lien direct et très simple avec le problème de Riemann-Hilbert. Ainsi, ce qui était au départ recette empirique, justifiée par l’expérience physique, est maintenant relié à un des grands problèmes des mathématiques, et non des moindres puisque le problème de Riemann-Hilbert est le 21ième de la liste des 23 problèmes proposés par David Hilbert au congrès international de Paris en 1900. Enfin, dans sa collaboration récente avec Matilde Marcolli, Alain Connes a trouvé la signification mathématique de cette correspondance de Riemann-Hilbert : cette dernière est reliée à la théorie de Galois motivique, introduite par Grothendieck. Elle fait apparaître un groupe de symétrie dont Pierre Cartier avait conjecturé l’existence sous le nom de « groupe de Galois cosmique » et qui est donc de nature arithmétique. Ainsi la géométrie non commutative relie la physique à la théorie des nombres. On entrevoit ainsi le lien entre ces deux univers mystérieux, celui des particules élémentaires et celui des nombres premiers, qui ont toujours fasciné Alain Connes. Pierre JULG Université d’Orléans et CNRS UMR 6628-MAPMO BP 6759 45067 Orléans Cedex 2 Pierre.Julg@univ-orleans.fr L’auteur remercie Alain Connes d’avoir relu son texte et de lui avoir suggéré quelques améliorations. Il remerci également Claire Anantharaman de lui avoir permis d’utiliser les posters qu’elle avait conçus en vue de la Journée Portes Ouvertes de l’Université d’Orléans. Résumé de Math Sup : déterminants I - Applications multilinéaires 1) Définition. Soient E1,..., En, F n + 1 espaces vectoriels . Soit f une application de E1 × ... × En dans F. f est n − linéaire ⇔ f est linéaire par rapport à chaque variable ⇔ ∀(ai)16i6n, ∀i ∈ J1, nK, l’application Ei → F xi 7→ f(a1, ..., ai−1, xi, ai+1, ..., an) est linéaire. Exemple. Dans R 3 euclidien orienté, f : R 3 × R 3 (x, y) 7→ x ∧ y est bilinéaire. Si E1 = ... = En = E et F = K, on obtient les formes n-linéaires sur E. 2/ Formes symétriques , antisymétriques , alternées. Définition. Soit f une forme n-linéaire sur E. 1) f est symétrique ⇔ ∀(x1, ..., xn) ∈ E n, ∀σ ∈ Sn, f(xσ(1) , ..., xσ(n)) = f(x1, ..., xn). 2) f est antisymétrique ⇔ ∀(x1, ..., xn) ∈ E n, ∀σ ∈ Sn, f(xσ(1) , ..., xσ(n)) = ε(σ)f(x1, ..., xn). 3) f est alternée ⇔ ∀(x1, ..., xn) ∈ E n, [(∃(i, j) ∈ J1, nK 2/ i 6= j et xi = xj) ⇒ f(x1, ..., xn) = 0]. Théorème. f est antisymétrique ⇔ ∀(x1, ..., xn) ∈ E n, ∀τ transposition de J1, nK, f(xτ(1) , ..., xτ(n)) = −f(x1, ..., xn). Démonstration. Soit σ ∈ Sn, on écrit σ = τ1 ◦ ... ◦ τk où les τi sont des transpositions et on sait que ε(σ) = (−1) k. Théorème. Soient E un K-espace vectoriel (K sous-corps de C) puis f une forme n-linéaire sur E. f alternée ⇔ f antisymétrique. Démonstration. ⇒ / Soit (x1, ..., xn) ∈ E n. Soient i 6= j puis τ = τi,j. 0 = f(x1, ..., xi + xj, ...xi + xj, ..., xn) = f(x1, ..., xi, ..., xi, ...xn) + f(x1, ..., xi, ..., xj, ...xn) + f(x1, ..., xj, ..., xi, ...xn) + f(x1, ..., xj, ..., xj, ...xn) = f(x1, ..., xi, ..., xj, ...xn) + f(x1, ..., xj, ..., xi, ...xn). Donc pour tout (x1, ..., xn) ∈ E n, pour toute transposition τ, f(xτ(1) , ..., xτ(n)) = −f(x1, ..., xn) et f est antisymétrique. ⇐ / Soit (x1, ..., xn) ∈ E n tel qu’il existe i 6= j tel que xi = xj = x. L’égalité f(x1, ..., xi, ..., xj, ...xn) = −f(x1, ..., xj, ..., xi, ...xn) s’écrit encore f(x1, ..., x, ..., x, ...xn) = −f(x1, ..., x, ..., x, ...xn) ou encore 2f(x1, ..., x, .., x, ...xn) = 0 ou enfin f(x1, ..., x, .., x, ...xn) = 0. II- Formes n-linéaires alternées sur un espace de dimension n . Définition de la forme déterminant dans une base 1) Théorème fondamental. Soit E un espace vectoriel de dimension finie n > 1. On note Λ∗ n(E) l’ensemble des formes n-linéaires alternées sur E. Théorème. 1) Λ∗ n(E) est un K-espace vectoriel de dimension 1. 2) Si B = (ei)16i6n est une base donnée de E, il existe une et une seule forme f n-linéaire alternée sur E telle que f(B) = 1. Définition. L’unique forme f n-linéaire alternée sur E telle que f(B) = 1 s’appelle la forme déterminant dans la base B et se note detB. Démonstration. Soit (x1, ..., xn) ∈ E n. Pour j ∈ J1, nK, posons xj = Xn i=1 xi,jei (où les xi,j sont dans K). Soit f une forme n-linéaire alternée sur E. En développant f(x1, ..., xn) par n-linéarité, on obtient une somme de n n termes du type xχ(1),1...xχ(n),nf eχ(1) , ..., eχ(n) où χ est une application quelconque de J1, nK dans lui-même. f est alternée et les termes correspondant aux applications χ telles que ∃i 6= j/ χ(i) = χ(j), sont nuls. Donc tous les termes pour lesquels χ n’est pas injective disparaissent. Maintenant, J1, nK étant un ensemble fini, χ est injective si et seulement si χ est une bijection de J1, nK sur lui-même ou encore une permuation de J1, nK. Il ne reste donc que les termes du type xσ(1),1...xσ(n),nf(eσ(1) , ..., eσ(n)) = ε(σ)xσ(1),1...xσ(n),nf(e1, ..., en) où σ est une élément quelconque de Sn. On a montré que nécessairement c Jean-Louis Rouget, 2009. Tous droits réservés. 1 http ://www.maths-france.fr ∀(x1, ..., xn) ∈ E n, f(x1, ..., xn) = f(e1, ..., en) X σ∈Sn ε(σ)xσ(1),1...xσ(n),n. Pour (x1, ..., xn) ∈ E n, posons ϕ(x1, ..., xn) = X σ∈Sn ε(σ)xσ(1),1...xσ(n),n. • ϕ est une forme n-linéaire sur E car linéaire par rapport à chaque variable. • ϕ est non nulle car ϕ(e1, ..., en) = X σ∈Sn ε(σ)δσ(1),1...δσ(n),n = 1. • ϕ est alternée. En effet, soit (x1, ..., xn) ∈ E n tel que xi = xj pour un certain couple (i, j) tel que i 6= j. Soit τ = τi,j. On sait que si An est l’ensemble des permutations paires, Anτ est l’ensemble des permutations impaires. Donc ϕ(x1, ..., xi, ..., xj, ..., xn) = X σ∈An xσ(1),1...xσ(n),n − X σ∈An xστ(1),1...xστ(n),n = X σ∈An xσ(1),1...xσ(i),i . . .xσ(j),j . . . xσ(n),n − X σ∈An xσ(1),1...xσ(j),i . . . xσ(i),j . . . xσ(n),n = X σ∈An xσ(1),1...xσ(i),i . . .xσ(j),j . . . xσ(n),n − X σ∈An xσ(1),1...xσ(j),j . . . xσ(i),i . . . xσ(n),n (car xi = xj) = 0. Finalement, Λ∗ n(E) = Vect(ϕ) avec ϕ 6= 0 et Λ∗ n(E) est un K-espace vectoriel de dimension 1. On a vu que ϕ(B) = 1 et que si f ∈ Λ∗ n(E), f = f(B)ϕ. Par suite, f(B) = 1 ⇔ f = ϕ. 2) Propriétés. Théorème. 1) detB(B) = 1. 2) detB′ = detB′ (B)detB. 3) detB(B′ ) × detB′ (B) = 1. 4) detB(B′ ) × detB′ (B′′) = detB(B′′). Démonstration. On applique : ∀f ∈ Λ∗ n(E), f = f(B)detB. 3) Applications. a) Théorème. Soit B une base de E de dimension finie n > 1 et B′ une famille de n vecteurs de E. B′ est une base de E si et seulement si detB(B′ ) 6= 0. Démonstration. Si B′ est une base, detB(B′ ) × detB′ (B) = 1 et en particulier detB(B′ ) 6= 0. Si B′ n’est pas une base, puisque card(B) = n, B′ est liée. Par suite, l’un des vecteurs de B′ est combinaison linéaire des autres vecteurs de B′ . Par n linéarité de detB et puisque detB est alternée, on a bien detB(B′ ) = 0. b) Orientation. Soient B et B′ deux bases de E 6= {0}. On définit la relation : « B′ a même orientation que B ⇔ detB(B′ ) > 0 ». La relation précédente est une relation d’équivalence à deux classes. On appelle arbitrairement l’une des deux classes, classe des bases directes et l’autre, classe des bases indirectes. L’espace E est alors orienté. III - Déterminant d’une matrice . Déterminant d’un endomorphisme . 1) Déterminant d’un endomorphisme. a) Définition. (Remarque. L’ordre de présentation 1) déterminant d’un endomorphisme 2) déterminant d’une matrice est un ordre qui permet de parvenir facilement à det(AB) = detA × detB). Soit f une forme n-linéaire alternée non nulle sur E (de dimension finie non nulle n) et soit u un endomorphisme de E. Soit fu définie par : ∀(x1, .., xn) ∈ E n, fu(x1, ..., xn) = f(u(x1), ..., u(xn)). fu est une forme n-linéaire alternée sur E. Donc il existe un scalaire λ tel que ∀(x1, ..., xn) ∈ E n, fu(x1, ..., xn) = λf(x1, ..., xn) ou encore telle que ∀(x1, ..., xn) ∈ E n, f(u(x1), ..., u(xn)) = λf(x1, ..., xn). Montrons que λ ne dépend que de u et pas de f. Soit g une forme n-linéaire alternée sur E. Puisque dimΛ∗ n(E) = 1, il existe un scalaire µ tel que g = µf. Donc pour (x1, ..., xn) ∈ E n g(u(x1), ..., u(xn)) = µf(u(x1), ..., u(xn)) = λµf(x1, ..., xn) = λg(x1, ..., xn). Définition. Soit u ∈ L (E). Le déterminant de u est le scalaire noté det(u) tel que ∀f ∈ Λ∗ n(E), ∀(x1, ..., xn) ∈ E n f(u(x1), ..., u(xn)) = (det(u)) × f(x1, ..., xn). c Jean-Louis Rouget, 2009. Tous droits réservés. 2 http ://www.maths-france.fr b) Propriétés. Théorème. Soit u ∈ L (E). 1) Pour toute base B de E, pour tout (x1, ..., xn) ∈ E n, detB(u(x1), ..., u(xn)) = (detu)detB(x1, ..., xn). 2) Pour toute base B de E, detB(u(B)) = det(u). Théorème. 1) det(Id) = 1. 2) ∀(u, v) ∈ (L (E))2 , det(u ◦ v) = (detu) × (detv). 3) ∀u ∈ L (E), (u ∈ GL(E) ⇔ detu 6= 0) et dans ce cas, det(u −1 ) = (detu) −1 . Démonstration. Soit B une base de E. 1) det(Id) = detB(B) = 1. 2) det(u ◦ v) = detB(u(v(B))) = (detu)detB(v(B)) = (detu)(detv). 3) Si u ∈ GL(E), (detu)(det(u −1 ) = det(u ◦ u −1 ) = det(Id) = 1 6= 0. Si detu 6= 0, detB(u(B)) = detu 6= 0 et u(B) est une base de E puis u ∈ GL(E). Théorème. ∀u ∈ L (E), ∀λ ∈ K, det(λu) = λ n(detu). Danger. En général, det(u + v) 6= detu + detv. 2) Déterminant d’une matrice carrée. a) Définition. Si A = (ai,j)16i,j6n ∈ Mn(K), le déterminant de A est le nombre detA = X σ∈Sn ε(σ)aσ(1),1...aσ(n),n. Notation. detA = |ai,j| 16i,j6n . b) Propriétés. Théorème. Soit u ∈ L (E) de matrice A dans une base B de E donnée. Alors det(u) = det(A). Théorème. Deux matrices semblables ont même déterminant Théorème. 1) L’application det : (Mn(K), ×) → (K, ×) A 7→ detA est un morphisme et l’application det : (GLn(K), ×) → (K∗ , ×) A 7→ detA est un morphisme de groupes. 2) ∀(A, B) ∈ (Mn(K))2 , det(AB) = (detA)(detB). 3) det(In) = 1. 4) ∀A ∈ Mn(K), [A ∈ GLn(K) ⇔ detA 6= 0] et dans ce cas det(A−1 ) = (detA) −1 . Le noyau du morphisme de groupes det : (GLn(K), ×) → (K∗ , ×) A 7→ detA c’est-à-dire l’ensemble des matrice carrées de déterminant 1 est un sous-groupe de (GLn(K), ×) noté SLn(K) (groupe spécial linéaire). Théorème. detA = det( tA). Démonstration. det( tA) = X σ∈Sn ε(σ)a1,σ(1) ...an,σ(n) . Soit σ un élément donné de Sn. Si on pose i1 = σ(1), ... , in = σ(n), alors σ −1 (i1) = 1, ..., σ −1 (in) = n. Le monôme a1,σ(1) ...an,σ(n) s’écrit aσ−1(i1),i1 ...aσ−1(in),in ou encore aσ−1(1),1 . . . aσ−1(n),n, après avoir remis dans l’ordre les n facteurs. Donc, det( tA) = X σ∈Sn ε(σ)a1,σ(1) ...an,σ(n) = X σ∈Sn ε(σ)aσ−1(1),1 . . . aσ−1(n),n = X σ∈Sn ε(σ −1 )aσ−1(1),1 . . . aσ−1(n),n = X σ′∈Sn ε(σ ′ )aσ′(1),1 . . . aσ′(n),n = det(A) car l’application σ 7→ σ −1 est une permutation de Sn (puisque application involutive de Sn dans lui-même). IV - Calculs de déterminants 1) Transposition. detA = det( tA) et donc toutes les règles portant sur les colonnes sont encore valables sur les lignes. 2) Matrices triangulaires. Le déterminant d’une matrice triangulaire est égal au produit de ses coefficients diagonaux (se montre directement en revenant à la définition). 3) Opérations élémentaires. a) ∀σ ∈ Sn, det(Cσ(1) , ..., Cσ(n)) = ε(σ)det(C1, ..., Cn). Quand on permute des colonnes, le déterminant est multiplié par la signature de la permutation. b) Si on ajoute à une colonne une combinaison linéaire des autres colonnes, le déterminant garde la même valeur. c Jean-Louis Rouget, 2009. Tous droits réservés. 3 http ://www.maths-france.fr c) det est n-linéaire et donc det(C1, ..., Ci+C ′ i , ..., Cn) = det(C1, ..., Ci, ..., Cn)+det(C1, ..., C′ i , ..., Cn) et det(C1, ..., λCi, ..., Cn) = λdet(C1, ..., Cn). Danger. det(A + B) 6= detA + detB en général et det(λA) = det(λC1, ..., λCi, ..., λCn) = λ ndet(C1, ..., Cn) = λ ndetA. 4) Calculs par blocs. Théorème. Si les Ai sont des matrices carrées, det   A1 × . . . × 0 . . . . . . . . . . . . . . . . . . × 0 . . . 0 Ap   = (detA1)(detA2)...(detAp). Démonstration. On le montre pour p = 2 puis on termine par récurrence. On montre que det A C 0 B = (detA)(detB) où A ∈ Mn(K) et B ∈ Mp(K). Soit ϕ : Mn(K) → K X 7→ det X C 0 B . ϕ « est » une application n-linéaire alternée des colonnes de X. Donc, il existe une constante k (indépendante de X) telle que ∀X ∈ Mn(K), ϕ(X) = k(detX) avec k = kdet(In) = ϕ(In) ce qui montre que ∀X ∈ Mn(K), det X C 0 B = detX × det In C 0 B . De même, l’application ψ : Mn(K) → K X 7→ det In C 0 Y « est » une application p-linéaire alternée des lignes de Y. Donc, il existe une constante k ′ (indépendante de Y) telle que ∀Y ∈ Mp(K), ψ(Y) = k(detY) avec k ′ = k ′detIp = ψ(Ip). Ceci montre que ∀X ∈ Mn(K), ∀Y ∈ Mp(K), det X C 0 Y = detX×detY ×det In C 0 Ip . Enfin le déterminant d’une matrice triangulaire est égale au produit de ses coefficients diagonaux et donc det In C 0 Ip = 1. 5) Développement suivant une ligne ou une colonne. Théorème. Soient mi,j le mineur de ai,j et Ai,j = (−1) i+jmi,j = cofacteur de ai,j. Alors, ∀(i, j) ∈ J1, nK 2 , detA = Xn k=1 ai,kAi,k (développement suivant la ligne i) = Xn k=1 ak,jAk,j (développement suivant la colonne j) Démonstration. Il suffit de démontrer la formule de développement suivant une colonne car detA = det( tA). Ensuite, il suffit de démontrer la formule de développement suivant la première colonne car alors, si on veut développer suivant la colonne j, on effectue la permutation des colonnes Cj → C1 → C2... → Cj−1 dont la signature est (−1) j−1 (signature d’un cycle de longueur j), puis en développant suivant la première colonne, on obtient detA = (−1) j−1 Xn k=1 ak,j(−1) k+1mk,j = Xn k=1 ak,jAk,j. Il reste à démontrer la formule de développement suivant la première colonne. C1 est somme de n colonnes du type   0 . . . 0 ai,1 0 . . . 0   et par n-linéarité du déterminant, detA = Xn i=1 detAi où detAi = 0 a1,2 . . . . . . a1,n . . . . . . . . . 0 ai−1,2 . . . . . . ai−1,n ai,1 ai,2 . . . . . . ai,n 0 ai+1,2 . . . . . . ai+1,n . . . . . . . . . 0 an,2 . . . . . . an,n . c Jean-Louis Rouget, 2009. Tous droits réservés. 4 http ://www.maths-france.fr Si i = 1, un calcul de déterminant par blocs fournit detA1 = a1,1A1,1. Si i > 2, on passe Li en L1, L1 en L2,..., Li−1 en Li. On obtient detAi = (−1) i−1 ai,1 ai,2 . . . . . . ai,n 0 a1,2 . . . . . . a1,n . . . . . . . . . 0 ai−1,2 . . . . . . ai−1,n 0 ai+1,2 . . . . . . ai+1,n . . . . . . . . . 0 an,2 . . . . . . an,n = Ai,1 (calcul par blocs) et finalement detA = Xn i=1 ai,1Ai,1. V - Comatrice. Inverse d’une matrice . La comatrice de la matrice carrée A de format n est la matrice, notée comA, dont le coefficient ligne i, colonne j, est le cofacteur de l’élément ai,j de A, c’est-à-dire si A = (ai,j)16i,j6n, alors com(A) = (Ai,j)16i,j6n. La transposée de la comatrice de A est appelée matrice complémentaire de A et est notée A˜ . Théorème. 1) ∀A ∈ Mn(K), AA˜ = AA˜ = (detA)In ou encore A( t comA) = (t comA)A = (detA)In. 2) ∀A ∈ Mn(K), (A ∈ GLn(K) ⇔ detA 6= 0) et dans ce cas, A−1 = 1 detA ( t comA). Démonstration. Le coefficient ligne i, colonne j de At comA vaut Xn k=1 ai,kAj,k. • Si i = j, cette expression n’est autre que le développement de detA suivant sa i-ème ligne et vaut donc detA. • Si i 6= j, Xn k=1 ai,kAj,k est le développement suivant la ligne j du déterminant déduit de detA en remplaçant la ligne j de detA par sa ligne i (et en ne modifant pas sa ligne i). Cette expression est donc nulle puisque égale à un déterminant ayant deux lignes identiques. Ensuite , il est clair que com( tA) = t (comA) (à partir de la définition de comA) et donc AA˜ = (t comA)A = com( tA) t ( tA) = t ( tAt (com( tA)) = t ((det( tA)In) = (detA)In. c Jean-Louis Rouget, 2009. Tous droits réservés. 5 http ://www.maths-france.fr

Voir également :

MATHEMATIQUESFINANCI ..> 29-Sep-2011 14:28 12K STATISTIQUEDESCRIPTI ..> 08-Oct-2011 11:48 15K approxautosemblabler ..> 29-Sep-2011 14:52 3.1k baseslogrecurrences.htm 29-Sep- 2011 14:57 21K calculfinanciervaleu ..> 29-Sep-2011 14:31 100K elementsmathfi.htm 29-Sep-2011 14:40 15K financementetemprunt ..> 29-Sep-2011 14:39 10K infocontextuelleweb.htm 29-Sep- 2011 15:02 81K informationpersonnal ..> 29-Sep-2011 15:04 156K interrogationflexibl ..> 29-Sep-2011 15:09 13K martingalepourlafina ..> 08-Oct-2011 11:48 47K mathematiquefinancie ..> 29 - Sep-2011 17:54 45K mathfiPGM110.htm 29-Sep-2011 14:44 10K mathphysiquestatisti ..> 29-Sep-2011 14:55 22K methodestatistiquege ..> 18-Oct-2011 17:05 116K metiersmathematiques ..> 29 -Sep-2011 14:51 114K probabilitesimulatio ..> 08-Oct-2011 11:48 96K rechercheinfobibliot ..> 29-Sep-2011 15:12 30K securiteinformatique ..> 29-Sep-2011 15:00 20K statistiqueelementai .. > 08-Oct-2011 11:48 101K

Pattern Recognition 40 (2007) 3277 – 3291 www.elsevier.com/locate/pr Signi?cant edges in the case of non-stationary Gaussian noise I. Abrahama , R. Abrahamb , A. Desolneux c,* , S. Li-Thiao-Te d a CEA/DIF, 91680 Bruyères le Chatel, France b Laboratoire MAPMO, Fédération Denis Poisson, Université d’Orléans, B.P. 6759, 45067 Orléans cedex 2, France c Laboratoire MAP5, Université René Descartes, 45 rue des Saints-Pères, 75270 Paris cedex 06, France d Laboratoire CMLA, ENS Cachan, 61 avenue du Président Wilson, 94235 Cachan cedex, France Received 8 June 2006; received in revised form 27 February 2007; accepted 28 February 2007 Abstract In this paper, we propose an edge detection technique based on some local smoothing of the image followed by a statistical hypothesis testing on the gradient. An edge point being de?ned as a zero-crossing of the Laplacian, it is said to be a signi?cant edge point if the gradient at this point is larger than a threshold s () de?ned by: if the image I is pure noise, then the probability of ?I (x )s () conditionally on I (x ) = 0 is less than . In other words, a signi?cant edge is an edge which has a very low probability to be there because of noise. We will show that the threshold s () can be explicitly computed in the case of a stationary Gaussian noise. In the images we are interested in, which are obtained by tomographic reconstruction from a radiograph, this method fails since the Gaussian noise is not stationary anymore. Nevertheless, we are still able to give the law of the gradient conditionally on the zero-crossing of the Laplacian, and thus compute the threshold s (). We will end this paper with some experiments and compare the results with those obtained with other edge detection methods. 2007 Pattern Recognition Society. Published by Elsevier Ltd. All rights reserved. Keywords: Edge detection; Signi?cant edges; Inverse problem; Statistical hypothesis testing 1. Introduction This work stems from analysis of physical experiments where a test object is submitted to a shock. More precisely, we surround the object with explosives and we monitor the shape of the hull as it ?rst collapses onto itself and then expands. The shape of inner hole(s) at the very beginning of the re-expansion phase is of particular interest. We use radiography to determine those shapes. Standard tomography cannot be carried out because the inverse Radon transform requires many radiographs (taken from different angles), whereas only one radiograph can be acquired in the lapse of the expansion phase. That is why the hulls are supposed to be radially symmetric at the beginning of the experiment and during the explosion. In that case, * Corresponding author. Tel.: +33 1 44 55 35 26. E-mail addresses: isabelle.abraham@cea.fr (I. Abraham), romain.abraham@univ-orleans.fr (R. Abraham), agnes.desolneux@math-info.univ-paris5.fr (A. Desolneux), lithiaote@cmla.ens-cachan.fr (S. Li-Thiao-Te). 0031-3203/$30.00 2007 Pattern Recognition Society. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.patcog.2007.02.015 a single radiograph is enough to perform the tomographic reconstruction (see [1,2], and Section 4.1). Thus, from now on, we will only consider radially symmetric objects. To describe such an object, it is suf?cient to provide the densities of the materials on a slice of the object that contains the symmetry axis. An example of studied object is given in Fig. 1. The ?rst step of the physical experiment is to take a radiography of this object (see Fig. 2(a)), then a tomographic reconstruction is computed (Fig. 2(b)) and ?nally an edge detection is performed (Fig. 2(c)). As we can see in Fig. 2(c), many detected edges do not correspond to real features. This is due to the high level of noise. Our goal is to perform an automated selection of the “true” edges. For that purpose, the edge detector will not be changed but we will compute also other signi?cant features that will allow us to select the “true” edges. Let us stress that our goal is not to develop a new edge detector but to give a method of selection of edges from the map given by an edge detector. We focus here on an edge detector which is certainly not optimal but where the computations are easy. However, we think that our method can be applied to other edge detectors3278 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 such as the Canny detector [3]. In his paper, Canny points out the fact that “it is very dif?cult to set a threshold so that there is small probability of marking noise edges while retaining high sensitivity.” And the solution he proposes to solve this problem is the use of hysteresis. Here, we will use the zero-crossings of the Laplacian to de?ne edge points. We will see that in this case, we can compute a threshold so that there is a small probability of selecting edges due to noise. The ideas used here mainly come from previous work of Desolneux et al. [4]. Informally speaking, they de?ne the notion of signi?cant edges by computing the probability that some edge-related events appear in an image of pure noise. When this probability is small enough, the edge is probably a feature of the image and not due to noise. Unfortunately, that method assumes that the noise is stationary (i.e., the law of the noise does not depend on the location in the image) which, as is easily seen in Fig. 2(b), is not the case in our study because of the tomographic inversion (see Section 4.5.1 for some examples of results obtained with their method). Moreover, their study is quite general and, apart from the stationarity, no assumption is made on the noise. In our case, as we deal with speci?c images, the noise is well-known and some statistical models can be used. Indeed, we may suppose that the noise on the radiograph (Fig. 2(a)) is a Gaussian white noise with mean zero and with a variance that can easily be estimated. After tomographic inversion, we still obtain Gaussian noise because that operation is linear, but the noise process is correlated and non-stationary. The edge detector we will use here consists in estimating the Laplacian at each point, and edge points are then de?ned as the Fig. 1. Slice of a studied object. Fig. 2. (a) Radiograph of the object of Fig. 1. (b) Tomographic reconstruction. (c) Edge detection on the reconstructed image. zero-crossings of the Laplacian. As already said, we will compute additional features that discriminate the signi?cant edges. The easiest feature to compute is a kind of contrast measurement C based on a gradient estimate. To be more precise, we consider an image I of pure noise (that is a realization of our model of noise after tomographic reconstruction), estimate the gradient and the Laplacian of I at a point (u, v) (with an abuse of notation, we will denote by ?I (u, v) and I (u, v) these estimates and by C(u, v) the contrast value) and we compute, for a ?xed > 0, the smallest value s () such that P(C(u, v)s ()|I (u, v) = 0), (1) where P(B|A) denotes the conditional probability of B given the event A. Then, we perform an edge detection on an image f of interest (where we also estimate ?f and f by the same method) and we keep the points (u, v) of the studied image f that satisfy: • f (u, v) = 0 (an edge is present at point (u, v)). • C(u, v)s () (this edge is signi?cant). From a statistical point of view, this consists in performing a hypothesis test. We consider a point (u, v) where an edge takes place (f (u, v) = 0) and we test the null hypothesis “the zerocrossing of the Laplacian is due to the noise.” The level of the test is arbitrarily chosen and related to the number of false detections allowed. It will be set to = 10-5 hereafter (in a 512 × 512 image, this corresponds to 512 × 512 × 10-5 2.6 false detections on average). Let us mention that the threshold value s () varies slowly with respect to . For instance, in the case of a white noise (see Section 3), the threshold value can be computed explicitly and is proportional to v - ln . This implies that when we increase the speci?city (less false detections), the threshold increases slowly and the method retains much of its sensitivity. When the null hypothesis is rejected, the edge is retained as it comes from a “true” feature of the image, whereas when the null hypothesis is accepted, the zero-crossing of the Laplacian may come from the noise and the edge is not signi?cant. Let us mention at this point that such statistical approaches have already been used for edge detection in [5–7]. In [8], Leclerc and Zucker develop a technique to capture discontinuities de?ned (in one dimension) as points such that the estimated left- and right-hand limits are statistically signi?cantlyI. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3279 different. The main dif?culty in their approach is that they need to have good estimates of the parameters of the ?tted curves (generally polynomials) on half-neighborhoods of points in order to apply some statistical test. To solve this problem, they propose a procedure to eliminate points where the estimates are not good, which are points whose half-neighborhoods overlap (unknown) discontinuities. Other methods usually employ estimates of the gradient based on ?nite differences which fail in our case because these are not robust enough to high levels of noise. Moreover, the noise is in most cases stationary. Let us also cite Ref. [9] where the authors have modi?ed the method of Ref. [4] to take into account the non-stationarity of some images, by a local noise estimate. Their work is still general and does not make any assumption on the noise structure. As we deal with speci?c experiments, the noise is always the same and well-known and we can take proper advantage of this knowledge. The problem of ?nding an optimal threshold is also a key point for image thresholding methods [10] and particularly for the binarization of document images [11]. In some cases, the thresholding is not stationary: this happens for instance in images with non-uniform illumination. Our framework here is different since we deal with a well-known noise, and the question we address is: given an edge detection method, what is the probability that an edge point is due to noise? The paper is organized as follows: in Section 2, we present the edge detector based on the estimates of the gradient and of the Laplacian. Then, in Section 3, our method is presented in the case of a Gaussian white noise. Of course, this does not correspond to our case but the computations are easier and demonstrate the performance of this framework. In Section 4.1, we will ?rst describe the tomographic inversion and the operators involved, and then describe the noise model we have to deal with. We will then apply the signi?cant edges detection method in the framework of this non-stationary noise. We will end the section with some experiments and comparisons with other methods. 2. Estimating the gradient and the Laplacian In this section, we introduce the method for edge detection. We consider that an image is a real-valued function (u, v) ? f (u, v) of two continuous real parameters u and v. Then, we say that there exists an edge at point (u, v) if the Laplacian of f is zero at this point. Moreover, the computation of the contrast function C will be based on the gradient of f (see the end of this section for the choice of this function). As the images are often very noisy, these derivatives cannot be estimated by usual ?nite differences. The method used here, sometimes known as Savitsky–Golay smoothing [12], consists in locally approximating the image by a polynomial. The derivatives of the polynomial are then identi?ed with those of the image. 2.1. An optimization problem Let (u, v) denote the point where we want to compute the ?rst- and second-order derivatives of the image f. We chose two parameters: d which is the maximum degree of the approximating polynomial and r which is the radius of the ball on which we perform the approximation. We denote by Br (u, v) the ball of radius r centered at point (u, v). We will simply write Br when the center of the ball is the origin (0, 0) of R 2 . We are then looking for a polynomial P of degree less than d such that E(P ) = Br (f (u + x , v + y ) - P (x , y )) 2 dx dy (2) is minimal among all polynomials of degree less than d. In other words, we are looking for the best approximation of f by a polynomial of degree less than d on the ball Br (u, v) in the sense of the L 2 -norm. This is an optimization problem where the unknowns are the coef?cients of the polynomial. As the problem is convex, there is a unique solution (given by the orthogonal projection of f on the space of polynomials of degree less than d) which is easily computed by solving the set of equations E ai = 0, where the ai ’s denote the coef?cients of the polynomial. Role of the ball radius: Two parameters are arbitrarily chosen in this method. The ?rst one is the ball radius the r. The larger the r is, the more effective the smoothing is. The in?uence of the noise is therefore attenuated with a large r but the location of the edge is then less precise. We must consequently make a balance between noise smoothing and edge detection accuracy. For instance, if we have a small level of noise or if the edges are very complicated (with high curvature), we must choose a small value for r. Role of the polynomial degree: The second parameter is the polynomial degree. Here again a large value of d gives a better approximation but does not smooth the noise enough. In fact, as we are, in a ?rst step, interested in the points where the Laplacian is zero, it appears that a second-order polynomial is enough. Of course, the estimate of the ?rst-order derivatives with a polynomial of degree 2 is not very good and highly depends on the size of the window Br . But we will see that this drawback can be useful for the choice of a contrast function. Choice of the L 2 -norm: Approximating a function in an L 2 -sense, although quite usual, is not always accurate as some oscillations may appear or local bad approximations are allowed (this usually occurs near the boundary of the domain). However, as we will only look for polynomials of degree 2, the oscillations do not appear and the L 2 - norm gives a good enough approximation together with easy computations. In what follows, the approximation is made with a polynomial of degree d =2, and the ?rst- and second-order derivatives of the image are identi?ed with those of the approximating polynomial.3280 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 2.2. Computations with a second-order polynomial Let us ?rst introduce some notations. In the following, we will set ?i, j ? N, bij (r) = Br x i y j dx dy. As the ball Br is symmetric, we have that bij (r) = 0 if i or j is odd and that bij (r) = bj i (r) for all i,j. In order to have simple expressions, we will also set b(r) = b20(r), (r) = - 2b20(r) b00(r) and (r) = 1 2 b40(r) + b22(r) - 2b 2 20 (r) b00(r) . Lemma 1. The gradient and the Laplacian of the polynomial P of degree 2 which is the best approximation of f on the ball Br (u, v) for the L 2 -norm, being, respectively, denoted by ?r f (u, v) = (( r f /x)(u, v), ( r f /y)(u, v)) and r f (u, v), are given by r f x (u, v) = 1 b(r) Br xf (u + x , v + y ) dx dy, r f y (u, v) = 1 b(r) Br yf (u + x , v + y ) dx dy, r f (u, v) = 1 (r) Br f (u + x , v + y ) × ((r) + x 2 + y 2 ) dx dy. Proof. We consider a polynomial of degree 2 which we write as P (x , y ) = ax x x 2 + ayy y 2 + axy xy + ax x + ay y + a0. The equations obtained by writing ?E(P ) = 0, where E(P ) is given by Eq. (2), are ? ? ? b40(r)ax x + b22(r)ayy + b20(r)a0 = Br x 2 f (u + x , v + y ) dx dy , b22(r)ax x + b40(r)ayy + b20(r)a0 = Br y 2 f (u + x , v + y ) dx dy , b22(r)axy = Br xyf (u + x , y + v) dx dy , b20(r)ax = Br xf (u + x , y + v) dx dy , b20(r)ay = Br yf (u + x , v + y ) dx dy , b20(r)ax x + b20(r)ayy + b00(r)a0 = Br f (u + x , v + y ) dx dy . Fig. 3. Object with an inhomogeneous material. We then obtain the following estimates for the derivatives: P x (0, 0) = ax = 1 b20(r) Br xf (u + x , v + y ) dx dy, P y (0, 0) = ay = 1 b20(r) Br yf (u + x , v + y ) dx dy, P (0, 0) = 2(ax x + ayy ) = 2 b40 + b22 - 2b 2 20 /b00 Br f (u + x , v + y ) × - 2b20 b00 + x 2 + y 2 dx dy . 2.3. Choice of the contrast function We would like to use a contrast function based on the estimates of the ?rst and second derivatives of the image f obtained in the previous section. The simplest contrast function we can choose is the norm of the gradient: C1(u, v) = ?r f (u, v). Indeed, the value of this norm tells how sharp the edge is. This contrast function performs reasonably well and will be used when the images we deal with are piecewise constant. However, in many cases, the objects we handle are not homogeneous and their images contain some slopes (see Fig. 3). In this case, the gradient norm is not a good contrast function. Indeed, let us consider an image with a constant slope with some noise (see Fig. 4). We would like to say that no edge is signi?cant in that case. However, the value of the gradient norm (which will be close to the value of the slope) will always be greater than the threshold value s when the noise level is small. In the latter case, we take advantage of the dependence of the ?rst-order derivative estimates with respect to the ball radius. Indeed, the estimates of the gradient in the case of the constant slope in Fig. 4 will not depend on the size of the window (see Fig. 4) whereas, when an edge (a discontinuity) occurs, the estimates do depend on that radius (see Fig. 5). So, we can useI. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3281 Fig. 4. A noisy constant slope: the gradient of the approximating polynomial does not depend on the radius r. Fig. 5. An edge on a noisy slope: the approximating polynomials with two different values of the radius r. as a contrast function the function C2(u, v) = ?r1 f (u, v) - ?r2 f (u, v), where r1 < r2 and ?r f denotes the value of the gradient estimate with a ball of radius r. 3. Signi?cant edges in the case of a Gaussian white noise 3.1. White noise and Wiener integral We recall here the de?nition and the main properties of a white noise in a continuous setting and of the Wiener integral. We refer to Refs. [13–15] for more on white noise and the Wiener integral. De?nition 1. A Gaussian white noise on R 2 of variance 2 is a random function W de?ned on the Borel sets A of R 2 of ?nite Lebesgue measure (denoted by |A|) such that • W(A) is a Gaussian random variable (r.v.) with mean 0 and variance 2 |A|. • If A1 n A2 = Ø, the r.v. W (A1) and W (A2) are independent and W (A1 ? A2) = W (A1) + W (A2). Such a function W exists but is not a true measure since the two-parameter process B(s , t ) := W ((0, s] × (0, t ]) (with the convention that (0, x] is [x , 0) when x is negative) is of unbounded total variation. This process is usually called the Brownian sheet. Nevertheless we can de?ne the so-called Wiener integral f dW for every function f in L 2 (R 2 ). We can also de?ne the derivatives of the Brownian sheet in the sense of Schwartz distributions (although the Brownian sheet is nowhere differentiable). Thus, we de?ne B (s, t ) ? = 2 B(s , t ) st and we have f dW = R 2 f (u, v)B(u, v) ? du dv a.s. for every function f in the Schwartz space. With a slight abuse of notations, we call B? a Gaussian white noise and we always denote by R 2 f (u, v)B(u, v) ? du dv the Wiener integral with respect to this white noise, for every function f ? L 2 . The main properties of this integral are: • For every f, the r.v. R 2 f (u, v)B(u, v) ? du dv is a Gaussian r.v. with mean 0 and variance 2 R 2 f (u, v) 2 du dv. • For every f, g, the random vector R 2 f (u, v)B(u, v) ? du dv, R 2 g(u, v)B(u, v) ? du dv is Gaussian with cross-covariance 2 R 2 f (u, v)g(u, v) du dv. We will use these properties to compute the laws of ?I and I . 3.2. Laws of the gradient and of the Laplacian We suppose here that our noise is a Gaussian white noise, of variance 2 . As we have already said, this case is not the one we are interested in and our method is probably out-performed by other standard methods in that case. The goal of this section is to present our method in a simple case where the computations are easy and can be carried out in a continuous setting. We will only focus here on the case of piecewise constant objects and therefore we will use the contrast function C1. Lemma 2. If the image I is a Gaussian white noise of variance 2 , then r I x , r I y , r I 3282 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 is a Gaussian random vector with mean zero and covariance matrix ? ? 2 b(r) 0 0 0 2 b(r) 0 0 0 V (r, ) ? ? , where V (r, ) = 2 2 (r) Br ((r) + x 2 + y 2 ) 2 dx dy. Proof. We compute the laws of the approximate derivatives of I when I = B? . We recall that these derivatives are given by r I x (u, v) = 1 b(r) Br xB (u ? + x , v + y ) dx dy, r I y (u, v) = 1 b(r) Br yB (u ? + x , v + y ) dx dy, r I (u, v) = 1 (r) Br B (u ? + x , v + y ) × ((r) + x 2 + y 2 ) dx dy. Because of the stationarity of B? , they have the same law as r I x (0, 0) = 1 b(r) Br xB (x , y ) ? dx dy, r I y (0, 0) = 1 b(r) Br yB (x , y ) ? dx dy, r I (0, 0) = 1 (r) Br B (x , y )( ? (r) + x 2 + y 2 ) dx dy. As we deal with Wiener integrals, we deduce that the random vector r I x , r I y , r I is a Gaussian random vector with mean zero. To compute its covariance matrix, let us recall that, if X and Y are random variables de?ned by X = Br h1(x, y)B (x , y ) ? dx dy, Y = Br h2(x, y)B (x , y ) ? dx dy, then we have Cov(X, Y ) = 2 Br h1(x, y)h2(x, y) dx dy. Consequently, we have for instance: Cov r I x , r I y = 2 b 2 (r) Br xy dx dy = 0. By some analogous calculations, we ?nally get the following covariance matrix for our Gaussian random vector: ? ? 2 b(r) 0 0 0 2 b(r) 0 0 0 V (r, ) ? ? , where V (r, ) = 2 2 (r) Br ((r) + x 2 + y 2 ) 2 dx dy . Let us recall that two uncorrelated components of a Gaussian random vector are independent random variables (which is not the case in general). Therefore, we immediately have the following properties: • The random variable ?r I 2 is the sum of two squared independent Gaussian random variables which have the same variance. It is therefore distributed as a 2 -law. More precisely, its law is 2 b(r) 2 (2), where 2 (2) denotes a 2 -law with two degrees of freedom. • The random variable r I is a Gaussian random variable with mean zero and variance V (r, ). • The random variables ?r I 2 and r I are independent. 3.3. Computation of the threshold Proposition 1. Let I be a Gaussian white noise and let s () be the threshold value such that P(?r Is ()|r I = 0). Then, s () = - 2 2 b(r) ln . Proof. To begin with, as the random variables ?r I 2 and r I are independent, we can forget the conditioning and only compute P(?r Is ()) = P(?r I 2 s () 2 ). As a consequence of Lemma 2, we have that the law of ?r I 2 is ( 2 /b(r)) 2 (2). Now, since the density of a 2 (2) law is the one of a ( 1 2 , 1) law, we have that the law of ?r I 2 is given by P(?r I 2 s 2 ) = +8 (b(r)/2 )s 2 1 2 e -t /2 dt = exp - b(r)s 2 2 2 . This ?nally leads to the announced threshold value s (). I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3283 Fig. 6. (a) The noisy image ( = 0.2). (b) The zero-crossings of the Laplacian with the contrast function C1 visualized in gray level. (c) The extracted signi?cant edges ( = 10-5 ). Fig. 7. (a) The noisy image ( = 0.4). (b) The zero-crossings of the Laplacian with the contrast function C1 visualized in grey-level. (c) The extracted signi?cant edges ( = 10-5 ). Remark. Choice of the ball radius. Conversely, the result of Proposition 1 can also be used in order to ?x the ball radius using some prior knowledge on the edges. Indeed, if we know the signal-to-noise ratio (denoted by SNR) and if we suppose that the edges are not too complicated so that the image can be locally considered as the characteristic function of a half-plane, we can compute the estimated gradient for a given ball radius a, say g(a) (for a step of height 1): g(a) = Ba x 1x 0 dx dy. Then, if we choose the ball radius a such that exp(- 1 2 b(a)SNR 2 g(a) 2 ) the edges are signi?cant. As we mentioned it before, the ball radius must be chosen close enough to this critical value in order not to smooth too much the edges. 3.4. Experiments We consider the piecewise constant object of Fig. 1 with some additive Gaussian white noise. Note that the noisy images here are obtained simply by summing the image of Fig. 1 and the pure white Gaussian noise image, and any tomographic reconstruction procedure is not applied unlike the experiments described in Section 4.5. The densities of the different materials of the object of Fig. 1 are: • 1 for the outer material. • 0.8 for the material inside the circle. • 0.3 for the other inner material. The standard deviation of the Gaussian noise is = 0.2 in the experiments of Fig. 6 and is = 0.4 in the experiments of Fig. 7. Both images have the same size 512 × 512 pixels. The experiments have been carried out with a ball of radius r = 12 pixels. The different images of Figs. 6 and 7, are, respectively: • (a) The noisy image. • (b) The zero-crossings of the Laplacian with the contrast function C1 visualized in gray level (the white color corresponds to high value for the contrast function C1). • (c) The extracted signi?cant edges ( = 10-5 ). In the case of a SNR large enough (Fig. 6), all the edges are well detected and the “false” edges are removed. However, the edges which have a high curvature are smoothed by our algorithm. This drawback is all the more important when the ball radius r is large (the in?uence of the value of this radius will be studied in the experiments of the next section).3284 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 When the noise level is rather large (Fig. 7), some edges of the image cannot be extracted from the noise (this happens when the contrast associated with this edge is close to the noise level). Estimating the noise variance: If the variance of the Gaussian noise is not readily available, standard methods can be used in order to estimate it directly from the image. For instance, as it is assumed that the noise is a centered Gaussian white noise, the image f we are studying can be written as f (x , y ) = f0(x, y) + (x, y), where f0 is a deterministic function and is the additive Gaussian noise (with mean 0 and unknown variance 2 ). Let h be the pixel size and let us consider the random variable V (x , y ) = f (x + h, y) + f (x - h, y) + f (x , y + h) + f (x , y - y ) - 4f (x , y ). Then, V (x , y ) = h 2 f0(x, y) + o(h 2 ) + (x + h, y) + (x - h, y) + (x, y + h) + (x, y - h) - 4 (x, y). Therefore, if we suppose h small, we have E[V 2 (x, y)] h 4 f0(x, y) 2 + 20 2 . Hence, as f0 is usually small, the quantity 2 = 1 20N2 1 x N 1 y N V 2 (x, y) is a good estimator of the noise variance. In any case, the variance is always over-estimated by this method, and consequently, the signi?cant edges detected when using 2 are also signi?cant for the variance 2 . 4. Signi?cant edges in the case of Gaussian white noise on a radiograph 4.1. Tomography Let us now turn to the more realistic case we are interested in. As we mentioned it in Introduction, we ?rst make a radiography of an object. Tomography is the inverse problem associated with reconstructing the initial object from its radiograph. This is now a well-known problem as it is the key tool in medical scanner imagery (or other medical imaging systems). To begin with, let us describe what radiography is from a mathematical point of view. The studied object is exposed to X-rays that traverse it and some of the X-photons are absorbed. In the output, we observe the quantity of X-photons that have not been absorbed by the material, and we thus measure in some sense the “mass” of material the ray went through. More precisely, if the object is described by its density (which is y x u v axis of symmetry Fig. 8. Radiography of a radially symmetric object. a function of the space coordinates), what can be measured at any point of the receptor is ray d, where “ray” means the straight line that goes from the source to the studied point of the receptor (we suppose that the X-ray source is just a point, which implies that the previous line is unique). We also assume that the X-ray source is far away from the object so that the rays are assumed to be all parallel. Then, to reconstruct an object from its radiographs, we must rotate the object and take a radiograph for every angle ? [0, ). This leads to the so-called Radon transform of the object, which is known to be invertible. This is the principle of the medical scanner. In our case, as the object is radially symmetric, if we rotate around its symmetry axis, all the radiographs are exactly the same. Consequently, a single radiograph of such an object is enough to perform the tomographic reconstruction. Indeed, if f (x , y ) denotes the density along a slice that contains the symmetry axis (see Figs. 1 and 8), then a radiograph of this object is given by g(u, v) = 2 +8 |u| f (x , v ) x v x 2 - u 2 dx. This is a linear transform and we will denote it hereafter by g = Hf . As we already said, this linear operator H is invertible and we in fact know explicitly its inverse on the space of continuously differentiable functions g: f (x , y ) = (H -1 g)(x, y) = - 1 +8 x 1 v u 2 - x 2 g u (u, y) du. (The proof of this formula and more details about the operators H and H -1 can be found in the book of Bracewell [16, pp. 244–250]).I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3285 Our assumption on the noise is that it is an additive Gaussian white noise on the radiograph (i.e., on g). But what we want is to study the object given by f. So we must transform the white noise by the operator H -1 . Unfortunately, because of the singularity of the integral at x = 0, we cannot apply the operator H -1 to a white noise B? , even in an L 2 -sense. Therefore, we will work in a discrete framework: the images f and g are naturally discretized (as they are numerical images). This leads to a discretization of the operator H, which we will still denote by H and which now may be viewed as a matrix. The discretization is made in such a way that the symmetry axis (x = 0) is settled between two pixels so that the previous singularity does not appear. This matrix is then invertible and we denote by H-1 its inverse which we can make now operate on a discrete Gaussian white noise. To give an idea of the ill-conditionedness of the matrix H, we can compute its condition number. For an image of size 512 × 512, we ?nd that this condition number is 377, which shows how H-1 ampli?es noise. 4.2. Law of the noise on the tomographic reconstruction Let us consider a ?eld = ( i,j )1 i p,1 j n of independent identically distributed (i.i.d.) random Gaussian variables with mean 0 and variance 2 . Let us de?ne I = (Ii,j ) the random ?eld obtained after tomographic reconstruction, i.e., after making H -1 operate on = ( i,j ). In fact, as the X-rays are supposed to be parallel, the reconstruction can be made line by line independently and therefore, if we consider the row vectors i = ( i,1 , . . . , i,n ) and I i = (Ii,1, . . . , Ii,n) then, if we denote M = (H t ) -1 (where H t is the transpose of H), we have I i = iM. (Notice that M is independent of i and that its size is n × n.) Consequently, the law of I is characterized by the following properties: • I = (Ii,j ) is a Gaussian random ?eld. • For i = k, I i and I k are independent. • For each i, the vector I i is a Gaussian random vector of mean 0 and covariance matrix = 2 Mt M, where Mt denotes the transpose of M. 4.3. Laws of the gradient and of the Laplacian The expressions obtained in Section 2 for the gradient and for the Laplacian of an image in a continuous setting are easily translated in the discrete framework we now deal with. Indeed, we have r I x (u, v) = 1 b(r) (i,j )?Br j I u+i,v+j , r I y (u, v) = 1 b(r) (i,j )?Br i I u+i,v+j , r I (u, v) = 1 (r) (i,j )?Br ((r) + i 2 + j 2 )Iu+i,v+j , where Br now denotes the discrete ball of radius r, i.e., Br = {(i, j ), i 2 + j 2 r 2 } and where the constants (r), (r), b(r), etc. are the discrete analogs of the constants in Section 2. With these estimates, the contrast functions C1 and C2 are easily computed. They are both of the form C(u, v) = C2 x (u, v) + C2 y (u, v), with Cx (u, v) = i,j j c ij Iu+i,v+j and Cy (u, v) = i,j i cij Iu+i,v+j , where the coef?cients cij are given by: (1) In the case of the contrast function C1, cij = 1 b(r) 1(i,j )?Br . (2) In the case of the contrast function C2 with two balls of radius r1 < r2, cij = 1 b(r1) 1(i,j )?Br 1 - 1 b(r2) 1(i,j )?Br 2 , where 1A denotes the characteristic function of the event A (its value is 1 if A is true and it is 0 otherwise). Therefore, the computations of the laws will be similar and they will be treated simultaneously using the coef?cients cij . When the contrast function C2 is used with two radii r1 < r2, we compute the Laplacian r I with the larger ball radius, that is with r = r2. Lemma 3. For both contrast functions C1 and C2, the vector (Cx (u, v), Cy (u, v), r I (u, v)) is a Gaussian random vector with mean 0 and covariance matrix of the form: ? ? 2 x 0 x , 0 2 y 0 x , 0 2 ? ? . In particular, we have that Cy is independent of (Cx , r I ).3286 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 Proof. The lemma is a consequence of the two following remarks. The ?rst one is that, in both cases for the contrast function, the coef?cients cij are symmetric: ci,j = c-i,j and ci,j = ci,-j . Thus, they satisfy that whenever k or l is odd then (i,j )?Br i k j l cij = 0. (3) The second remark is that the random vectors I i and I k are independent if i = k. And we thus have E[Ii,j Ik,l ] = 0 if i = k, (j, l) if i = k, where E[X] denotes the expectation (i.e., the mean value) of the random variable X. We can now compute the covariance matrix. For instance, let us start with E[CxCy ] = (i,j,k,l) j kc ij ckl E[Iu+i,v+j Iu+k,v+l ], = (i,j,l) j i c ij ci l(v + j , v + l ), = (j,l) j(v + j , v + l ) i i cij ci l = 0. Similar computations give E[Cyr I ] = 0 and 2 x := E[C 2 x ] = (i,j,l) j l c ij ci l(v + j , v + l ), 2 y := E[C 2 y ] = (i,j,l) i 2 cij ci l(v + j , v + l ), 2 := E[(r I ) 2 ] = 1 2 (r) (i,j,l)?r ((r) + i 2 + j 2 ) × ((r) + i 2 + l 2 )(v + j , v + l ), x , := E[Cxr I ] = 1 (r2) (i,j,l)?r 2 j c ij ((r) + i 2 + l 2 )(v + j , v + l ), where we have set r = {(i, j , l ) such that (i, j ) ? Br and (i, l) ? Br }. 4.4. Computation of the threshold As the ?rst- and second-order derivatives are no longer independent, we must compute the conditional law of the contrast function knowing that r I = 0. Proposition 2. Let C be one of the two contrast functions. Then, the random variable C 2 is distributed, conditionally on {r I =0}, as the sum of the square of two independent Gaussian random variables, with mean zero and respective variance 2 y and 2 x|=0 = 2 x 2 - 2 x , 2 , that is a Gamma law with parameters 1 2 and 1 2 ( 2 y + 2 x|=0 ). Remark. The threshold value s () de?ned by P(Cs ()|r I = 0) can no longer be computed explicitly but a numerical approximation is easy to get as the Gamma density is well known. Proof. Cy is independent of the pair (Cx , r I ). Thus, conditionally on {r I = 0}, the random variables Cy and Cx are still independent and the conditional law of Cy is the Gaussian distribution with mean 0 and variance 2 y . Now, if D2 := 2 x 2 - 2 x , = 0, then the law of the pair (Cx , r I ) has a density which is given by fx ,(t1, t2) = 1 2D e -1/2(t1,t2)(t1,t2) t where is the inverse of the covariance matrix, i.e., = 1 D2 2 -x , -x , 2 x . Let us recall that, if f denotes the Gaussian density of r I , then the law of Cx conditionally on r I =0 has a density given by fx ,(t1, 0) f(0) and so is Gaussian with mean zero and variance 2 x|=0 = D2 2 · This result is still valid when D = 0 since it implies that Cx and r I are proportional and thus the law of Cx conditionally on r I = 0 is Gaussian with mean 0 and variance 0 (it is not random anymore). 4.5. Experiments 4.5.1. Case of a piecewise constant object To begin with, we still study the piecewise constant object of Fig. 1 described in Section 3.4. Let us recall that this image represents a slice of the object that contains the symmetry axis. The three-dimensional object is obtained by rotation around the vertical axis that goes through the middle of the image. In that case, we will use the contrast function C1, which is simply the norm of the gradient. The experiments of Fig. 9 correspond to a ball radius r =12 pixels. The size of the images is 512 × 512 pixels.I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3287 Fig. 9. (a) Reconstructed object from a single noisy radiograph, (b) contrast value at the zero-crossings of the Laplacian for the contrast function C1 and (c) signi?cant edges. -300 -200 -100 0 100 200 300 0.02 0.06 0.10 0.14 0.18 0.22 0.26 0.30 0.34 0.38 0.42 Fig. 10. Numerical computation of the threshold as a function of the distance from the symmetry axis for the noise level = 4 in the case r = 12. We start with the image of the radiograph obtained after the application of matrix H to our initial image. Then a Gaussian white noise (with = 4) is added to this radiograph. Then tomographic inversion (application of the matrix H -1 ) is performed. This gives the image of Fig. 9(a). As we already mentioned it, the noise is not stationary, it is now correlated and its variance depends on the distance from the symmetry axis. For instance, if the standard deviation of the Gaussian white noise on the radiograph is = 4, the variance of the noise on the tomography is about 2 2 = 32 near the axis, 0.02 2 = 0.32 at a distance of 65 pixels from the axis and 8.10-3 2 = 0.128 at the edge of the image located on the right at 200 pixels from the axis. The threshold value can be numerically computed. As an illustration, we display in Fig. 10 the threshold value for the contrast C1 as a function of the distance from the symmetry axis (range from -256 to +256), in the case of a noise level = 4 and of a ball radius r = 12. Some edges are not declared signi?cant by our algorithm near the symmetry axis; the noise is too pronounced there to distinguish the true edges from the noise. In other words, the threshold value near the symmetry axis (see Fig. 10) can be higher than the contrast associated with small differences in the density of the object. In those images, an edge with a low contrast is signi?cant if it is far enough from the symmetry axis (see also Fig. 11). Even when the edges are signi?cant, a high level of noise leads to a decrease in the location accuracy and increased raggedness of the edges. Moreover, some details are lost because of the smoothing due to the size of the ball. Since the edges separate two materials, one included in another, they must be closed curves. Currently, a human operator has to close them manually. Our method yields open edges. This does not imply that there is no frontier between the materials but only that the noise level is too high to give an accurate position of the edge. Therefore, we can either close the curves manually, or with standard curve completion methods, but this will not tell which closure is better (i.e., the closest to the real shape). Several methods are usual for curve completion. Let us mention • Methods based on graph exploration [17]. • Methods inspired by neurobiology [18]. • Snakes with ?xed tips [19]. Choice of the ball radius: Let us compare the results obtained with different ball radii (see Fig. 11). When the ball radius is small, the edges are more accurate but some are not signi?cant: the smoothing of the noise is not enough to get rid of it. On the contrary, when the radius is large, most of the edges are detected but small details are lost because of this smoothing. The choice of the ball radius depends essentially on the complexity of the edges (or the degree of accuracy we want). If we know (from prior knowledge) that the edges are very smooth, a large radius gives very good results (see the outer edge of our studied object). On the contrary, if we look for very complicated objects, we must choose smaller radii but some edges (with low contrast) may become non-signi?cant. The choice of radius r = 12 corresponds to an appropriate balance between the size of typical details and an effective smoothing of the noise in our case. Another strategy is to use both values r = 6 and 12. The ?rst value gives accurate edges but some true edges are not detected, and the value r = 123288 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 Fig. 11. Signi?cant edges obtained with different ball radius: from left to right: r = 6, 12 and 20. gives more signi?cant edges but smoothed. A dif?cult question is then how to merge both pieces of information. Let us also mention that the “inverse” use of the formula of Proposition 1 in order to estimate a good value for the ball radius (as it was explained in the remark of Section 3.3) would lead to a value of the ball radius that depends on the distance from the symmetry axis (large values when close to the axis, small ones when far from it). Comparison with other methods: We will give here the results obtained with two other methods which have both the advantage of directly providing closed curves. • The ?rst method is the one introduced in Ref. [4]. One keeps only the meaningful level lines of the image, which are de- ?ned by: the minimum of the norm of the gradient along the level line is larger than a threshold T (). This threshold is computed from the gradient histogram of the image. The meaning of this de?nition is that such curves have a probability less than to appear in a pure noise image (with same gradient histogram as the original image). The results obtained with this method are shown in Fig. 12. In the ?rst row: we smooth the image of Fig. 9(a) by convolution with a Gaussian kernel with respective standard deviation 2 and 4 pixels. And then, in the second row, we display the respective meaningful level lines. This experiment clearly shows that, since the noise model is not adapted to the image (in particular, the non-stationarity is not taken into account), many false contours are detected. • The second method is the famous Mumford–Shah segmentation for piecewise constant images [20]. Given an observed image g0 de?ned on a domain D, one looks for the piecewise constant approximation g of g0 that minimizes the functional E (g) = D |g - g0| 2 + Length(K(g)), where Length(K(g)) is the one-dimensional measure of the discontinuity set of g (which is a set of curves denoted by K (g)) and is a parameter which weights the second term of the functional. The results obtained with this method are shown in Fig. 13 for three different values of , respectively; = 42 882, 12 220 and 11 621 (the size of the image is 512 × 512). The implementation used here is the regionmerging algorithm described in Ref. [21]. The main drawFig. 12. First row: the image of Fig. 9(a) is smoothed by convolution with a Gaussian kernel with respective standard deviation 2 and 4 pixels. Second row: the meaningful level lines of each image. back of this method is that there is no veri?cation that the obtained contours are not due to the noise, and moreover the user has to give the value of (or of the ?nal number of regions). We end this subsection with the results of Canny detector on the original image and on the images obtained after Gaussian ?ltering with respective standard deviation = 2 and 4. The results are shown in Fig. 14. On these edge maps, no thresholding has been performed. The next step would be to compute the signi?cant edge points. According to the de?nition of a signi?cant edge point, the question would be: I being a pure noise image (with known parameters), what is the law of |?I |(x) conditionally on x being a Canny edge point? Such a computation is dif?cult to perform and should be part of some future work. The obtained result would probably be close to the ones obtained with our method (see for instance Fig. 11).I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3289 Fig. 13. Results obtained with the Mumford–Shah segmentation for piecewise constant images, for three different values of (see text). From left to right, the number of regions in the segmented image is, respectively, 3, 6 and 7. Fig. 14. From left to right: results of the Canny edge detector for, respectively, the original image (Fig. 9(a)), the image smoothed by convolution with a Gaussian kernel of width = 2 (middle) and = 4 (see ?rst row of Fig. 12). Fig. 15. (a) Inhomogeneous object, (b) its noisy radiograph ( = 4) and (c) tomographic reconstruction. 4.5.2. Case of an inhomogeneous material Let us turn now to a more realistic case: the materials are not homogeneous and consequently the object is no longer piecewise constant (see Fig. 15(a)). The object is now composed of an inhomogeneous material (with constant gradient) with homogeneous objects inside. In the ?rst experiment, the noise level is =4. The results (Figs. 15 and 16) are similar to the ones obtained for a homogeneous object in the previous section. The contrast function used here is C1. With the contrast function C2, the edges are not signi?cant anymore in this case because the noise level is too large. In the second experiment (Fig. 17), the noise level is = 0.4. It is lower than the one used in the previous experiments, and this is why the artifacts of tomographic reconstruction are less pronounced here. The noise level here is such that the Fig. 16. (a) Contrast value at the zero-crossings of the Laplacian for the contrast function C1 (r =12) and (b) signi?cant edges of the inhomogeneous object.3290 I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 Fig. 17. (a) Inhomogeneous object, (b) its noisy radiograph ( = 0.4) and (c) tomographic reconstruction. Fig. 18. (a) Signi?cant edges with the C1 contrast function (r = 12): there are many false detections. (b) Signi?cant edges with the C2 contrast function (r1 = 6 and r2 = 12): only the “true” edges are obtained. threshold computed with the contrast function C1 is lower than the gradient of the material. In that case the use of the contrast function C1 fails as many detected edges remain signi?cant. This is illustrated by Fig. 18(a). Fig. 18(b) gives the signi?cant edges obtained with the contrast function C2 with two ball radii r1 = 6 and r2 = 12. With this contrast function, we eventually get only the “true” edges. Choice of the ball radii: Here again, we must discuss how the radii are chosen. The discussion concerning the choice of the radius of the outer ball is exactly the same as for the contrast function C1. That is why we kept the value of 12. The inner ball radius must not be too close to the outer one as the difference between the two gradients would be too small. But this radius must also be not too small in order to still perform a good smoothing of the noise. The choice of the half of the outer radius is usually a good compromise. 5. Discussion and conclusion In this paper, we have introduced a method for the detection of signi?cant edges in noisy images when the parameters of the noise are known. The general idea of the method is to make a statistical test to decide whether an edge point is due to the noise or not. The edge detector we have used here is the zero-crossings of an estimated Laplacian of the image. We did not use a Gaussian ?lter to smooth the image and estimate its Laplacian (as it is done in Marr–Hildreth detector [22]), but rather local polynomial approximations of the image. The reason for this is that such an estimator is more robust to noise. We think that future extensions of our work should mainly focus on: • performing the same kind of computations for Marr–Hildreth detector for Canny detector and more generally for “derivative-based” edge detector (including non-isotropic ones); • being able to merge the information obtained at different scales (which are the different ball radii in our case). Indeed, as we have shown it in the experimental section, it is sometimes necessary to use the gradient estimated at two different scales to obtain good results; • extending the framework to tomography without rotational symmetry, and more generally to other situations of nonstationary noise (including possibly non-Gaussian ones). References [1] I. Abraham, R. Abraham, J.-M. Lagrange, F. Lavallou, Méthodes inverses pour la reconstruction tomographique X monovue, Revue Chocs 31 (2005) (chocs@cea.fr). [2] J.M. Dinten, Tomographie à partir d’un nombre limité de projections: régularisation par des champs markoviens, Ph.D. Thesis, Université ParisSud, 1990. [3] J. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Mach. Intell. 8 (6) (1986) 679–698. [4] A. Desolneux, L. Moisan, J.-M. Morel, Edge detection by Helmholtz principle, J. Math. Imaging Vision 14 (2001) 271–284. [5] R. Touzi, A. Lopes, P. Bousquet, A statistical and geometrical edge detector for SAR images, IEEE Trans. Geosci. Remote Sensing 26 (1988) 764–773. [6] P. Qiu, S. Bhandarkar, An edge detection technique using local smoothing and statistical hypothesis testing, Pattern Recognition Lett. 17 (1996) 849–872. [7] D. Marimont, Y. Rubner, A probabilistic framework for edge detection and scale selection, in: Sixth International Conference on Computer Vision, 1998. [8] Y.G. Leclerc, S.W. Zucker, The local structure of image discontinuities in one dimension, IEEE Trans. Pattern Anal. Mach. Intell. 9 (3) (1987) 341–355. [9] F. Cao, P. Musé, F. Sur, Extracting meaningful curves from images, J. Math. Imaging Vision 22 (2005) 159–181. [10] M. Sezgin, B. Sankur, Survey over image thresholding techniques and quantitative performance evaluation, J. Electron. Imaging 13 (1) (2004) 146–165.I. Abraham et al. / Pattern Recognition 40 (2007) 3277 – 3291 3291 [11] O.D. Trier, T. Taxt, Evaluation of binarization methods for document images, IEEE Trans. Pattern Anal. Mach. Intell. 17 (3) (1995) 312–315. [12] A. Savitsky, M.J.E. Golay, Smoothing and differentiation of data by simpli?ed least squares procedures, Anal. Chem 36 (1964) 1627–1639. [13] J.B. Walsh, An introduction to stochastic partial differential equations, in: Ecole d’été de Probabilités de Saint-Flour XIV 1984, Lecture Notes in Mathematics, vol. 1180, Springer, Berlin, 1986. [14] T. Hida, Brownian Motion, Applications of Mathematics, vol. 11, Springer, Berlin, 1980. [15] T. Hida, H.H. Kuo, J. Potthoff, L. Streit, White Noise. An In?nite Dimensional Calculus, Mathematics and Its Applications, vol. 253, Kluwer Academic Publishers, Dordrecht, 1993. [16] R.N. Bracewell, The Fourier Transform and its Applications, second ed., Mc Graw-Hill, New York, 1978. [17] R. Deriche, J.P. Cocquerez, G. Almouzni, An ef?cient method to build early image description, 9th International Conference on Pattern Recognition, Roma, 1988. [18] P. Gaussier, J.P. Cocquerez, Utilisation des réseaux de neurones pour la reconnaissance de scènes complexes: simulation d’un système comprenant plusieurs aires corticales, Traitement du signal 6 (1991) 441 –466. [19] L.D. Cohen, Etudes des modèles de contours actifs et autres techniques de traitement d’image, Ph.D. Thesis, University Paris-Sud, 1990. [20] D. Mumford, J. Shah, Boundary detection by minimizing functionals, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, San Francisco, 1985. [21] G. Koep?er, C. Lopez, J.-M. Morel, A multiscale algorithm for image segmentation by variational method, SIAM J. Numer. Anal. 31 (1) (1994) 282–299. [22] D. Marr, E. Hildreth, Theory of edge detection, Proc. R. Soc. London B 207 (1980) 187–217. About the Author—ISABELLE ABRAHAM is currently a researcher at the CEA (Centre Energie Atomique) in Bruyeres le Chatel (near Paris). She is working on image processing and tomographic reconstruction. About the Author—ROMAIN ABRAHAM is a professor of mathematics at the University of Orleans. He is working on stochastic processes (like Brownian snakes), and also on some problems in image analysis. About the Author—AGNES DESOLNEUX is a researcher at the CNRS (Centre National de la Recherche Scienti?que) at University Paris 5. She is working on stochastic models for images, and the geometrical analysis of images. About the Author—SEBASTIEN LI-THIAO-TE is currently a PhD student at the Ecole Normale Superieure de Cachan. He is working with Bernard Chalmond and Benno Schwikowski (Institut Pasteur, Paris) on applied mathematics for proteomic.