UNIVERSITE PARIS IX DAUPHINE U.F.R. MATHEMATIQUES DE LA DECISION

THESE

pour l’obtention du titre de

DOCTEUR EN SCIENCES

spécialité : Mathématiques Appliquées arrêté du 23 novembre 1988

Sujet :

THEORIE DE LA DISCRIMINATION DE TEXTURES. AXIOMATISATION ET ALGORITHMES NUMERIQUES.

Candidat : Christian LOPEZ

Directeur de recherche : Jean-Michel MOREL

Jury :

Président :  Monsieur Yves MEYER
Rapporteur :  Monsieur Pietro PERONA
Rapporteur :  Monsieur Michel RASCLE
Examinateur :  Monsieur Gilles AUBERT
Examinateur :  Monsieur Michel BARLAUD
Examinateur :  Monsieur Pierre-Louis LIONS
Examinateur :  Monsieur Jean-Michel MOREL
Examinateur :  Monsieur Claude PUECH
Invité :  Monsieur Jean-Paul BOUJOT

Année de soutenance : 1993 L’université n’entend donner aucune approbation ni improbation aux opinions émises dans les thèses : ces opinions doivent être considérées comme propres à leurs auteurs. à mes parents,

à mes amis.

Yves Meyer a accepté de présider ce jury. Je le remercie de l’honneur qui m’est ainsi fait.

Pour avoir guidé mes recherches tout au long de la préparation de cette thèse et pour maintes discussions passionantes, je tiens à remercier Jean-Michel Morel. Grâce à lui, j’ai découvert la “Recherche”.

Je tiens à remercier Pietro Perona et Michel Rascle pour avoir accepté de rédiger le rapport de cette thèse.

Gilles Aubert, Michel Barlaud, Jean-Paul Boujot, Pierre-Louis Lions et Claude Puech participent à ce jury. Qu’ils soient assurés de ma sincère gratitude.

Pour une collaboration enrichissante, je remercie l’équipe de traitement d’image de Dauphine : Jean Pierre D’Alès, Francine Catté, Antonin Chambolle, Albert Cohen, Thierry Cohignac, Françoise Dibos, Frédéric Eve, Jacques Froment, Frédéric Guichard, Georges Koepfler, Bertrand Leroy, Agnès Tourin, sans oublier les membres extérieurs Luis Alvarez, Vicente Casalles et Tomeu Coll.

Pour m’avoir permis d’utiliser dans les meilleures conditions le support informatique indispensable à mes recherches, je tiens à remercier Michel Vanbreugel.

Je voudrais exprimer ma reconnaissance à Josette Levy et Magali Marc-Dibildos pour m’avoir facilité toutes les formalités administratives.

Je souhaite remercier tous les membres du CEREMADE pour leur acceuil chaleureux.

Enfin, je n’oublie pas l’UFR MD et notamment : Francine Catté, Annie Charles, Gisèle Fritz, Etienne Lanery, H. Pham, Suzanne Roverato, Patrick Saint-Pierre ainsi que Sylvie Goujard, Gérard Roé et les secrétaires.

Table des matières

Introduction. 1

1 Formalisation mathématique et équations fondamentales de la pyramide visuelle. 3 1.1 Bref aperçu historique. . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 La pyramide visuelle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2.1 Quel est l’input de la pyramide visuelle? . . . . . . . . . . . . . 7 1.2.2 Quel est l’ output de la pyramide visuelle? . . . . . . . . . . . . 8 1.2.3 Quels principes de base doit vérifier la pyramide visuelle ? … 8

1.3 Axiomatisation de l’analyse multiéchelle d’images et de formes et classification des modèles selon leur equation fondamentale. . . . . . . . . . 11 1.3.1 Analyse multiéchelle d’images. . . . . . . . . . . . . . . . . . . . 11 1.3.2 Axiomatisation de l’analyse multiéchelle de formes. . . . . . . . 13

1.4 Segmentation multiéchelle. . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.5 Une application : la discussion d’une théorie psychophysique, la Théorie

des Textons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

A. Le “théorème fondamental” de l’analyse d’images.

B. Morphologie Mathématique Affine Invariante.

2 Un algorithme multiéchelle de segmentation d’images par méthode variationnelle. 31 2.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 Principes généraux des procédés de segmentation. . . . . . . . . . . . . 32

2.2.1 Formalisation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.2 Choix des propriétés et algorithmes de segmentation. . . . . . . 34

2.3 Propriété de compacité des segmentations obtenues par croissance de régions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.1 Propriétés topologiques des segmentations affines 1-normales. . 38 2.3.2 Estimations sur des segmentations obtenues par croissance de

régions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.4 Un algorithme pyramidal de construction de segmentations affines 2normales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.4.1 Description de l’algorithme. . . . . . . . . . . . . . . . . . . . . 43 2.5 Application à la discrimination de textures. . . . . . . . . . . . . . . . 47 2.6 Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3 Axiomatisation de l’analyse de formes et application à l’hyperdiscrimination de textures. 57 3.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 Les quatre principes de l’analyse multiéchelle de formes. . . . . . . . . 58 3.3 L’équation fondamentale de l’analyse de formes. . . . . . . . . . . . . . 61 3.4 L’équation fondamentale de l’analyse d’images. . . . . . . . . . . . . . . 63 3.5 La Théorie des Textons de Julesz. . . . . . . . . . . . . . . . . . . . . . 64 3.6 Comment définir des densités de textons. Une approche axiomatique. . 66 3.7 Vérification de la théorie de Julesz. . . . . . . . . . . . . . . . . . . . . 67 3.8 Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

A. Preuve du lemme de renormalisation.

4 Axiomatisation de l’analyse de formes et application à l’hyperdiscrimination de textures. 83 4.1 Les axiomes de Julesz. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2 Formalisation des textons . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 Les équations fondamentales de l’analyse de formes et d’images. . . . . 86 4.4 Confirmation expérimentale de la Théorie des Textons. . . . . . . . . . 87

5 Analyse Numérique de l’équation fondamentale du traitement d’images. 91 5.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.2 Les quatre principes de l’analyse multiéchelle de formes et l’équation

fondamentale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.2.1 L’équation fondamentale de l’analyse de formes. . . . . . . . . . 97 5.2.2 L’équation fondamentale du traitement d’images. . . . . . . . . 97 5.2.3 Nécessité de propriétés d’invariance et de stabilité des schémas

numériques de l’équation fundamentale. . . . . . . . . . . . . . . 98

5.3 Schémas aux différences finies quasilinéaires. . . . . . . . . . . . . . . . 99 5.3.1 Schéma de Luis Alvarez. . . . . . . . . . . . . . . . . . . . . . . 100 5.3.2 Schéma de Luis Alvarez et de Frédéric Guichard. . . . . . . . . 100

5.4 Schémas morphologiques (d’ensemble). . . . . . . . . . . . . . . . . . . 102 5.4.1 Un schéma théorique. . . . . . . . . . . . . . . . . . . . . . . . . 102 5.4.2 Filtres médians itérés et mouvement dépendant de la courbure. 104

5.5 Schémas morphologiques purs (évolution de courbes). . . . . . . . . . . 105 5.5.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.5.2 Schémas quasilinéaires de diffusion intrinsèque. . . . . . . . . . 105 5.5.3 Un schéma géométrique d’évolution de courbes. . . . . . . . . . 107

5.6 Résultats expérimentaux. . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.6.1 Le schéma de diffusion anisotropique de Luis Alvarez. . . . . . . 109 5.6.2 Le schéma aux différences finies invariant d’Alvarez et Guichard. 110 5.6.3 ”Morphologisation” du schéma d’Alvarez et Guichard. . . . . . 111 5.6.4 Morphologie sur la grille. . . . . . . . . . . . . . . . . . . . . . . 111 5.6.5 Schémas purement morphologiques (évolution de courbes)…. 112

5.7 Conclusion et conjectures. . . . . . . . . . . . . . . . . . . . . . . . . . 113

A. Appendice. iv

Contents

Introduction. 1

1 Mathematical formalization and fundamental equations of the visual pyramid. 3 1.1 A short story of the subject. . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 The visual pyramid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2.1 What is the input of the visual pyramid ? . . . . . . . . . . . . 7 1.2.2 What is the output of the visual pyramid ? . . . . . . . . . . . . 8 1.2.3 What basic principles must obey the visual pyramid ? . . . . . . 8

1.3 Axiomatization of multiscale shape and image analysis and classification of the models depending on their fundamental equation. . . . . . . . . 11 1.3.1 Multiscale image analysis. . . . . . . . . . . . . . . . . . . . . . 11 1.3.2 Axiomatization of multiscale shape analysis. . . . . . . . . . . . 13

1.4 Multiscale segmentation. . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.5 An example : the discussion of a psychophysical theory, the Texton

Theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

A. The “fundamental theorem” of image analysis.

B. Affine Invariant Mathematical Morphology.

2 A multiscale algorithm for image segmentation by variational method. 31 2.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 General principles of segmentation devices. . . . . . . . . . . . . . . . . 32

2.2.1 Formalization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.2 Choice of properties and algorithms for segmentation. . . . . . . 34

2.3 Compactness properties of segmentations obtained by “merging”…. 36 2.3.1 Topological properties of affine 1-normal segmentations. . . . . . 38 2.3.2 Estimates on segmentations obtained by merging. . . . . . . . . 39

2.4 A pyramidal algorithm constructing 2-normal affine segmentations… 42

2.4.1 Description of the algorithm. . . . . . . . . . . . . . . . . . . . . 43 2.5 Application to texture discrimination. . . . . . . . . . . . . . . . . . . . 47 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3 Axiomatization of shape analysis and application to texture hyper-

discrimination. 57 3.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 The four principles of Multiscale Shape Analysis. . . . . . . . . . . . . 58 3.3 The fundamental equation of Shape Analysis. . . . . . . . . . . . . . . 61 3.4 The fundamental equation of Image Processing. . . . . . . . . . . . . . 63 3.5 The ”Texton Density” Doctrine of Julesz. . . . . . . . . . . . . . . . . . 64 3.6 How to define Texton Densities. An Axiomatic Approach. . . . . . . . 66 3.7 Experimental verification of Julesz Theory. . . . . . . . . . . . . . . . . 67 3.8 Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

A. Proof of the renormalization lemma.

4 Axiomatization of shape analysis and application to texture hyperdiscrimination. 83 4.1 Julesz axioms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2 Texton formalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 The fundamental equations of shape and image analysis. . . . . . . . . 86 4.4 Experimental confirmation of Julesz Theory. . . . . . . . . . . . . . . . 87

5 Numerical Analysis of the fundamental equation of image processing. 91 5.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.2 The four principles of Multiscale Shape Analysis and the fundamental

equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.2.1 The fundamental equation of Shape Analysis. . . . . . . . . . . 97 5.2.2 The fundamental equation of Image Processing. . . . . . . . . . 97 5.2.3 Invariance and stability requirements for numerical schemes of

the fundamental equation. . . . . . . . . . . . . . . . . . . . . . 98

5.3 Quasilinear finite difference schemes. . . . . . . . . . . . . . . . . . . . 99 5.3.1 Luis Alvarez’s scheme. . . . . . . . . . . . . . . . . . . . . . . . 100 5.3.2 Luis Alvarez, Frederic Guichard’s scheme. . . . . . . . . . . . . 100

5.4 Morphological (set evolution) schemes. . . . . . . . . . . . . . . . . . . 102 5.4.1 A theoretical scheme. . . . . . . . . . . . . . . . . . . . . . . . . 102 5.4.2 Iterated median filters and curvature motion. . . . . . . . . . . 104

5.5 Purely morphological (curve evolution) schemes. . . . . . . . . . . . . . 105 5.5.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.5.2 Quasilinear intrinsic diffusion schemes. . . . . . . . . . . . . . . 105 5.5.3 A geometrical curve evolution scheme. . . . . . . . . . . . . . . 107

5.6 Experimental results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.6.1 The Alvarez anisotropic diffusion scheme. . . . . . . . . . . . . 109 5.6.2 The Alvarez-Guichard rotationally invariant finite difference scheme.110 5.6.3 ”Morphologization” of the Alvarez-Guichard scheme. . . . . . . 111 5.6.4 Set morphology on a grid. . . . . . . . . . . . . . . . . . . . . . 111

5.6.5 Purely morphological (curve evolution) scheme. . . . . . . . . . 112 5.7 Conclusion and conjectures. . . . . . . . . . . . . . . . . . . . . . . . . 113

A. Appendix. viii

Introduction.

Le problème de la segmentation d’images a donné lieu à beaucoup de travaux dans les trente dernières années. Rappelons que segmenter une image c’est la “découper en régions homogènes séparées par des bords ou coupées par des failles”. La littérature propose d’innombrables algorithmes de segmentation. Cette abondance a été dénoncée et on a vu de plus en plus de formulations mathématiques qui permettaient de classer les anciens algorithmes et d’en créer d’autres plus performants (il suffit de penser aux méthodes variationnelles ou stochastiques). La segmentation de textures semblait particulièrement difficile à appréhender et, là aussi, de très nombreux algorithmes étaient publiés. Il semblait alors intéressant de voir ce qui les rapprochaient.

C’est dans ce cadre que s’inscrit le travail présenté ici. L’objectif est de formaliser mathématiquement la discrimination de textures et de faire une étude des algorithmes nécessaires pour l’expérimentation numérique.

Le premier chapitre constitue une introduction à la pyramide visuelle. Il peut être rapporté à une bibliographie dont les premiers articles datent des années 1960. En posant les principes généraux d’une analyse multiéchelle d’images, on aboutit, soit à une famille d’équations aux dérivées partielles géométriquement invariantes, soit à une formulation variationnelle. Nous donnons ensuite un programme en trois étapes permettant de discriminer les textures.

Le deuxième chapitre contient une étude approfondie du modèle variationnel, il s’agit d’étudier la fonctionnelle de Mumford et Shah avec approximation constante par morceaux. Nous introduisons pour cela une classe de segmentations, dites 2-normales que l’on peut relier aux segmentations déterminées par des algorithmes de croissance de régions introduits par Pavlidis. Ensuite, en suivant le programme en trois étapes du chapitre précédent mais en choisissant pour analyse multiéchelle l’analyse en ondelettes, nous parvenons à discriminer des textures naturelles telles que des peaux de reptiles, des tissages, des matériaux... L’analyse en ondelettes permet dans ce cas d’extraire deux caractéristiques importantes: l’échelle et l’orientation.

Dans les chapitres 3 et 4, nous traduisons l’axiomatique informelle de Julesz. Nous montrons que l’analyse des textons revient à faire une analyse de formes qui fait intervenir la courbure. Pour éviter de traiter de nombreuses formes en parallèle, nous en déduisons une analyse multiéchelle d’images géométriquement invariante qui fait apparaître la courbure des isophotes de l’image. En suivant le programme en trois étapes du premier chapitre, il est alors possible de vérifier la théorie de Julesz sur ses propres exemples de textures et en même temps notre formalisation.

Le cinquième chapitre décrit l’analyse numérique des équations d’évolution obtenues dans les chapitres 4 et 5. Nous insistons sur le fait que les algorithmes doivent vérifier au mieux les propriétés d’invariance qui ont permis d’aboutir aux équations. Nous notons aussi le lien entre le filtre médian et le schéma numérique de Merriman-Bence-Osher.

Chapitre 1

Formalisation mathématique et équations fondamentales de la pyramide visuelle.

Ce chapitre a été écrit en collaboration avec Jean-Michel Morel et va paraître dans Matapli.

Résumé. Une image numérique se modélise comme une fonction réelle u_{0}(x) définie sur le plan (le niveau de gris en x). Le pivot de la théorie de la vision et du traitement d’images est le concept d’analyse multiéchelle (ou ”scale space”). L’analyse multiéchelle associe à u(0) = u_{0} une séquence d’images simplifiées u(t, x) qui dépend d’un paramètre, t > 0 abstrait : l’échelle. L’image u(t, x) est l’analyse de u_{0} à l’échelle t (cf. fig. 1,2). La formalisation du scale space s’est accélérée dans les dix dernières années ; une dizaine de théories de l’analyse multiéchelle des images, des ”formes” et des ”textures” ont été proposées et des travaux récents ont permis de formaliser le champ tout entier. Ainsi, il suffit de quelques principes pour caractériser et regrouper plus de dix théories en analyse d’images, de formes et de textures. Ces principes sont

• la causalité (un concept de théorie de la vision qui se ramène au principe du maximum et à l’hypothèse d’une architecture pyramidale),

• l’invariance euclidienne et/ou affine, qui traduit le fait que le traitement ne dépend pas de la distance et de l’orientation de l’image analysée (fig. 2),

• l’invariance morphologique qui signifie que le traitement de l’image ne dépend pas d’un changement de l’échelle de contraste.

La caractérisation et la classification des nombreuses théories de l’analyse multiéchelle d’images et de formes sont obtenues en identifiant leurs équations aux dérivées partielles sousjacentes. La caractérisation axiomatique permet du même coup de démontrer que certaines d’entre elles sont identiques. Une autre conséquence de la réduction à des E.D.P. est qu’on peut améliorer la plupart des algorithmes et en proposer d’autres plus invariants. Parmi les théories multiéchelle qui sont maintenant entièrement axiomatisées et unifiées, nous pouvons citer

-le raw primal sketch de Hildreth et Marr,

-l’ entropy scale space de Kimia, Tannenbaum et Zucker,

-la théorie des textons de Julesz,

-le dynamic shape de Koenderink et Van Doorn,

-le curvature primal sketch de Mackworth et Mocktarian, Asada et Brady, -la morphologie mathématique de Matheron, Serra et l’école de Fontainebleau, -la diffusion anisotropique de Perona et Malik, -l’affine scale space of curves de Sapiro et Tannenbaum, -le scale space de Witkin, Koenderink,...

La classification de ces théories multiéchelle ne conduit pas à les conserver toutes, mais plutôt à proposer l’unique théorie jouissant simultanément de toutes les bonnes propriétés de stabilité et d’invariance qu’elles présentent.

Cette analyse multiéchelle la plus invariante (à la fois affine et morphologique) peut être définie par une simple équation aux dérivées partielles (qui est une équation de la chaleur ”intrinsèque”),

\frac{∂u}{∂t} = |Du|(t.div(Du/|Du|))^{1/3}, u(0, x) = u_{0}(x)

où u(t, x) désigne l’image analysée à l’échelle t et au point x. (L’équation admet une unique solution de viscosité au sens de Crandall-Ishii-Lions [C.I.L.]) Cette équation (Affine Morphological Scale Space, A.M.S.S.) traite séparément tous les ensembles de niveau d’une image et réalise donc une opération morphologique. De plus, chaque ensemble de niveau de l’image est traité par l’équation comme une forme. Par conséquent, on a une théorie commune et unique pour l’analyse de formes et l’analyse d’images.

Une formalisation multiéchelle de la segmentation d’images peut être développée sur des principes analogues et conduit à des algorithmes de segmentation multiéchelle d’images. Une segmentation multiéchelle associe à une image initiale u_{0} une séquence u(t, x), K(t, x), où u(t, x) est l’image simplifiée à l’échelle t et K(t, x) les bords des régions de l’image à l’échelle t. Des dizaines d’algorithmes de segmentation ont pu être regroupés en un seul grâce à la formalisation.

En guise d’application de la méthode axiomatique, nous montrons comment les formalisations multiéchelle précédentes permettent de discuter la théorie psychophysique de discrimination des textures (fig. 4,5). Le résultat est inattendu ...

PLAN

1. Bref aperçu historique.

2. La pyramide visuelle.

2.1 Quel est l’input de la pyramide visuelle? 2.2 Quel est l’output de la pyramide visuelle? 2.3 Quels principes de base doit vérifier la pyramide visuelle?

3. Axiomatisation de l’analyse multiéchelle d’images et de formes et classification des modèles selon leur équation fondamentale.

3.1 Analyse multiéchelle d’images. 3.2 Axiomatisation de l’analyse multiéchelle de formes.

4. Segmentation multiéchelle.

5. Une application : La discussion d’une théorie psychophysique, la Théorie des Textons.

Appendices A. Le ”théorème fondamental” de l’analyse d’images. B. Morphologie Mathématique Affine Invariante. Planches 1,2,3,4,5.

1.1 Bref aperçu historique.

La théorie de la Vision et le traitement d’images (Computer Vision) essaient de répondre à la question suivante :

Comment une information locale due à la lumière arrivant sur la rétine (ou sur un capteur) peut-elle être transformée en une perception globale des objets ?

A partir des années soixante, la communauté étudiant la vision a beaucoup utilisé les ordinateurs pour exploiter de nouvelles possibilités d’expérimentation sur des images numériques. Les nouvelles technologies ont permis de réaliser des mesures précises sur les performances du système visuel humain et ont ouvert la voie aux premières expériences de Vision par ordinateur. Les développements simultanés de la psychophysique et du traitement d’images ont convergé vers une nouvelle doctrine : l’existence d’une vision bas niveau. L’histoire de cette doctrine est bien expliquée dans le livre de David Marr ”Vision” ([Marr]) et nous allons juste donner quelques faits.

D’une part, plusieurs expériences psychophysiques dues à Bela Julesz et son école, ont prouvé que la reconstruction de l’environnement spatial à partir d’une information binoculaire était une action automatique, un réflexe ne dépendant d’aucun apprentissage. Julesz a aussi étudié la perception préattentive des textures et a prouvé l’existence d’un processus permettant de discriminer des textures sans connaissance a priori ni apprentissage (fig. 4,5). Le processus de discrimination est rapide et parallèle. Julesz et son école l’ont décrit d’abord par un formalisme statistique, puis géométrique.

Ces expériences, ainsi que les expériences neurobiologiques de Hubel et Wiesel, ont donné une preuve d’existence, dans les premières millisecondes du processus de perception, d’une série d’opérations parallèles, rapides et irréversibles appliquées à l’information rétinienne et donnant déjà une information très riche et très utile pour la compréhension de l’image.

1.2 La pyramide visuelle.

C’est cette série d’opérations que nous appelons ”pyramide visuelle”. On doit la penser comme un algorithme, mais pas au sens de Turing; dans un sens plus général où un algorithme est une boîte noire transformant l’input en un output d’une manière déterministe par un processus physique : cette machine est supposée fonctionner physiquement dans le cerveau. Il y a deux questions à distinguer : comment cette machine fonctionne-t-elle dans le cas du cerveau et quel est le traitement abstrait appliqué à l’information ? (Pour reprendre des termes dus à Chomsky on distingue ”performance et compétence”). En effet, la seconde question est simplement une question mathématique, tandis que la première est très liée à la neurobiologie.

Nous allons nous intéresser au problème mathématique et le traiter d’une manière assez générale en répondant aux trois questions suivantes

a) Quel est l’ input de la pyramide visuelle ?

b) Quel est son output ?

c) Quels principes doit vérifier le traitement effectué par la pyramide visuelle?

1.2.1 Quel est l’input de la pyramide visuelle ?

Le modèle le plus simple en traitement d’images consiste à définir une image comme une fonction de luminance u_{0}(x) en chaque point x d’un domaine du plan. Ce domaine, qui peut être le plan lui-même, est un modèle de la rétine ou d’une surface photosensible. Dans ce qui suit, nous allons prendre le plan pour simplifier. Bien sûr, la donnée de u_{0}(x) n’est pas absolue en théorie de la perception. Elle doit être considérée comme un élément d’une classe d’équivalence. Si y est un vecteur du plan, u_{0}(x − y), qui représente l’image translatée de y, est aussi un élément de la classe d’équivalence. De la même manière, u_{0}(Rx) où R est une isométrie du plan est encore un élément de la même classe d’équivalence car l’analyse d’image ne doit pas dépendre d’un repère orthonormé particulier (l’analyse doit être géométriquement intrinsèque). On peut imposer plus généralement que u_{0}(Ax), où A est une transformation affine du plan, appartienne aussi à la même classe d’équivalence. En effet, une image plane peut être perçue de la même manière par un observateur quelle que soit sa position par rapport à l’image (penser à un tableau dans une exposition). Par conséquent, l’input u_{0}(x) est supposé équivalent à n’importe quelle de ses anamorphoses u_{0}(Ax) où A est une transformation affine inversible du plan (effet de perspective pour un objet éloigné de l’oeil) (voir fig.2).

Enfin, u_{0}(x) ne donne généralement pas d’information sur le nombre de photons émis et reçus sur le capteur. Par conséquent, la classe d’équivalence sera g(u_{0}(x)), où g représente n’importe quelle fonction de contraste dépendant du capteur. Cette dernière hypothèse, qui consiste à dire que seuls les isophotes comptent, est un des principes de la morphologie mathématique [Mat,Ser]. C’est pourquoi nous l’appellerons hypothèse morphologique. Il faut penser à deux vues d’une même scène prises avec deux appareils différents (ou avec le même appareil à des instants différents) : une des images se déduit de l’autre par un changement de contraste g.

En résumé, une image est une classe d’équivalence de fonctions u_{0}(Ax) où A est n’importe quelle isométrie dans le modèle euclidien, A est une transformation affine inversible dans le modèle affine. Une image peut se définir aussi comme une classe de fonctions g(u_{0}(x)) où g est n’importe quelle fonction strictement monotone et continue dans le modèle morphologique. On va associer ces modèles et considérer un modèle affine morphologique c’est à dire une classe d’équivalence de la forme (g(u_{0}(Ax))).

1.2.2 Quel est l’output de la pyramide visuelle ?

Partant de l’information locale de luminance, chaque couche de la pyramide visuelle est supposée fournir une information bas niveau de plus en plus globale de l’image. Quelle que soit son utilisation, la plupart des modèles définissent l’output comme

- une image régulière (à partir de laquelle des caractéristiques fiables (features) peuvent être extraites par des opérateurs différentiels) (voir fig.1)

- ou bien une segmentation, c’est à dire, soit une décomposition de l’image en régions homogènes (”segmentation forte”), avec des bords ou un ensemble de points frontière ou edge map (voir fig.3).

Dans les deux cas, l’output dépend de deux variables : une variable x qui représente le centre d’un voisinage spatial et une variable t qui peut être identifiée ou corrélée avec

- la ”hauteur” de la pyramide visuelle (ou la distance par rapport à la première couche, i.e. la rétine). Cette distance correspond au temps biologique entre l’arrivée sur la rétine et la première arrivée sur une couche donnée.

-Le degré de globalité de l’information locale dans la couche considérée, c’est à dire, la taille du voisinage dans la rétine u_{0} qui influe sur la valeur de u(t, x).

En résumé, l’output de la pyramide visuelle est

- soit une image multiéchelle u(t, x)

- soit une carte de bords multiéchelle K(t, x), où t est un paramètre qui peut être identifié à un temps d’analyse ou à une mesure de globalité spatiale de l’information de u(t, x). La base de la pyramide est l’image originale u(0, x) = u_{0}(x).

1.2.3 Quels principes de base doit vérifier la pyramide visuelle?

La causalité.

Ces principes proviennent du fait que l’on modélise la vision préattentive : aucun feedback n’est acceptable ; l’information visuelle est traitée en parallèle à travers une séquence de filtres. Ainsi, ce qui se passe aux grandes échelles n’agit pas sur les petites échelles : la pyramide provoque un flux des petites échelles vers les grandes échelles (fine to coarse). De plus, nous supposons que l’output à l’échelle t peut être calculé à partir de l’output à l’échelle t − h avec h > 0. (Prendre en compte que la pyramide visuelle est une succession de filtres à travers lesquels de l’information nouvelle est constamment traitée.)

Pour formaliser cette relation, nous appelons T_{t} la transformation qui associe à l’image originale u_{0} son image filtrée à l’échelle t, T_{t}u_{0}. (De la même manière, nous notons aussi T_{t} la transformation qui associe à un ensemble K de bords un autre ensemble T_{t}K de bords à l’échelle t). Cette transformation est obtenue en composant des filtres de transition que nous notons T_{t+h,t} et nous avons obtenu le premier principe régissant la vision préattentive

Structure pyramidale T_{t+h} = T_{t+h,t}T_{t}.

De plus, on supppose que l’opérateur T_{t+h,t} agisse toujours localement. En d’autres termes, (T_{t+h,t}u_{0})(x) (resp. T_{t+h,t}K(x)) doit essentiellement dépendre des valeurs de u_{0}(y) (resp. K(y)) quand y se trouve près de x.

Nous allons donner deux versions formelles de cette hypothèse de localité. Donnons tout de suite son interprétation ”physique” : si les éléments de base de la pyramide sont supposés être des ”neurones”, cela signifie qu’un neurone est d’abord influencé par ses voisins. Un argument permet de justifier cette remarque : seuls les neurones proches peuvent s’influencer sans délai. Terminons par une exigence intuitive qui est appelée en traitement d’image ”causalité”. Puisque la pyramide visuelle est supposée donner de plus en plus d’information globale sur l’image elle-même et ses caractéristiques, il est clair que lorsque l’échelle croit, il ne devrait pas y avoir création de nouvelles structures durant l’analyse multiéchelle : l’image et les bords à l’échelle t′ > t doivent être plus simples que ceux à l’échelle t. L’hypothèse de causalité doit, bien sûr, être formalisée. Sa formalisation a été étudiée par Hummel [Hum], Koenderink [Koen1,2], Yuille [Yui], Witkin [Wit], Perona et Malik [PeMa] dans le cadre du traitement d’image, par Kimia, Tannenbaum et Zucker [KTZ] dans le cadre de l’analyse de forme et par Muerle et Allen [MuAl], Brice et Fennema [BrFe], Horowitz [Hor], Pavlidis [Pav] dans des travaux sur la segmentation d’image.

Le résultat de la discussion dans le cas du traitement d’image est que la causalité doit être formalisée comme un principe de comparaison : Si une image u est localement plus claire qu’une autre image v, alors cet ordre doit être conservé au cours de l’analyse (principe de comparaison). En termes formels,

Principe de comparaison locale

si u(y) > v(y) pour y dans un voisinage du point x et y ≠ x, alors on doit avoir pour h assez petit

(T_{t+h,t}u)(x) ≥ (T_{t+h,t}v)(x)

Dans le cas de la détection de contour, il y a plusieurs formalisations, mais la plus simple établit qu’il n’y a pas de bord créé lorsque l’échelle croît, c’est à dire que T_{t′}K est inclus dans T_{t}K si t′ > t.

L’invariance morphologique et affine.

Outre l’exigence de causalité, nous devons garder à l’esprit que la pyramide agit sur des classes d’équivalence de fonctions de la forme g(u_{0}(Ax)), où g est une fonction monotone continue et A est soit une isométrie soit une transformation affine inversible du plan. Par conséquent, l’ output ne devrait pas dépendre de u_{0} mais de sa classe d’équivalence. Ainsi, les opérateurs de transition T_{t+h,t}u doivent commuter avec les ”perturbations” g et A. Dans le cas du changement de contraste g, on peut facilement le traduire sous la forme

Axiome morphologique gT_{t+h,t} = T_{t+h,t}g,

ce qui signifie que le changement de contraste et l’analyse multiéchelle peuvent être appliqués dans n’importe quel ordre. Si A est une isométrie, la même forme de relation doit être vraie. Notons Au la fonction définie par Au(x) = u(Ax). Alors nous posons

Invariance euclidienne AT_{t+h,t} = T_{t+h,t}A

Examinons maintenant le cas d’une transformation affine A quelconque. La relation de commutation ne peut plus être aussi simple car A peut réduire ou agrandir l’image (penser au cas où A est un zoom défini par Au(x) = u(λx) avec λ une constante positive). Puisque le zoom a changé l’échelle de l’image, nous pouvons uniquement imposer une propriété de commutation faible

Invariance affine AT_{t′+h′,t′} = T_{t+h,t}A

où t′ et h′ dépendent de t, h et A. Il a été prouvé dans [AGLM] (et c’est un premier résultat mathématique) que l’on peut toujours imposer t′ = |det(A)|^{1/2}t et h′ = |det(A)|^{1/2}h (modulo un changement d’échelle t → ϕ(t)).

En résumé, la pyramide visuelle T_{t} doit vérifier

-Causalité : T_{t+h} = T_{t+h,t}T_{t} , T_{t,t} = Id, et (T_{t+h,t}u)(x) > (T_{t+h,t}v)(x) si u(y) > v(y) pour y dans un voisinage de x et y ≠ x.

Dans le cas de l’analyse multiéchelle de bords, cette dernière hypothèse est remplacée par T_{t}K ⊂ T_{t}K′ si K ⊂ K′.

-Invariance morphologique : gT_{t+h,t} = T_{t+h,t}g quel que soit le changement de contraste g.

-Invariance euclidienne : AT_{t+h,t} = T_{t+h,t}A quelle que soit l’isométrie A du plan.

-Invariance affine : AT_{t′+h′,t′} = T_{t+h,t}A avec t′ = |det(A)|^{1/2}t et h′ = |det(A)|^{1/2}h. Il faut noter que l’invariance affine implique l’invariance euclidienne.

-Linéarité : T_{t+h,t}(au + bv) = aT_{t+h,t}(u) + bT_{t+h,t}(v) Nous ajoutons cette propriété optionnelle car elle a beaucoup été utilisée en ”Computer Vision” [Wit,Koe].

Il y a donc cinq axiomes principaux et nous allons voir qu’ils permettent de classer et caractériser complétement les théories du traitement multiéchelle d’images et de formes et qu’ils unifient et améliorent certaines d’entre elles.

1.3 Axiomatisation de l’analyse multiéchelle d’images et de formes et classification des modèles selon leur équation fondamentale.

1.3.1 Analyse multiéchelle d’images.

D’abord, il est bien reconnu en traitement d’images (Koenderink [Koe2], Hummel [Hum], Yuille [Yui]) que si une analyse multiéchelle est causale, si elle vérifie l’invariance euclidienne et si elle est linéaire, elle est régie (modulo un changement d’échelle t → ϕ(t)) par l’équation de la chaleur

\frac{∂u}{∂t} = Δu

(C’est à dire que la fonction u(t, x) = (T_{t}u)(x) vérifie l’equation de la chaleur). Ce modèle de l’analyse multiéchelle est dû à Hildreth et Marr (le ”raw primal sketch” [Marr]) et Witkin et Koenderink (le scale space [Wit], [Koe]).

Que se passe-t-il si on élimine l’axiome de linéarité qui fixait un seul modèle ?

Comme l’ont remarqué Perona et Malik [PeMa] dans leur théorie non linéaire du scale space, (Anisotropic diffusion [PeMa]) on peut obtenir des équations non linéaires

\frac{∂u}{∂t} = F (D^{2}u, Du, u, x, t)

L’étude complète des modèles non linéaires a été menée dans [AGLM] : si l’analyse multiéchelle d’image T_{t} est causale alors elle est une solution de viscosité de l’équation précédente et la fonction F est croissante par rapport à son premier argument (pour une démonstration rapide, voir Appendice A). Le cas particulier de l’équation de la chaleur correspond au cas F (A, p, u, x, t) = trace(A). Par conséquent, tous les modèles multiéchelle peuvent être classés et de nouveaux modèles plus invariants peuvent être proposés. Par exemple, si l’analyse est causale et vérifie l’invariance euclidienne et morphologique, alors elle est régie par une équation de la forme

\frac{∂u}{∂t} = |Du|F (div(Du/|Du|), t)

où div(Du/|Du|)(x) peut être interprété comme la courbure de la ligne de niveau de l’image u(t, x) passant par x et F (s, t) est une fonction croissante par rapport à la variable s.

Un cas particulier important apparaît lorsque F est une fonction constante : si F = 1 (ou F = −1) alors l’équation devient \frac{∂u}{∂t} = |Du| (resp. \frac{∂u}{∂t} = −|Du|), qui correspond à la dilatation morphologique si le signe est ”+” et à l’érosion morphologique si le signe est ”-”. La dilatation et l’érosion sont les opérations de base d’une des premières théories du traitement d’images, dues à Matheron et son école, la ”Morphologie Mathématique” ([Mat], [Ser]). La dilatation à l’échelle t est définie par

D_{t}u_{0}(x) = sup_{y∈B(x,t)} u_{0}(y),

où B(x, t) est un ensemble centré en x : l’élément structurant (en général une boule de rayon t). Par exemple, si u_{0} est la fonction caractéristique d’un ensemble X, D_{t}u_{0} est la fonction caractéristique du t-voisinage de X. Pour l’érosion, on remplace le ”sup” par un ”inf”. Comme nous l’a fait remarquer Michel Rascle, un autre exemple important vérifiant l’axiomatique multiéchelle morphologique et donc une E.D.P. parabolique (dégénérée) est la famille des ”zooms” de rapport t : u(t, x) = u_{0}(tx). Le théorème précédent s’applique et l’équation est

\frac{∂u}{∂t} = \frac{1}{t}(∇u.x)

Noter que le zoom d’une image discrétisée sur une grille n’est pas une opération numériquement triviale... Les exemples précédents ne sont que des cas particuliers. En fait, il existe beaucoup d’autres possibilités d’analyse morphologique ! Si nous posons F (s, t) = s par exemple, nous obtenons la mean curvature motion (mouvement par courbure moyenne)

\frac{∂u}{∂t} = t.|Du|div(Du/|Du|)

Cette équation provient d’une reformulation faite par Osher et Sethian [OsSe] d’un modèle de géométrie différentielle étudié par Grayson, Gage-Hamilton [Gray, GaHa]. De plus, cette équation est très proche du modèle de ”diffusion anisotropique” de Perona et Malik et d’une équation permettant restauration d’images due à Osher et Rudin [OsRu] :

\frac{∂u}{∂t} = div(Du/|Du|)

Il y a un modèle encore plus invariant et nouveau : il est prouvé dans [AGLM] qu’il existe une unique analyse multiéchelle d’image causale et vérifiant à la fois l’invariance affine et l’invariance morphologique. Son équation est

\frac{∂u}{∂t} = |Du|(t.div(Du/|Du|))^{1/3}

Nous appellerons ce modèle A.M.S.S. (Affine Morphological Scale Space). Nous allons mieux comprendre cette équation dans le cadre de l’analyse multiéchelle de formes. Sapiro et Tannenbaum [SaTa1,2], qui ont découvert indépendamment cette équation, lui ont donné une interprétation géométrique remarquable.

1.3.2 Axiomatisation de l’analyse multiéchelle de formes.

Cette analyse multiéchelle prend pour input une ”silhouette”, c’est à dire un ensemble fermé du plan dont le bord est une courbe de Jordan fermée x_{0}(s) et essaie de faire évoluer la courbe. Une adapation facile de la formalisation précédente peut être faite pour axiomatiser l’analyse de forme. Lopez et Morel [LoMo] ont montré que quatre principes de base engendrent un unique modèle pour l’analyse multiéchelle de formes. Ces principes sont -Principe pyramidal : T_{t+h} = T_{t+h,t}T_{t}.

-Principe d’inclusion locale des formes : Si une forme est incluse localement dans une autre, cela reste vrai lorsqu’on applique T_{t+h,t} pour h petit.

-Invariance affine : L’analyse de forme commute avec les transformations affines du plan.

-Régularité : L’évolution du rayon d’un disque est régulière par rapport à t.

Le seul modèle d’analyse de forme obtenu à partir de ces principes est donné par une équation d’évolution portant sur le bord de la forme. Notons curv(x) la courbure en un point x du bord et n⃗(x) le vecteur unitaire normal en ce point (on suppose le bord de classe C^{2} au voisinage de x), alors l’équation d’évolution est

ẋ = t^{1/3}curv(x)^{1/3}n⃗(x), x(0, s) = x_{0}(s)

Cette équation est équivalente à l’équation fondamentale de l’analyse d’images (A.M.S.S.)

\frac{∂u}{∂t} = |Du|(t.div(Du/|Du|)^{1/3}

puisqu’il est prouvé dans [AGLM] que les isophotes d’une image analysée avec la seconde équation évoluent suivant la première. Finalement, les deux modèles sont équivalents à l’affine invariant scale space (ASS) de courbes découvert par Sapiro et Tannenbaum. Le modèle de Sapiro et Tannenbaum permet une interprétation géométrique simple de l’équation fondamentale. Considérons deux manières de paramétrer une courbe de Jordan régulière x(s) :

• soit de manière euclidienne-invariante en imposant la relation |x_{s}| = 1.

• soit de manière affine-invariante en imposant la relation |det(x_{s}, x_{ss})| = 1.

Paramétrons la courbe analysée par son ”abscisse affine” s. L’abscisse affine s est une paramétrisation affine invariante d’une courbe ( car |det(Ax_{s}, Ax_{ss})| = | det A||det(x_{s}, x_{ss})| quelle que soit la transformation linéaire A ). Alors l’affine invariant scale space est équivalent à l’équation de la chaleur intrinsèque (affine)

\frac{∂x(s, t)}{∂t} = t^{1/3} \frac{∂^{2}x(s, t)}{∂s^{2}}, x(0, s) = x_{0}(s)

Cette dernière équation est donc aussi équivalente aux précédentes! Pour finir, donnons une liste de théories de l’analyse de forme.

Les théories qui sont axiomatisées par la présentation précédente sont le dynamic shape de Koenderink-Van Doorn [KoDo], le curvature primal sketch de Mackworth-Mocktarian [MaMo], Asada-Brady [AsBr], l’ entropy scale space de Kimia, Tannenbaum et Zucker [KTZ] (ces théories sont de fait régies par la Mean Curvature Motion et sont donc équivalentes entre elles) et finalement l’affine scale space de Sapiro-Tannenbaum [SaTa1,2]. Il y a des théorèmes d’existence, d’unicité et de régularité des solutions pour l’évolution euclidienne et affine de courbes; ils sont dus à Gage, Hamilton, Grayson, Sapiro et Tannenbaum ([GaHa], [Gra], [SaTa]...).

1.4 Segmentation multiéchelle.

Etant données une image u_{0} et une carte initiale de bords K_{0}, où K_{0} est un ensemble de longueur de Hausdorff finie du plan, on cherche à définir une séquence (K(t)) de segmentations multiéchelles, et une séquence d’images filtrées u(t) dont K(t) sera l’ensemble des bords. Nous allons nous contenter de décrire une méthode de segmentation multiéchelle pour laquelle l’axiomatisation est presque complète. Les principes que l’on impose sont les suivants.

On note T_{t}K_{0} = K(t).

-Causalité K(t′) ⊂ K(t) si t < t′ et T_{t,s}T_{s} = T_{t}

-Invariance euclidienne T_{t+h,t}R = RT_{t+h,t} pour toute isométrie R.

-Invariance par changement d’échelle Si Z_{λ} désigne le ”zoom” de rapport λ alors

T_{t}Z_{λ} = Z_{λ}T_{λ^{1/2}t}

Il nous faut aussi un axiome fixant les relations entre les bords K(t) à l’échelle t et l’image u(t) à l’échelle t. On choisit ici le principe le plus simple que Mumford et Shah [MuSh] appellent principe de bande dessinée (cartoon).

-Bande dessinée : u(t) est localement constante et égale à la moyenne de u_{0} sur IR^{2} \ K(t).

On remarquera que le principe de causalité et le principe de bande dessinée fixent le type d’algorithme utilisé : il s’agit d’une croissance de régions. On appelle région de l’image à l’échelle t toute composante connexe de IR^{2} \ K(t). Comme pour t′ > t, K(t′) ⊂ K(t), on voit que les régions à l’échelle t′ sont des unions de régions à l’échelle t (c’est une autre manière d’exprimer la causalité). Donc, pour décrire complétement l’analyse multiéchelle, il suffit de donner un critère de fusion des régions. Le critère le plus simple de tous, et qui est effectivement invariant par similitude (isométrie et homothétie) est le critère de Mumford et Shah : la fusion doit faire décroître l’énergie suivante

E(u(t), K(t)) = ∫_{IR^{2}\K(t)}(u(t, x) − u_{0}(x))^{2}dx + t.longueur(K(t))

Bien que l’usage de cette énergie ne soit pas totalement axiomatisé, il faut remarquer qu’elle représente un point d’aboutissement d’une centaine de théories et d’algorithmes de traitement d’images. Citons notamment la théorie stochastique de la segmentation de Geman-Geman [Gem], les modèles variationnels de Mumford et Shah [MuSh], Terzopoulos [Terz], et les algorithmes de croissance de régions de Brice-Fenema [BrFe], Muerle-Allen [MuAl], Pavlidis [Pav] (une classification de ces algorithmes est présentée dans un livre à paraître [MoSo]). L’algorithme qui se déduit des considérations précédentes consiste à fusionner récursivement les régions de la segmentation (K(t), u(t)) de manière à faire diminuer le plus possible le critère (avec t′ à la place de t). La théorie des segmentations minimisant l’énergie a été faite par Mumford et Shah [MuSh], Morel-Solimini [MoSo], Koepfler-Lopez-Morel [KLM] et elle a l’avantage de donner des résultats d’existence et de régularité que l’image u_{0} soit réelle ou vectorielle (Segmentation multicanaux). Une version affine invariante du même algorithme a été définie et testée dans [GoBa].

Une dernière remarque concernant l’invariance morphologique. La segmentation recherche des régions les plus homogènes possibles et des bords les plus contrastés possibles. Cela contredit l’invariance morphologique car celle-ci exclut tout critère fondé sur le contraste. Toutefois, l’image initiale n’a pas à être nécessairement dépendante du contraste. Par exemple, la segmentation peut porter sur une image multicanaux (vectorielle) dont tout canal résulte d’un filtrage calculant des invariants morphologiques. Il en résulte (nous allons voir un exemple avec la formalisation de la théorie des textons) que l’on peut rajouter l’invariance morphologique aux autres propriétés de l’algorithme de segmentation défini précédément.

1.5 Une application : la discussion d’une théorie psychophysique, la Théorie des Textons.

Certaines théories psychophysiques de la vision bas niveau sont axiomatisées. C’est le cas de la théorie des textons de Julesz ([Jul1, Jul2, JuKr]), qui tente de formaliser la façon dont la perception humaine est capable de discriminer des textures de manière préattentive en une trentaine de millisecondes. Pour comprendre le concept intuitif de texture, il suffit de regarder les paires de textures dans la figure 4. Même si on est bien incapable de décrire géométriquement ou statistiquement les deux textures, la discrimination des deux textures est très nette pour la perception humaine. Inversement, deux textures différentes peuvent être indiscriminables. C’est le cas pour la figure 5 où une région de ”10” est entourée de ”5”. Bien que ces formes soient individuellement faciles à distinguer, quand elles composent une texture, la discrimination devient très difficile. Ce sont des faits de ce type que Julesz a tenté de formaliser dans la théorie des ”textons”.

Pour résumer briévement la doctrine de Julesz, disons que Julesz suppose que la perception analyse par un filtrage local les courbures et orientations des formes présentes dans les textures, puis en effectue des moyennes locales. Ces moyennes s’appellent les ”densités de textons”. Donc à l’image u_{0} est associée une double suite d’images u_{i}(t, x) qui sont des densités de textons du type i à l’échelle t. Ensuite, un processus de segmentation agissant sur cette suite d’images (segmentation multicanaux) est effectué. Il permet le calcul des bords entre textures différentes et les régions homogènes.

La refonte de la théorie de Julesz dans le cadre axiomatique multiéchelle que nous avons développé plus haut est évidente.

1) On applique l’équation fondamentale (AMSS) à l’image : c’est la seule habilitée par ses invariances à calculer une courbure et une orientation multiéchelle.

2) Les densités de textons en sont déduites étant les seuls invariants intégraux calculés à partir de la courbure et de l’orientation multiéchelle.

3) Les densités de textons étant calculées, on effectue une segmentation multiéchelle de l’image. On obtient alors un dispositif rigoureux permettant de dire s’il y a discrimination ou pas pour la théorie des textons : on arrète l’algorithme de croissance de régions décrit au paragraphe précédent (Region Growing) quand il n’y a plus que deux régions. Si ces deux régions sont celles prédites par la théorie, il y a discrimination.

Par conséquent, la théorie des textons peut être testée numériquement, et le résultat est très différent de celui attendu. En effet, les premières tentatives de formalisation était liées au filtrage linéaire multiéchelle suivi de mécanismes nonlinéaires compatibles avec les connaissances physiologiques. Une des tentatives les plus réussies est due à Malik et Perona [MaPe] car leur algorithme réalise les mêmes performances que celles de la vision humaine dans le domaine de la discrimination de textures.

Mais en fait, une implémentation rigoureuse de la théorie des textons permet de réaliser l’hyperdiscrimination ! Julesz avait créé des paires de textures qui étaient indiscriminables en vision préattentive mais elles ne le sont pas pour l’algorithme que nous avons décrit (figure 5). Ceci permet d’affirmer que dans le cas de la discrimination de textures, l’analyse des textons est inhibée à une échelle beaucoup plus petite que leur taille. Ainsi, les axiomes de la théorie des textons (en particulier la notion de Δ-voisinage) peuvent être reformulés mathématiquement.

A. Le “théorème fondamental” de l’analyse d’images

Soit (T_{t+h,t}) avec t > 0 et h ≥ 0 une famille d’opérateurs agissant sur un espace de fonctions C^{2} à croissance lente telle que l’on ait la causalité, le principe de comparaison local et le principe de régularité suivant : ”Pour une image très régulière, l’évolution est régulière au moins pendant une durée h petite”. Précisons cette dernière hypothèse ; Considérons la fonction Q définie sur IR^{2} par y → Q_{x}(y) = \frac{1}{2}(A(y − x), (y − x)) + (p, (y − x)) + c où A est une matrice symétrique, p un vecteur de IR^{2} et c un réel. Considérons l’application h → (T_{t+h,t}Q_{x})(x). On suppose que cette application est dérivable en h = 0 et que F (A, p, c, x, t) = \frac{∂T_{t+h,t}Q_{x}}{∂h} (x) est continue par rapport à la matrice symétrique A.

On peut alors conclure qu’en tout point x où u(t, x) = T_{t}u(x) est de classe C^{2}, u vérifie

\frac{∂u}{∂t} = F (D^{2}u, Du, u, x, t)

et l’application A → F (A, p, c, x, t) est croissante.

Preuve : Par hypothèse u est régulière au voisinage du point x

u(t, y) = u(t, x) + (Du, y − x) + \frac{1}{2}D^{2}u(y − x, y − x) + o(|y − x|^{2})

Posons

Q_{ϵ}(y) = u(t, x) + (Du, y − x) + \frac{1}{2}D^{2}u(y − x, y − x) + ϵ.|y − x|^{2}

Alors on a la double inégalité suivante au voisinage de x

Q_{−ϵ}(y) ≤ u(t, y) ≤ Q_{ϵ}(y)

Les principes de causalité et de comparaison locale permettent d’écrire

(T_{t+h,t}Q_{−ϵ})(x) ≤ (T_{t+h,t}u(t))(x) ≤ (T_{t+h,t}Q_{ϵ})(x)

d’autre part, on a aussi

Q_{−ϵ}(x) = Q_{ϵ}(x) = u(t, x) = (T_{t,t}Q_{−ϵ})(x) = (T_{t,t}Q_{ϵ})(x)

on obtient ensuite en prenant la différence des deux relations précédentes et en divisant par h > 0 et en faisant tendre h vers 0

\frac{∂T_{t+h,t}Q_{−ϵ}}{∂h} (x) ≤ lim inf \frac{T_{t+h,t}u(, x) − T_{t,t}u(x)}{h} ≤

lim sup \frac{T_{t+h,t}u(, x) − T_{t,t}u(x)}{h} ≤ \frac{∂T_{t+h,t}Q_{ϵ}}{∂h} (x).

Le principe de régularité et la continuité de F (on regarde la limite lorsque ϵ tend vers 0) permettent de conclure. En reprenant le même raisonnement avec

Q_{A}(y) = \frac{1}{2}(A(y − x), (y − x)) + (p, y − x) + c Q_{B}(y) = \frac{1}{2}(B(y − x), (y − x)) + (p, y − x) + c

on obtient F (A, p, c, x, t) ≥ F (B, p, c, x, t) si A ≥ B.

Remarque : Pour la simplicité, nous avons montré que l’équation (A.M.S.S.) est vraie ponctuellement quand u est de classe C^{2}. En fait, on montre par les mêmes techniques si T_{t} vérifie les axiomes et si u est une solution uniformément continue alors elle est solution de viscosité de (A.M.S.S.) (voir [A.G.L.M.], [GiGo]). Inversement, une équation du type précd́ent a toujours une solution de viscosité [CIL, EvSp, CGG].

B. Morphologie Mathématique Affine invariante.

Quelle est la relation du modèle AMSS avec les opérations de la morphologie mathématique? Définissons d’abord une ”ouverture affine” et une ”fermeture affine” et posons

(O_{t}u)(x) = inf_{B∈C_{t}(x)}(sup_{y∈B}(u(y))

(F_{t}u)(x) = sup_{B∈C_{t}(x)}(inf_{y∈B}(u(y))

où C_{t} est l’ensemble de tous les ensembles convexes dont la surface est t^{3/2}. Si C_{t} est réduit à la boule de rayon t, les définitions ci-dessus correspondent aux fermetures et ouvertures classiques de la morphologie mathématique [Serr]. Lasry et Morel [LaMo] ont prouvé que si on itére alternativement la version affine invariante de ces deux filtres avec un petit temps δt, alors le filtrage asymptotique vérifie l’équation fondamentale (AMSS)!

Ces idées sont utilisées pour définir des algorithmes aussi invariants que le modèle continu (voir [CEGLM]). C’est un algorithme de ce type qui est utilisé dans la figure 2.

En résumé, du point de vue algorithmique, on peut obtenir l’équation (A.M.S.S.) en itérant des filtres morphologiques : l’ouverture et la fermeture affines. Les filtres morphologiques ont été introduits en traitement d’images par Matheron, et l’école de Fontainebleau à cause de leurs remarquables propriétés d’invariance, qui les rendent indispensables pour l’analyse de formes. L’équation fondamentale du traitement d’images, permet donc d’effectuer une infinité d’inf-sup-convolution d’un seul coup avec des ellipses de toutes excentricités. (D’où l’invariance affine !) La théorie unifiée peut donc aussi être appelée morphologie mathématique affine-invariante et elle aura un impact important sur les applications (robotique, reconnaissance des formes, télédetection, ...).

Références:

[AGLM1] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomes et équations fondamentales du traitement d’images. Report 9216, CEREMADE. Université Paris Dauphine, 1992.

[AGLM2] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomatisation et nouveaux opérateurs de la morphologie mathematique. C.R. Acad. Sci. Paris, t.315, Série I, p. 265-268, 1992.

[ALM] L. Alvarez, P.L. Lions and J.M. Morel. Image selective smoothing and edge detection by nonlinear diffusion (II). SIAM Journal of numerical analysis 29, pp 845-866. 1992.

[Ang] S. Angenent. Parabolic equations for curves on surfaces I, II. University of Wisconsin-Madison Technical Summary Reports, 89-19, 89-24, 1989.

[AsBr] H. Asada and M.Brady. The curvature primal sketch. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 8, No. 1, January. 1986.

[Bar] G. Barles. Remarks on a flame propagation model. Technical Report No 464 INRIA Rapports de Recherche, December 1985.

[BaGo] C. Ballester et M. Gonzalez. Affine invariant multiscale segmentation by variational method. à paraître.

[BaSo] G. Barles and P.E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equation, Asymp. Anal., to appear.

[BrFe] C. Brice and C. Fennema. Scene analysis using regions. Artificial Intelligence, 1, 205-226, 1970.

[BrMa] W. Brockett and P. Maragos. Evolution Equations for continuous-scale morphology. ICASSP March 23-26, San Francisco, 1992.

[CEGLM] T. Cohignac, F. Eve, F. Guichard, C. Lopez, J.M. Morel. Numerical Analysis of the fundamental equation of image processing. Soumis.

[CGG] Y-G. Chen, Y. Giga and S. Goto. Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations, Preprint, Hokkaido University, 1989.

[CIL] M.G. Crandall, H. Ishii and P.L. Lions. User’s guide to viscosity solution of second order partial differential equation. Preprint, CEREMADE, 1990.

[DMS] G. Dal Maso, J.M. Morel and S. Solimini. Une approche variationnelle en traitement d’images : résultats d’existence et d’approximation. C.R. Acad. Sci. Paris, t.308, Série I, p. 549-554, 1989.

[DiKo] F. Dibos and G. Koepfler. Propriété de régularité des contours d’une image segmentée. C.R. Acad. Sci. Paris, t. 313, Série I, p. 573-578, 1991.

[Enn] J. Enns. Seeing textons in context. Perception and Psychophysics, 39 (2), 143-147, 1986.

[EvSp] L.C. Evans and J. Spruck. Motion of level sets by mean curvature, I., preprint.

[FRKV] L. Florack, B. ter Haar Romeny, J.J. Koenderink and M. Viergever. Scale and the differential structure of images. Image and vision computing. Vol 10, n. 6 July/August 1992.

[FRKV] L. Florack, B. ter Haar Romeny, J.J. Koenderink and M. Viergever. General Intensity transformations and Second Order Invariants. Procedings of the 7th Scandinavian Conference on Image Analysis. Aalborg, august 13-16, 1991.

[FMZ] D.Forsyth, J.L.Mundy, A.Zisserman, al. Invariant descriptors for 3-D object recognition and Pose. IEEE Transaction of pattern analysis and machine intelligence. Vol 13 N.10 1991.

[GaHa] M. Gage and R.S. Hamilton. the heat equation shrinking convex plane curves. J. Differential Geometry 23, pp. 69-96, 1986.

[Gray] M. Grayson. The heat equation shrinks embedded plane curves to round points. J. Differential Geometry 26, pp. 285-314, 1987.

[GeGe] S. Geman and D. Geman. Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE PAMI 6, 1984.

[HaSh] R.M. Haralick and L.G. Shapiro. Image segmentation techniques. Computer Vision Graphics and Image Processing, 29, 100-132, 1985.

[JuBe] B. Julesz and J.R. Bergen. Textons, the fundamental elements in preattentive vision and perception of textures. Bell System Tech. J. 62(6) 1619-1645, 1983.

[JuKr] B. Julesz and B. Kroese. Features and spatial filters. Nature, 333, 302-303, May 1988.

[Jul1] B. Julesz. Textons, the elements of texture perception, and their interactions. Nature, 290, 12 March 1981.

[Jul2] B. Julesz. Texton gradients : the texton theory revisited. Biol. Cyber. 54, 245-251, 1986.

[Kim] B.B. Kimia. Toward a computational theory of shape. Ph. D. Dissertation, Department of Electrical Engineering, McGill University, Montreal, Canada, August 1990.

[KLM] G. Koepfler, C. Lopez, J.M. Morel. A multiscale algorithm for image segmentation by variational method . Accepted at SIAM journal of Numerical Analysis.

[KMS] G. Koepfler, J.M. Morel and S. Solimini. Segmentation by minimizing a functional and the “merging” methods. Proceedings ‘GRETSI Colloque’ in Juan-les-Pins, September 1991 (France).

[KoDo1] J.J. Koenderink and A.J. van Doorn. Dynamic shape. Biol. Cyber. 53 : 383-396, 1986.

[KoDo2] J.J. Koenderink and A.J. van Doorn. Representation of Local Geometry in the Visual System. Biol. Cyb. 55 : 367-375, 1987.

[Koe1] J.J. Koenderink. Solid Shape. MIT Press, cambridge, MA 1990.

[Koe2] J.J. Koenderink. The brain a geometry engine. Psychol. Res., Vol 52 (1990) pp. 122-127.

[KTZ] B.B. Kimia, A. Tannenbaum, and S.W. Zucker. On the evolution of curves via a function of curvature, 1: the classical case. J. of Math. Analysis and Applications. Vol. 163, No 2, Jan. 15 1992.

[KWT] M. Kass, A. Witkin and D. Terzopoulos. Snakes: active contour models. 1st Int.Comp.Vis.Conf. 1987, IEEE n^{∘}777.

[LaMo] J.M. Lasry and J.M Morel. A monotone consistent theoretical scheme for the fundamental equation of image processing. preprint.

[LoMo] C. Lopez, J.M. Morel. Axiomatisation of shape analysis and application to texture hyperdiscrimination. Proceedings of the Trento Conference on Surface Tension and Movement by Mean Curvature. De Gruyter Publishers. Berlin 1992.

[LSW] Y.Lamdan, J.T.Schwartz and H.J.Wolfson. Object recognition by affine invariant matching. In Proc. CVPR 88, 1988.

[MaMo1] A. Mackworth and F.Mokhtarian. Scale-Based description and recognition of planar curves and two-dimensional shapes. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 8, No. 1, January. 1986.

[MaMo2] A. Mackworth and F.Mokhtarian. A theory of multiscale, curvature-based shape representation for planar curves. IEEE Trans. Pattern Anal. Machine Intell. 14, pp 789-805, 1992.

[MaPe] J. Malik and P. Perona. Preattentive texture discrimination with early vision mechanisms. Journ. of the Opt. Society. of America. A, 923-932, vol. 7, num 5, 1991.

[Mara1] P. Maragos. Tutorial on advances in morphological image processing and analysis. optical engineering, july 1987, Vol. 26 No. 7.

[Marr1] D. Marr. Analysing natural images : a computational theory of texture vision. Cold Spring Harbor Symposium on Quantitative Biology, XL, 647-662, 1976.

[Marr2] D. Marr. Vision. Freeman and Co. 1982.

[Math] G. Matheron. Random Sets and Integral Geometry. John Wiley, N.Y. 1975.

[MoSo1]J.M. Morel and S. Solimini. Segmentation of images by variational methods: a constructive approach. Rev. Matematica de la Universidad Complutense de Madrid, Vol.1, n^{∘}1,2,3, 169-182, 1988.

[MoSo2] JM. Morel and S. Solimini. Segmentation d’images par méthode variationnelle : une preuve constructive d’existence. C. R. Acad. Sci. Paris Série I, 1988.

[MoSo] J.M. Morel, S. Solimini. Variational methods in image segmentation. A paraître, Birkhauser, 1993.

[MuAl] J.L. Muerle and D.C. Allen. Experimental evaluation of techniques for automatic segmentation of objects in a complex scene. In Pictorial Pattern Recognition, G.C. Cheng and al. eds., pp 3-13, Thompson, Washington, 1968

[MuSh1] D. Mumford and J. Shah. Optimal Approximations by Piecewise Smooth Functions and Associated Variational Problems. Communications on Pure and Applied Mathematics. vol. XLII No.4, 1989.

[MuSh2] D. Mumford and J. Shah. Boundary detection by minimizing functionals. Image Understanding, 1988, Ed. S. Ullman and W. Richards.

[OsSe] S. Osher and J. Sethian. Fronts propagating with curvature dependent speed : algorithms based on the Hamilton-Jacobi formulation. J. Comp. Physics, 79, 12-49, 1988.

[Pav] T. Pavlidis. Segmentation of pictures and maps through functional approximation. Comp.Gr. and Im.Proc. 1, 360-372, 1972.

[PaLi] T. Pavlidis and Y.T. Liow. Integrating Region Growing and Edge Detection.

Proc. of the IEEE Conf. Comp. Vision and Patt. Recognition 1988.

[PeMa] J. Malik and P. Perona A scale space and edge detection using anisotropic diffusion. Proc. IEEE Computer Soc. Workshop on Computer Vision, 1987.

[SaTa1] G. Sapiro and A. Tannenbaum. On affine plane curve evolution. Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 821. February 1992, submitted.

[SaTa2] G. Sapiro and A. Tannenbaum. Affine shortening of non-convex plane curves . Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 845. February 1992, submitted.

[Ser] J. Serra Image analysis and mathematical morphology, Vol 1. Academic Press 1982.

[Tre] A. Treisman. Preattentive processing in vision. Comput. Vision, Graphics, Image Process., 31(2), 156-177, 1985.

[VoPo1] H. Voorhees and T. Poggio. Detecting textons and texture boundaries in natural images. Proc. Intern. Conf. Comput. Vision, IEEE, 250-258, 1987.

[Wit] Witkin, A. P. Scale-space filtering. Proc. of IJCAI, Karlsruhe 1983, 1019-1021.

[Yui] A. Yuille. The creation of structure in dynamic shape. Proc. of the second international conference on computer vision, Tampa 1988, pp. 685-689.

[YuPo] A. Yuille and T. Poggio. Scaling theorems for zero crossings. IEEE Transactions on Pattern Analysis and Machine Intelligence, 8, January 1986.

[Zuc] S.W. Zucker. Region growing: Childhood and Adolescence (Survey). Comp. Graphics and Image Proc. 5, 382-399, 1976. fig.1 Analyse Multiéchelle.

Auteurs : L. Alvarez - F. Guichard.

L’analyse multiéchelle associe à une image u(x) une succession d’images de plus en plus lisses u(x, t). u(x, 0) est l’image de départ u(x), et t est l’échelle. u(x, t) vérifie l’E.D.P. ([AGLM])

\frac{∂u}{∂t} = |Du|(t.curv(u))^{\frac{1}{3}}

Du est le gradient de u, curv(u) la courbure de la ligne de niveau de u passant par x. Colonne de gauche : En haut : l’image originale 128 × 128. Au milieu : l’image originale bruitée avec un bruit d’écart type 80. En bas : les valeurs de 70 % des pixels sont changées aléatoirement. Colonne de droite : Extrait de l’analyse multiéchelle de chacunes des images, (échelle 10, temps CPU : 21s). L’analyse multiéchelle des images est un passage du local au global ; d’où une très grande stabilité face au bruit. fig.2 Reconnaissance des formes projectivement invariante grâce à l’équation fondamentale du traitement d’images.

Auteurs: Cohignac - Eve - Guichard - Lopez - Morel.

Quelle implémentation numérique faire d’une équation qui doit servir à la reconnaissance de formes? Les exigences d’une telle implémentation sont d’un type nouveau. On désire, après un grand nombre d’itérations sur une donnée aussi complexe qu’une image, avoir un résultat qui permette la reconnaissance. L’algorithme d’analyse multiéchelle doit donc vérifier strictement les mêmes invariants que ceux vérifiés par l’équation fondamentale. Cette planche teste, pour une application au Ministère de la Défense, l’invariance projective. Une figure subit trois anamorphoses de valeurs propres (1/2,2), (1/3,3) et (1/4,4) (voir le quadrant 1) avant d’être analysée par l’équation fondamentale (quadrants 5,6,7,8). L’anamorphose inverse permet de vérifier une parfaite invariance de l’algorithme (quadrants 2,3,4). fig.3 Segmentation d’image à l’aide de l’energie de Mumford et Shah. Auteurs: Koepfler - Morel. Résultat obtenu lorsqu’on décide de garder 100 régions. On peut voir l’image originale, la carte des bords, la ”bande dessinée” et la superposition de l’image originale et de la carte des bords. fig.4 Segmentation de textures naturelles.

Auteurs: Koepfler - Lopez - Morel.

A gauche : deux textures de Brodatz, à droite : trois différentes textures de Brodatz (image 256 x 256). Dessous : estimation des bords (Segmentation d’une image filtrée en 18 canaux différents). Il faut noter que la même texture est répétée deux fois à gauche. Bien que qu’une frontière soit apparente, l’algorithme reconnait l’identité ”en moyenne” des textures de haut en bas et segmente correctement. fig.5 Vérification de la théorie de Julesz et algorithmes hyperdiscriminants.

Auteurs: Lopez - Morel.

Julesz avait synthétisé des textures qui sont indiscriminables pour le système visuel humain et il pensait que cette indiscriminabilité était conforme à sa théorie des ”textons”. En réalité, l’équation fondamentale du traitement d’images, conforme à sa théorie, permet de discriminer ces textures! Premier quadrant (512x512): une paire de textures de Julesz indiscriminables en vision préattentive. Deuxième quadrant: image obtenue après évolution par l’équation à l’échelle de 7 pixels. Troisième et quatrième quadrants: courbure négative et positive à l’échelle de 7 pixels. Cinquième quadrant: segmentation obtenue en appliquant l’algorithme de croissance de régions aux canaux de courbures multiéchelle (Le dernier quadrant illustre la nécessité d’extraire une courbure à plusieurs échelles car la courbure négative de l’image originale ne permet pas de conclure.).

A multiscale algorithm for image segmentation by variational method.

This chapter has been written in collaboration with Georges Koepfler and Jean Michel Morel and has been accepted at SIAM Journal of Numerical Analysis.

Abstract. Most segmentation algorithms are composed of several procedures: split and merge, small region elimination, boundary smoothing, …, each depending on several parameters. The introduction of an energy to minimize leads to a drastic reduction of these parameters. We prove that the most simple segmentation tool, the “region merging” algorithm, made according to the simplest energy, is enough to compute a local energy minimum belonging to a compact class and to achieve the job of most of the tools mentioned above. We explain why “merging” in a variational framework leads to a fast multiscale, multichannel algorithm, with a pyramidal structure. The obtained algorithm is O(n ln n), where n is the number of pixels of the picture. We apply this fast algorithm to make grey level and texture segmentation and we show experimental results.

2.1 Introduction.

The aim of this paper is to describe a fast and universal image segmentation algorithm. Properties of this algorithm, are :

It is multiscale and pyramidal. In other terms, it will not only compute one segmentation, but a hierarchy of segmentations from fine to coarse scales. Moreover the coarser segmentation will be deduced from the finer by “merging” operations, with a pyramidal structure for the computation. This corresponds to Marr’s remark [Marr] that textures “live” at several scales. Therefore, a discrimination algorithm should give different kinds of segmentations which depend on the scale.

As a consequence of this pyramidal structure, the computation time is in practice proportional to the size of the datum. Moreover, this algorithmic structure will make it accessible to efficient hardware implementation.

The algorithm is universal, that is, does not depend on any a priori knowledge on the statistics of the image. Texture discrimination is achieved according to a universal criterion depending on very few parameters. More precisely, a picture is defined by a certain number of “channels” (grey level, colour levels, integrodifferential channels obtained by fast wavelet transform). In our segmentation algorithm, the only parameters are the weights attached to each channel, that is, the importance given to each channel as a segmentation criterion. However these parameters can be fixed once for all with reliable results.

The algorithm has been constructed by making a synthesis of several theories : the textons theory of Julesz [Ju], the energy methods in image segmentation introduced by Geman and Geman [GemG], Blake and Zisserman [BlakZ], Mumford and Shah [MumS1] , the Wavelet transform theory of Meyer [Me] , Mallat [Mall] and Cohen [Coh] , which unifies the theory of recursive filtering and pyramidal schemes.

The segmentations provided by the algorithm will be proved to have a large range of good topological and numerical properties, including compactness of the set of approximate solutions, convergence of minimizing sequences made of finer and finer solutions, smoothness of the locally optimal solutions, completeness of the multiscale representation, a priori estimates on the size of the regions of the segmentations.

2.2 General principles of segmentation devices.

2.2.1 Formalization.

We define an image g as a scalar function, defined on the image domain Ω (generally a rectangle). The function g may also be vectorial, in the case where it has several channels for characterizing textures, histograms, colors, … That does not change anything in the theorems and proofs that we later state. The only hypothesis is that these channels have been defined in order to be good indicators of the similarity – or difference – of points of the picture, and therefore good indicators for the autosimilarity of regions. We seek for a segmentation, that is, a partition of this rectangle into a finite set of regions, each of which corresponds to a part of the image where g is as constant as possible. Moreover, we wish to compute explicitly the region boundaries and of course control their regularity and location. More precisely we will adopt the following principles.

The first principle is that the boundary detection problem must follow some universal rules. In other terms, we admit the possibility of a universal boundary detection device, definable and analyzable independently from the kind of channels (grey level, colour, texture, …) to be used as input for the segmentation ([BecSI, MaliP2, Ju]). This first principle allows us to try to get a complete mathematical understanding of the segmentation problem, by considering the grey level segmentation. It is the simplest case of boundary detection, and it eliminates every early discussion about the concept of texture. (By “boundary in an image”, we mean the boundary in the topological sense : it is the boundary of an homogeneous region of the image. Thus boundaries are different from edges obtained by some local filtering.)

A second principle which we adopt in the following is that an algorithm for boundary detection must be scale and space invariant. By space invariance, we mean that all points of the analyzed picture will be treated on the same way, and therefore a segmentation device must be translation and rotation invariant. Since boundaries may be present at every “scale”, the scale invariance only means that we shall consider multiscale segmentation algorithms, depending on a scale parameter λ which we do not try to estimate and leave to the user’s choice.

Finally we shall adopt a principle without which no discussion about segmentation can even start, and which we call comparison principle. It states that given two different segmentations of a datum, we are always able to decide which of them is considered as better than (or equivalent to) the other. Thus we assume the existence of some total ordering over all possible segmentations, and this can be simply achieved only if this ordering is reflected by some real functional E such that if E(K_{1}) < E(K_{2}), then the segmentation K_{1} has to be considered “better” than the segmentation K_{2}. For instance, this principle is verified by segmentation devices based on some Gibbs energy functional like in [GemG]. It is not verified by the region growing methods based on thresholds [Pav1, Pav2, PavL, Zu], nor by the edge detection devices [Marr, MaliP1].

Since, by the comparison principle, all of these criteria must be taken into account in the functional E, we see that this functional necessarily contains terms which control :

- the autosimilarity of each region with respect to the chosen channels. Those channels must be as constant as possible on each region.

- the size, location and regularity of the boundaries.

Denote by K the whole boundary set of the segmentation. According to these principles, a natural energy functional for segmentation will contain two terms, a two dimensional term for the autosimilarity of the regions which will roughly speaking measure the variance of g on each connected component of Ω \ K and a one dimensional one for controlling the length, and eventually the adequacy of location of the boundaries. Such a generic justification is developed e.g. in [MumS1, MumS2]. By the space invariance principle, those terms will be integral terms, the two dimensional one being an integral with respect to the Lebesgue measure and the one dimensional an integral with respect to the Hausdorff measure or “length” along the boundaries. Clearly, the weights given to the terms of the functional will be left to our choice. Consider for example the energy used by Mumford and Shah :

E(u, K) = ∫_{Ω\K} |∇u|^{2}dx + ∫_{Ω}(u − g)^{2}dx + ∫_{K} dσ ,

Both first integrals are bidimensional and the third one is with respect to the uniform measure dσ supported by K. This energy means that in a “good” segmentation (u, K), the curves of K should be the boundaries of homogeneous regions in the image and u a sort of mean or, more generally, a regularized version of g in the interior of such areas (as illustration see the cartographical application in figure 2.1, u is choosen to be piecewise constant). The third term gives a control on the length, and therefore on the regularity of the boundaries of K ( see below ). This kind of functional represents of course a compromise between accuracy of the regions and parsimony of the boundaries: in case of grey level segmentation, as noticed in Zucker [Zu] or Haralick and Shapiro [HarS], a pure region growing would simply put together the pixels with similar grey levels. But this generates very nonsmooth boundaries and therefore also “small” or “thin” regions. So if no control is made on the boundaries, one needs additional criteria and thresholding to achieve a decent segmentation. A functional like the above functional is designed to avoid this kind of mixed methods : one hopes having all the criteria put together in the same functional. For instance, the functional used for the so called “snakes” [KasWT, AmiTW] uses a more sophisticated third term for controlling both location (the boundary, or “snake”, is forced to be close to edges) and smoothness of the boundary. Many early functionals in image analysis had only two dimensional energy terms (which were coupled with some thresholding criteria [Pav1, Zu]).

2.2.2 Choice of properties and algorithms for segmentation.

The algorithmic description given by Pavlidis and Liow [PavL] is an excellent illustration of the level of sophistication to which the segmentation methods have arrived in Computer Vision, by superposing several tricky devices arisen in the two last decades. Now, the question which is naturally raised by the coexistence of energy functionals on one side, and of segmentation devices on the other is whether they match in any aspect, and how.

Our purpose is to classify the properties which are sought by those devices and to decide which of them are basic, and which can be deduced. Beyond being a way to attain certainty, mathematical proofs are a powerful tool for obtaining such a classification of properties. Once this hierarchy among tools has been established, one may hope to simplify the algorithms and to make them completely transmissible and reproducible. Our aim in this paper is to apply this method and to classify terms like “small or thin region elimination”, “splitting”, “merging”, “energy minimizing”, “region growing”, which have proved to be so necessary in most discussions about segmentation, all this following the ideas developed above. In the following we shall focus on the simplest model of Mumford and Shah [MumS1]. According to this model the segmentation (u, K) should be obtained by minimizing the functional

E(u, K) = ∫_{Ω\K}(u − g)^{2}dxdy + λℓ(K),  (2.1)

where K is a union of boundaries in Ω with Hausdorff length ℓ(K), and u is piecewise constant on Ω \ K. The constant λ represents the scale parameter of the functional and measures the amount of boundary : if λ is low, a lot of boundaries are allowed and we get a “fine” segmentation. As λ increases, the segmentation gets coarser and coarser. Similar functionals have been recently introduced by several authors in order to modelize physical phenomena like phase transition [FoT] and liquid crystals [DeG, AmbDeG, CarDeGL] .

In what follows, we shall focus on functional (2.1) because it is the simplest one containing most of the algorithmic difficulties of such functionals. Moreover, it is the only one for which a complete mathematical analysis is available. Indeed, Mumford and Shah have proved in [MumS1] that the boundary set obtained by this functional has the following geometric property : either the points of K are regular (at least C^{1}) or the singular points are of two types, namely, triple points where three branches meet with 120^{∘} angles and boundary points where K meets the boundary of Ω at a 90^{∘} angle. Moreover, the boundaries of the segmentation verify variational inequalities of the kind

(u^{+} − g^{−})^{2} − (u^{−} − g^{−})^{2} ≤ curv(x) ≤ (u^{−} − g^{+})^{2} − (u^{+} − g^{+})^{2},  (2.2)

where u^{+} and u^{−} denote the values of u(x) on both sides of the boundary.

In [MorS1, MorS2] Morel and Solimini have given an elementary and constructive proof of the above theorem. Of course, the fact that the crossing points must be ternary proves that such a functional gives no hope of accurate segmentations. But this is, in our view, secondary, because one can improve the functional by adding or modifying boundary terms (of the kind used for snakes) or surface terms (for instance by imposing that u is a linear, quadratic or cubic spline, …). This would not change the nature of the proofs and methods which we discuss here and which apply to most energy functionals mentioned above. Note however, that the above mentioned general functional of Mumford and Shah has a different structure. As proved in [BlatM, CarLPP, MumS1], it allows open boundaries and therefore the merging methods which will be discussed here are not adequate. Moreover, the proof of the “Mumford and Shah conjecture”, which states that this functional would lead to a finite set of smooth boundaries is not achieved. The best available result [DalMMS1, DalMMS2] is that the minimizing boundaries (which in some weak sense exist [Amb]) can be approximated by a finite number of curves. No regularity result is available, and we see it difficult to comment now mathematically the first and interesting numerical simulations [BlakZ] of this problem. An attempt to do so is however presented in [DibK]. Let us therefore return to our “simplest functional” (2.1). It is well known that functionals of this kind may have many local minimizers. In the case of the functional (2.1), it is easy to give explicit examples. Moreover, a reasonable conjecture is that to find a minimizer, or even to prove a posteriori that some computed local minimum is global leads to a NP-complete problem.

One has thus to choose between two strategies :

• The global minimization by simulated annealing methods, which leads to huge computations, but ensures that in some asymptotic sense, the global minimum is attained [GemG].

• Another way is then to define some concept of local minimum which should be more accessible to fast computations and verify the same properties which were seeked for the global minimum.

For instance, in the case of simulated annealing, a recent tendency is to define some faster, parallel processes, but which do not pretend to find a global minimum anymore [Az1, Az2]. The homotopy method of Blake and Zisserman [BlakZ] also seeks for “good” local minima. (It seems to be close to a Γ-convergence device suggested by the mathematical school of De Giorgi [AmbT] and developed by Tom Richardson in his Phd. Dissertation at MIT [Ric]).One needs a theory to do that in the deterministic approach to which the segmentation devices mentioned above belong. In the next section we shall give some basic definitions and notations concerning our “simplest” energy functional. In sections 2.3.1 and 2.3.2, we give a precise account of topological properties and the proofs of the announced results concerning segmentations obtained by “merging”. The definition and properties of what we could call the “simplest possible” recursive merging algorithm which will be described in section 2.4.1, there will also be shown results obtained by running this algorithm on grey level images and images containing “textures”.

2.3 Compactness properties of segmentations obtained by “merging”.

Before presenting the main result, let us fix a functional and topological glossary.

• The function g. It is a bounded measurable real function on a rectangle Ω (more generally for g : Ω ↦ IR^{N}, N ≥ 2 the same result holds). Without loss of generality we assume that |g(x, y)| ≤ 1.

• The energy E(K) : note that given the boundaries K, the corresponding minimal u is completely defined by the fact that its value on each connected component of Ω \ K is equal to the mean value of u on this connected component. Thus we shall always assume in the following that to each K is associated this unique u. Therefore we shall write E(K) instead of E(u, K).

• Regions, or connected components of Ω \ K: we shall denote them by (O_{i})_{i}.

• Common boundary of two regions O_{i} and O_{j} : we denote it by ∂(O_{i}, O_{j}). It is contained in K. If i = j, ∂O_{i} denotes the boundary of O_{i}.

• Two dimensional measure of O_{i} : denoted by |O_{i}|.

• Isoperimetric inequality in IR^{2} and Ω : denote by O a region in IR^{2}. Then one

has

ℓ(∂O) ≥ 2\sqrt{π} · \sqrt{|O|} .

In the case of a domain Ω with smooth boundary, the same kind of inequality holds for the relative boundary of O in Ω, ∂O ∩ Ω, with a smaller constant C :

ℓ(∂O ∩ Ω) ≥ C · \sqrt{|O|} .

(Since there is no ambiguity, Ω being fixed, we still denote by ∂O = ∂O ∩ O the relative boundary of O in Ω.) By ℓ(∂O) we mean the length of the boundary of O. For a general definition of the length of the boundary (“perimeter”) of a region in the plane, see [Fed, §4.5.1] or [Si-L, §14].

Definition 1 A segmentation K will be called normal if every subsegmentation K′ of K verifies E(K′) > E(K).

By a subsegmentation of K, we mean a segmentation obtained by merging an arbitrary number of adjacent regions. The normality of a segmentation means that no merging operation of a set of regions of the segmentation K can decrease the energy. This definition implies in particular that if K is normal, no boundary between two regions of the segmentation may be removed without increasing the energy E.

Definition 2 A segmentation K will be called 2-normal if for every pair of regions O_{i} and O_{j}, the new segmentation K′ obtained by merging these regions verifies E(K′) > E(K).

2-normality, with the same definition, but under the name of “optimality” has been introduced by Pavlidis [Pav1].

We shall consider only segmentations having the following properties, which are easy to check for computationally defined segmentations.

a) The number of regions is finite. In other terms, Ω \ K has a finite number of connected components. b) No region has internal boundaries. In other terms, the interior of the closure of

each region O is equal to the interior of O. Indeed, if this is not the case, one can remove the internal boundaries without increasing the energy. This property is called the “1-normality” (compare to the definitions above).

c) The segmentations are made of piecewise affine curves. (This assumption can be changed into piecewise C^{1}, or piecewise Lipschitz without changing the proofs.)

Let us introduce the framework in which we will study the topological properties of affine 1-normal segmentations.

• A curve c from [0, 1] into Ω has tips c(0) and c(1), all other points of the range of c are interior points.

• Geometrical support of a curve c : denoted by S(c), it is the range of the curve c. For several curves it is the union of the ranges. For instance K = ⋃_{i} S (c_{i}).

• A segmentation is said to be piecewise affine – resp. C^{1} – if the corresponding c_{i} are piecewise affine – resp. C^{1} –. A polygon is a connected region of an affine segmentation.

• Geometrical crossings : all the points of K where either a curve meets the interior of another, either three curves at least have a common tip, or a curve meets ∂Ω.

• Geometrical curve : subset of K which is in the geometrical support of a curve c such that its boundary is contained in the set of geometrical crossings and that no one of its interior points is a geometrical crossing.

• A Jordan curve is a continuous curve, such that for all σ, σ′ ∈]0, 1[ : c(σ) ≠ c(σ′) if σ ≠ σ′; if c(0) = c(1) the Jordan curve is said to be closed.

2.3.1 Topological properties of affine 1-normal segmentations.

Let us recall a classical result on Jordan curves (for a proof see for example [Ale, Th])

Lemma 1 Every closed Jordan curve c divides the plane in exactly two connected components, one is bounded “enclosed by c” and one unbounded.

We will say that a closed Jordan curve c encloses a part of the plane if this part is included in the bounded connected component delimited by c.

By using the preceding lemma one easily proofs the following (see[MorS1, MorS2]).

Lemma 2 Take an affine 1-normal segmentation K having more than one region, then there are two regions enclosed, O and O′, such that ∂(O, O′) is a Jordan curve.

Lemma 3 Let K be an affine 1-normal segmentation with α regions. Then K is the union of α − 1 affine Jordan curves having no common segments.

Proof. Let c_{1} be an affine Jordan curve as given by lemma 2 . Let us define a new segmentation K_{1} obtained by removing the segments of c_{1}. We merge exactly two regions of K. The segmentation K_{1} is also affine and 1-normal, so we can iterate the process and define a series of segmentations K_{i} such that every segmentation has a region less than the preceding one and has been obtained by removing an affine Jordan curve. The process stops if there is only one region left, Ω.

Lemma 4 Let α be the number of regions of an affine 1-normal segmentation, β the number of geometrical curves and γ the number of geometrical crossings. Then

γ ≤ 2 · (α − 1) and β ≤ 3 · (α − 1) − 2 .

Proof. Consider the proof of the preceding lemma (3). The curve c_{1} contains at most two geometrical crossings of K_{1} : its tips if it is an open Jordan curve. So passing from K = K_{0} to K_{1} eliminates at most two geometrical crossings. By iterating the process the first inequality follows. Let us now analyze how the number of geometrical curves decreases when passing from K_{0} to K_{1}. If c_{1} is a closed Jordan curve it is the only geometrical curve to disappear. If c_{1} is an open Jordan curve there are two geometrical crossings at its tips. If one of these crossings meets more than four segments of K the other geometrical curves are not modified. But if there are only three segments at this crossing it will vanish when c_{1} is removed and the other two geometrical curves will merge into one. Removing c_{1} decreases β of at most 3 and iterating the process yields the second inequality.

2.3.2 Estimates on segmentations obtained by merging.

Lemma 5 Let (u, K) be a 2-normal segmentation. Then every pair of regions O and O′ verifies

ℓ(∂(O, O′)) ≤ \frac{osc(g)^{2}}{λ} · min(|O|, |O′|) ,

where osc(g) = sup(g) − inf(g) is the oscillation of g.

Proof. Let K′ = K \ ∂(O, O′). By 2-normality we have

0 ≤ E(K′) − E(K) = ∫_{O∪O′}(u_{ij} − g)^{2} − ∫_{O}(u_{i} − g)^{2} − ∫_{O′}(u_{j} − g)^{2} − λℓ(∂(O, O′)),

where u_{i}, u_{j} and u_{ij} are the mean values of g on O, O′ and O ∪ O′. We get, supposing for example that |O| ≤ |O′|,

λ · ℓ(∂(O, O′)) ≤ ∫_{O} [(u_{ij} − g)^{2} − (u_{i} − g)^{2}]
≤ min(|O|, |O′|) osc(g)^{2}

Lemma 6 For every region O of a 2-normal segmentation, denote by N(O) the number of neighbouring regions. Then

N(O) ≥ \frac{Cλ}{\sqrt{|O|} osc(g)^{2}} ,

where C is the isoperimetric constant in Ω.

Proof. Call O_{j} a neighbouring region of O. The fact that O and O_{j} cannot be merged without increasing the energy E implies that

λ ℓ (∂ (O, O_{j})) ≤ |O| osc(g)^{2}.

Thus by adding these inequalities for all neighbours of O, we obtain

λ ℓ(∂O) ≤ N(O)|O| osc(g)^{2}.

We conclude by applying to O the isoperimetric inequality in Ω.

Lemma 7 Let α be the number of regions of a 2-normal affine segmentation. Then

α ≤ \frac{288 |Ω| osc(g)^{4}}{C^{2}λ^{2}} .

Proof. The union of all regions O_{i} is equal to Ω and therefore ∑_{i}|O_{i}| = |Ω|. Thus

the number of O_{i} verifying |O_{i}| ≤ \frac{2}{α} |Ω| is greater than \frac{α}{2}.

Let us now apply lemma 6 to all of these O_{i}s. Each one of them has at least

\frac{Cλ}{\sqrt{|O_{i}|} osc(g)^{2}} ≥ Cλ · \sqrt{\frac{α}{2 |Ω|}} · \frac{1}{osc(g)^{2}}

neighbouring regions. Consequently the number β of common boundaries, and therefore of geometrical curves, verifies

β ≥ \frac{α}{4} · Cλ \sqrt{\frac{α}{2 |Ω|}} · \frac{1}{osc(g)^{2}} = Cλ2^{− \frac{5}{2}} · \frac{α^{\frac{3}{2}}}{\sqrt{|Ω|}} · \frac{1}{osc(g)^{2}} .

By lemma 4, we have

β ≤ 3α,

this implies that

Cλ2^{− \frac{5}{2}} \frac{α^{\frac{3}{2}}}{\sqrt{|Ω|} osc(g)^{2}} ≤ 3α

and thus we obtain

α ≤ \frac{288 |Ω| osc(g)^{4}}{C^{2}λ^{2}} .

Remark 1 This is the same estimate (with a bigger constant) than in Mumford and Shah [MumS1, Th.5.2]. However, in the mentioned paper, this estimate is obtained for a global minimum. The fact that we get this estimate in the case of 2-normal segmentations indicates an analogy of structure between global minima and this kind of local minimum.

Corollary 1 The set of 2-normal affine segmentations has the following compactness property : for every sequence K_{n} of such segmentations, there exists a subsequence converging to a segmentation K such that

E(K) ≤ lim inf_{n} E(K_{n}).

K is not necessarily 2-normal, but has anyway a 2-normal subsegmentation with still less energy.

Proof. The proof of the announced compactness property is based on the fact that the number of edges of any 2-normal segmentation is now bounded from above by the preceding estimates (cf. lemmas 4 and 7). By the Ascoli-Arzela theorem, each one of the edges can be supposed to converge to a limit edge if we extract an ad hoc subsequence. The limit segmentation is then defined by these limit edges. When passing to the limit, we know by Fatou’s lemma that neither the integral part of the energy nor the length of the edges can increase, and therefore this limit segmentation has an energy smaller than the inf limit of the energies of the sequence. For more (technical) details see [MorS1].

Remark 2 : Elimination of small regions ([Zu, HarS ]). It is easy to deduce from the preceding proof a lower bound on the area of each region of the segmentation. Indeed, take any region O of the segmentation. By Lemma 6, the number of neighbouring regions O_{j} is at least

N(O) ≥ \frac{Cλ}{\sqrt{|O|} osc(g)^{2}} .

Thus the number α of regions of the segmentation is at least

\frac{Cλ}{\sqrt{|O|} osc(g)^{2}} .

By using the upper bound for α, given in lemma 7 we get

\frac{288 |Ω| osc(g)^{4}}{C^{2}λ^{2}} ≥ \frac{Cλ}{\sqrt{|O|} osc(g)^{2}} .

Thus the area of O is bounded from below by a positive constant only depending on g, λ and Ω. Therefore a merging method based on the minimizing of the energy E(K) will spontaneously eliminate the small regions. This process was considered as heuristic and parameter-dependent in [HarS].

Remark 3 : Elimination of thin regions ([Zu, HarS ]). It is also easy to deduce from the above estimates that the regions are not too “thin”, that is, verify an “inverse isoperimetric inequality”. Indeed, each region O verifies for a constant C depending only on g and Ω,

|O|^{\frac{1}{2}} ≥ C · ℓ(∂O) .

This last result can be deduced from remark 2, since the surface of O is bounded from below and ℓ(∂O) from above (as the length of the segmentation is under control). The estimate obtained can be proved with a better constant C (see [MorS1]). Thus the devices based on the elimination of thin regions (see [PapJ]) for instance, and many clustering algorithms [Pav2]) can be considered as implicit in the search of an optimal 2-normal segmentation and do not depend anymore on extra threshold parameters.

Remark 4 : Smoothing of the boundaries ([PavL ]). The 2-normal segmentations have no chance of having boundaries smooth everywhere. However, since their length is under control, a classical geometric measure theory (see for example [Si-L, th.14.3]) asserts that they are almost everywhere C^{1}. Thus the presence of noise, for instance, can alter the regularity of these boundaries but not make them increase indefinitely, as it is the case for some region growing devices. What can be done in order to restore this regularity and how can such a regularization device be deduced from the energy to be minimized ? The answer is in relation (2.2), which asserts that the curvature of an energy minimizing boundary is controlled. This equation shows that the length term of our “simplest energy”, coupled to its bidimensional contrast measuring term, is enough to ensure that the boundaries are analogous to snakes. They tend to keep at places where g has a jump and they remain smooth, since from the bound on the curvature follows a bound on the first derivative [MorS1]. The merging of regions should therefore be completed by a merging of the boundaries in the same sense : imposing therefore that this boundary is smoothed according to the criterion on curvature imposed by the energy. This idea is implicit in [PavL]. As noticed by these last authors, that imposes to treat region boundaries like snakes, but the relation (2.2) proves that still this tool can be deduced from the “simplest energy” that we have considered.

2.4 A pyramidal algorithm constructing 2-normal affine segmentations.

We now consider the problem of defining and computing a 2-normal segmentation. Notice that not all 2-normal segmentations are equally interesting : for instance, the empty segmentation, where Ω is the single region is clearly a 2-normal segmentation. If the scale parameter λ is very large, it is however also a reasonable segmentation since one “pays” a too large energy amount for having any boundary. Now, it is obvious from the definition that the empty segmentation is 2-normal for every λ, which certainly proves that the assertion that a segmentation is 2-normal is not enough to ensure that it is “good”. But if we follow the main idea of the region growing methods [Zu], we shall see that what they compute is precisely a 2-normal subsegmentation of a fine initial segmentation, obtained by recursive merging.

Assume that the datum g is defined on a rectangle. This rectangle is divided in small squares of constant size (the pixels) and g is assumed to be constant on each pixel. Here are the properties which we require for the segmentations computed by a region growing algorithm, defined as an application associating to g and λ a segmentation (u, K).

a) “Correctedness” (Fixed point property) : Assume that g is piecewise constant on some regions of the rectangle. Then there exists a value λ_{0} of the parameter λ such that for every λ < λ_{0}, the segmentation (u, K) obtained by the algorithm verifies u = g and K is the union of the boundaries of the areas where g is constant. This property has been proved to be asymptotically true for the segmentations which are global minima of the energy E as λ tends to zero [Ric]. But we impose it here as a nonasymptotic property.

b) “Causality” (Pyramidal segmentation property) : If λ > λ′, then the boundaries provided by the algorithm for λ are contained in those obtained for λ′ and the regions of the segmentation associated to λ are the unions of some of the regions obtained for λ′.

The last property ensures that a fast pyramidal algorithm can be implemented, computing a hierarchy of segmentations from fine to coarse scales. Moreover the coarser segmentation will be deduced from the finer by “merging” operations, with a pyramidal structure for the computation. Note that, as a consequence of the fixed point property, if λ is very small, the computed segmentation is attained with (u_{0}, K_{0}) where u_{0} = u and K_{0} consists of all the boundaries of all the pixels, and therefore coincides with the global minimum as λ is zero. We shall call this segmentation, where each pixel is a region, the “trivial segmentation”. The recursive merging algorithm which we use verifies all the above mentioned properties.

2.4.1 Description of the algorithm.

The criterion. The decision to proceed to a merging of two regions O_{i} and O_{j} depends on the sign of E(ũ, K \ ∂(O_{i}, O_{j})) − E(u, K). The algorithm looks for a decrease of global energy by merging these regions. This is the criterion of 2-normal segmentations introduced in the discussion of section 2. The simplified Mumford and Shah model is implemented by choosing the following energy functional :

E(u, K) = ∫_{Ω} ∥u − g∥^{2} + λℓ(K).

Here g is a vector valued function, whose components are different channels, defined on the rectangle Ω, u is the approximating vector function, and K is the set of boundaries with total length ℓ(K). As in the piecewise constant case, u is the mean value of g on each region O. Thus the above functional is just denoted by E(K), i.e. to obtain (u, K) one needs to know only K. Then the merging criterion is

E (K \ ∂ (O_{i}, O_{j})) − E(K) = \frac{|O_{i}| · |O_{j}|}{|O_{i}| + |O_{j}|} · ∥u_{i} − u_{j}∥^{2} − λ · ℓ (∂ (O_{i}, O_{j})) ,

where |.| is the area measure and u_{i} the approximation of g on O_{i}. When g is scalar the norm ∥.∥ is just the absolute value. For multichannel data a weighted norm ∥.∥ is used. It is specific to each application and to the meaning of the different channels. This will be emphasized in the next section.

To obtain the necessary data for evaluating the criterion the following information has to be used: Suppose g = ^{t}(g^{1}, …, g^{n}), then to each region O we associate its area |O| and n channels c_{O}^{l} = ∫_{O} g^{l} , (l = 1, …, n). These yield the values for u restricted to O

: u_{O} = ^{t}(u_{O}^{1}, …, u_{O}^{n}) by simply computing u_{O}^{l} = \frac{c_{O}^{l}}{|O|} to get the mean value.

The channels of region O_{new} obtained by merging O_{i} and O_{j} are given by |O_{new}| = |O_{i}| + |O_{j}| and c_{O_{new}}^{l} = c_{O_{i}}^{l} + c_{O_{j}}^{l} , (l = 1, …, n).

Thus a merging of two regions only implies adding the corresponding channels and updating the data structure.

The algorithm.

(i) Take the segmentation (u_{0}, K_{0}) resulting from the initialization and λ_{i} = λ_{1} as a scale parameter.

(ii) Scan the list of regions and for every candidate region look for the adjacent region which yields the best merging score (i.e. the maximal energy decrease). If such a region exists proceed to merge and update the data. The next region in the list becomes a candidate for merging. For λ = λ_{i} fixed, repeat the scanning of the picture until no merging is possible. After this step, a 2-normal subsegmentation (u_{i}, K_{i}) of the initial segmentation for scale parameter λ_{i} is achieved.

(iii) For every λ_{i}, i=1,...,L, calculate a 2-normal segmentation by iterating step (ii). The algorithm stops if there is just one region left or after computing a 2-normal segmentation (u_{L}, K_{L}) for the last scale parameter λ_{L}.

Experimental results. In figure 2.1 we show in the top left corner a satellite image which has been segmented with respect to a single channel, the grey level: in the upper right image we show the piecewise constant reconstruction u, in the lower left picture Chapitre 2 45 the boundary set K is represented and in the lower right corner we have drawn the boundaries upon the original picture.

Figure 2.1 : Scale λ = 2048, 117 regions (original picture : GRECO TdSI). Representation of the original picture g, the piecewise constant approximation u, the boundaries K and the original picture with boundaries.

2.5 Application to texture discrimination.

We follow the ideas of David Marr [Marr] concerning the “raw primal sketch”. According to this theory the grouping process of the human visual system should be based on the detection of local features which from the mathematical point of view are simply differential or semilocal integrodifferential operators like the derivatives or convolution by Gaussians. The preattentive texture discrimination depends on the fact that at least one of the features has a bigger or lower local density ( see [MaliP2]).

It is important to note that most of the channels introduced by the algorithms of [MaliP2], [VoP], [BovCG], … are linear operators which, according to the Wavelet theory, can be computed at each scale in linear time by a pyramidal scheme (also called in this context quadrature mirror filter [Coh], [EsG]).

The wavelet transform allows to compute a multiscale gradient in every direction in linear time (by a quadrature mirror filter). Since, as recalls David Marr, textures “live independently at different scales”, it will be much easier to get discrimination of most textures with a fast multiscale analysis. One can for instance use the Daubechies filters [Daub1, Daub2], which can be computed with quadrature mirror filters having very few coefficients (eight coefficients are enough for a smooth wavelet with two zero moments). In a work in collaboration with Jaffard and Journé, Daubechies has proved a fast algorithm for a variant of Gabor filtering ([DaubJJ]). Fast Gabor filtering may also provide cheap channels for our segmentation algorithm, which prove very useful for almost periodic textures. Note however that the human eye is by no means perfect in texture discrimination and that it can reasonably be argued that a discrimination algorithm could be more accurate and discriminate textures which are not directly accessible to the human eye. In our opinion, a successful computational model of texture discrimination should do better than the human eye in some cases, and therefore be used as an enhancement operator (see figure 2.4).

Let us summarize the assumptions on which we base our segmentation algorithm :

• We assume, following the theory of Beck and Julesz (see [Bec, Ju]) and some recent experimental and computational confirmation of it [BovCG, MaliP2, VoP] that a reasonable number of channels defined on the image is enough in order that for any pair of preattentively different textures, at least one of the channels will help to discriminate them.

• We are not able now to say which channels are necessary or sufficient in order to get a discrimination similar to that of the human eye. Now, our universal segmentation algorithm is designed for fast specific applications, but also as a fast and robust experimental device allowing progress in the question of which channels are necessary in order to match the human visual performance.

Examples.

Let g = (g_{1}, …, g_{n}) be the initial data calculated from the picture g̃ and u = (u_{1}, …, u_{n}) the piecewise constant estimate ( it turns out that u_{i} is the mean value of g_{i} on each connected component). Using the functional

∫_{Ω} ∥u − g∥_{pond.}^{2} + λ · ℓ(K) ,

where ∥.∥_{pond.} is a weighted norm, we calculate a 2-normal segmentation.

The initial datum g for the experiments related below is given in our experiments by an oversampled Haar-Wavelet transform of the picture g̃. More precisely the image is convolved with a bank of linear filters F_{k}, followed by half-wave rectification.

R_{2k} = (g̃ ⋆ F_{k})^{+} (x, y) ; R_{2k+1} = (g̃ ⋆ F_{k})^{−} (x, y)

It should be noticed that all the filters F_{k} are of zero mean and separable:

H_{a}^{1}(x, y) = χ_{[−a,a]}(x) (χ_{[0,a]}(y) − χ_{[−a,o]}(y))
H_{a}^{2}(x, y) = χ_{[−a,a]}(y) (χ_{[0,a]}(x) − χ_{[−a.0]}(x))
H_{a}^{3}(x, y) = (χ_{[0,a]}(x) − χ_{[−a,o]}(x)) (χ_{[0,a]}(y) − χ_{[−a,o]}(y))

where χ is the caracteristic function on IR and a = 2^{j}, 1 ≤ j ≤ J, J is the decomposition order of the analysis. The R_{i}s are then filtered by Gaussians in order to obtain texton densities. The size of these Gaussians corresponds to the “Δ-neighbourhoods” of Julesz. Let us recall the definition of Julesz [Ju]: “ The Δ-neighbourhood is the area, in which differences in texton densities are determined. Textons are formed only if the adjacent elements lie within the Δ-neighborhood.” Finally we obtain the g_{i}s, the components of our initial datum g, as the filtered versions of the R_{i} channels.

The aim of our experiences is to get a verification of the Julesz doctrine : we stop our region growing if the desired number of regions is reached ( e.g. if there are two textures we proceed until we have a partition in two regions of our image) . If the discrimination is successful ( i.e the two regions correspond to the textures’ location) one can say that the used channels are able to discriminate the given textures ( it is important to notice that we don’t use the grey-channel information ).

The experimentations used Brodatz pictures (see [Bro]) and a synthetic image which illustrates the need to keep all the channel involved in the discrimination process.

2.6 Conclusion.

We have proved that the minimizing of the “simplest” segmentation energy entails the implicit realization of properties sought by most segmentation devices. Indeed, Figure 2.2 : Above two Brodatz textures, below the boundary calculated (Δneighbourhood of 16 pixels, level of wavelet transform J = 3) .

the most primitive segmentation tool, the “merging”, applied to the simplest possible segmentation energy, is enough to ensure a compact and therefore small set of possible segmentations, with no small regions and no thin regions. Uniform a priori estimates for the size and number of the regions can be given for all segmentations obtained by exhaustive “merging”. Moreover, the region growing method associated with the recursive merging is enough to retrieve all piecewise constant functions. Such a merging method is not accurate enough to obtain smooth boundaries, but it controls anyway their length. The big advantage of this method is its velocity. \underline{Aknowledgment.} We thank Yves Meyer, Pietro Perona and Lenny Rudin for valuable conversations and remarks.

Figure 2.3 : Left part : two Brodatz textures, right part : three different Brodatz textures. Below each the boundary calculated (same parameters : Δ-neighbourhood of 16 pixels, level of wavelet transform J = 3) . Notice that in this case the same texture is repeated twice left. Unless the boundary between them is visible, the algorithm recognizes the identity of the texture and segments correctly.

Figure 2.4 : Above the textures : on the left hand one dominant channel; on the right hand the same dominant channel with a weak second channel. Below the boundary calculated (same parameters: Δ-neighbourhood of 16 pixels, level of wavelet transform J = 3) . This experiment illustrates the power of the algorithm for discrimination based on nondominant channels. Most texture analysis devices tend to base it on the dominant channel. If the dominant channel is the same in the two regions, no discrimination [BovCG,MaliP2,Vop2] is possible. It is important to notice that a variational method yields an easy formalization for segmentation based on non-necessary dominant channel.

52 Chapitre 2

Références

[Ale] P.S. Aleksandrov. Combinatorial Topology. Vol. 1, chap. 2, Graylock Press NY, 1956.

[Amb] L. Ambrosio. A compactness theorem for a special class of functions of bounded variation. Boll. Un. Mat. Ital. (to appear).

[AmbDeG] E. De Giorgi, and L. Ambrosio. Un nuovo tipo di funzionale del calcolo delle variazioni. To appear in Atti Acad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Natur.

[AmbT] L. Ambrosio, and V. M. Tortorelli. Approximation of functionals depending on jumps by elliptic functionals via G-convergence. Preprint. Scuola Normale Superiore di Pisa.

[AmiTW] A. Amini, S. Tehrani and T. E. Weymouth. Using dynamic programming for minimizing the energy of active contours in the presence of hard constraints. Second International Conference on Computer Vision (Proceedings of), ICCV’88, IEEE n^{∘}883.

[Az1] R. Azencott. Synchronous Boltzmann machines and Gibbs fields; Learning algorithms. To appear, Springer Verlag 1989, (NATO Series, Les Arcs Congress Proceedings).

[Az2] R. Azencott. Boltzmann Machines : high order interactions and synchronous learning. Submitted to IEEE PAMI Trans.

[Bec] J. Beck. Organization and representation in perception. Chapter in Textural se gmentation, Erlbaum, Hillsdale, NJ, 1982.

[BecSI] J. Beck, A. Sutter and R. Ivry. Spatial frequency channels and perceptual grouping in texture segregation. Computer Vision, Graphics and Image Processing, 37, 299-325, 1987.

[BlakZ] A. Blake and A. Zisserman. Visual Reconstruction. MIT Press, 1987.

[BlatM] J. Blat, and JM. Morel. Elliptic problems in image segmentation and their relation to fracture theory. Proceedings of the Int. Conf. on nonlinear elliptic and parabolic problems, Nancy 88 (Longman 89).

[BovCG] A.C. Bovik, M. Clark and W. S. Geisler. Computational texture analysis using localized spatial filtering. IEEE 1987.

[Bro] Ph. Brodatz. Textures for artists and designers. Dover Publications Inc. NY 1966.

[CarDeGL] M. Carriero, E. DeGiorgi and A. Leaci. Existence theorem for a minimum problem with free discontinuity set. Arch. Rat. Mech. Anal., vol. 108, 1989.

[CarLPP] M. Carriero, D. Leaci, E. Pallara and E. Pascali. Euler conditions for a minimum problem with discontinuity surfaces. Preprint Universita degli studi di Lecce (Italy)

[Coh] A. Cohen. Ondelettes, analyses multirésolutions et filtres miroir en quadrature. Ann. Inst. Henri Poincaré, Analyse non linéaire, Vol. 7, n^{∘}5, p.439-459, 1990.

[CohDF] A. Cohen, I. Daubechies and J.C. Feauveau. Biorthogonal bases of compactly supported wavelets. Preprint Bell Labs, 1989.

[CongT] G. Congedo and I. Tamanini. On the existence of solutions to a problem in image segmentation. Preprint, Dipartimento di Matematica, Univ. di Trento, 38050 Italy.

[DalMMS1] G. Dal Maso, J.M. Morel and S. Solimini. Une approche variationnelle en traitement d’images : résultats d’existence et d’approximation. C.R. Acad. Sci. Paris, t.308, Série I, p. 549-554, 1989.

[DalMMS2] G. Dal Maso, J.M. Morel and S. Solimini. A variational method in image segmentation: existence and approximation results. Acta Matematica Vol. 168, 1992.

[Daub1] I. Daubechies. Orthonormal bases of compactly supported wavelets. Comm. Pure and Applied Math., vol. 41, p.909-996, 1988.

[Daub2] I. Daubechies. Orthonormal bases of compactly supported wavelets. Part II, II, Preprint Bell Labs, 1989.

[DaubJJ] I. Daubechies, S. Jaffard and J.L. Journé. A simple Wilson orthonormal basis with exponential decay. Preprint Bell Labs, 1989.

[DeG] E. De Giorgi. Free discontinuity problems in calculus of variations. To appear, Proc. of the meeting in J.L. Lions’s honour, Paris, 1988.

[DeGCL] E. DeGiorgi, M. Carriero and A. Leaci. Existence theorem for a minimum problem with free discontinuity set. Arch.Rat.Mech.Anal., Vol.108, 1988, p.195-218.

[DibK] F. Dibos and G. Koepfler. Propriété de régularité des contours d’une image segmentée. C.R. Acad. Sci. Paris, t. 313, Série I, p. 573-578, 1991.

[EsG] D. Esteban and C. Galand. Application of QMF to split-band voice coding schemes. Proc. of Int. Conf. Acoust. Speech, Signal Processing, Harford, 1977.

[Fed] H. Federer. Geometric Measure Theory. Springer-Verlag,1969.

[FoT] I. Fonseca and L. Tartar. The gradient theory of phase transitions for systems with two potential wells. Carnegie Mellon Univ. To appear.

[FroM] J. Froment and J.M. Morel. Analyse multiéchelle, vision stéréo et ondelettes. Exposé n^{∘}5 du Séminaire d’Orsay “Les ondelettes en 89”, ed. P.G. Lemarié, Springer Verlag 1990.

[GemG] S. Geman and D. Geman. Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE PAMI 6, 1984.

[HarS] R.M. Haralick and L.G. Shapiro. Image segmentation techniques. Computer Vision Graphics and Image Processing, 29, 100-132, 1985.

[Ju] B. Julesz. Texton gradients : the texton theory revisited. Biological Cybernetics, 54 : 245-251, 1986.

[KasWT] M. Kass, A. Witkin and D. Terzopoulos. Snakes: active contour models. 1st Int.Comp.Vis.Conf. 1987, IEEE n^{∘}777.

[KoeMS] G. Koepfler, J.M. Morel and S. Solimini. Segmentation by minimizing a functional and the “merging” methods. Proceedings ‘GRETSI Colloque’ in Juan-les-Pins, September 1991 (France).

[Koepf] G. Koepfler. Formalisation et analyse numérique de la segmentation d’images. Rapport de thèse. 1991, université Paris Dauphine.

[MaliP1] J. Malik and P. Perona A scale space and edge detection using anisotropic diffusion. Proc. IEEE Computer Soc. Workshop on Computer Vision, 1987.

[MaliP2] J. Malik and P. Perona. Preattentive texture discrimination with early vision mechanisms. Journ. of the Opt. Society. of America. A, 923-932, vol. 7, n^{∘}5, 1991.

[Mall] S.Mallat. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans. on Pattern Analysis and Machine Intelligence, 674-693, July 1989.

[Marr] D. Marr. Vision. Freeman and Co. 1982.

[Me] Y. Meyer. Ondelettes, fonctions splines et analyses graduées. Rapport CERE-MADE n^{∘}8703, 1987.

[MorS1] J.M. Morel and S. Solimini. Segmentation of images by variational methods: a constructive approach. Rev. Matematica de la Universidad Complutense de Madrid, Vol.1, n^{∘}1,2,3, 169-182, 1988.

[MorS2] JM. Morel and S. Solimini. Segmentation d’images par méthode variationnelle : une preuve constructive d’existence. C. R. Acad. Sci. Paris Série I, 1988.

[MumS1] D. Mumford and J. Shah. Optimal Approximations by Piecewise Smooth Functions and Associated Variational Problems. Communications on Pure and Applied Mathematics. vol. XLII No.4, 1989.

[MumS2] D. Mumford and J. Shah. Boundary detection by minimizing functionals. Image Understanding, 1988, Ed. S. Ullman and W. Richards.

[PapJ] N.T. Pappas and N.S. Jayant. An adaptive clustering algorithm for image segmentation. Proc. of ICCV’88, IEEE n^{∘} 883, 310-315, 1988.

[Pav1] T. Pavlidis. Segmentation of pictures and maps through functional approximation. Comp.Gr. and Im.Proc. 1, 360-372, 1972.

[Pav2] T. Pavlidis. Structural Pattern Recognition. Springer, New York 1977.

[PavL] T. Pavlidis and Y.T. Liow. Integrating Region Growing and Edge Detection. Proc. of the IEEE Conf. Comp. Vision and Patt. Recognition 1988.

[Ric] T.J. Richardson. Scale independant piecewise smooth segmentation of images via variational methods. PhD. Dissertation, Lab. for Information and Decision Systems, MIT, Cambridge MA 02139, February 1990.

[Si-L] L. Simon. Lectures on Geometric Measure Theory. Centre for Math. Analysis. Australian Nat. Univ., vol. 33, 1983.

[Th] C. Thomassen. The Jordan-Schönflies theorem and the classification of surfaces. Am.Math.Monthly, 116-130, February 1992.

[VoP] H. Voorhees and T. Poggio. Computing texture boundaries from images. Nature, 333, 364-367, 1988.

[Zu] S.W. Zucker. Region growing: Childhood and Adolescence (Survey). Comp. Graphics and Image Proc. 5, 382-399, 1976.

Axiomatization of shape analysis and application to texture hyperdiscrimination.

Ce chapitre a été écrit en collaboration avec Jean Michel Morel.

Abstract. In this paper ^{1}, we first prove that, under four simple axioms, the multiscale analysis of shapes is given by a curvature motion equation. In a second part, we illustrate the advantages of such an axiomatic analysis to discuss a psychophysical theory of early vision due to B. Julesz : the texture preattentive discrimination theory. The result was unexpected: We prove that the Julesz axiomatic is too good for human vision and leads to a hyperdiscrimination algorithm.

3.1 Introduction.

Shape analysis methods try to extract from digital shapes local geometric ”features” like holes, corners, orientation, inflexion points ... on which a further ”intelligent” analysis could be based. In order to extract such features, the local analysis implies intuitively that one should use differential operators. However, because of perturbations due to the nature of the picture as well as to the recording conditions, pure differential operators are completely unreliable. Therefore, the vision research has defined as a paradigm the concept of multiscale analysis. Multiscale analysis aims at a smoothing of the shape. This smoothing allows the reliable extraction of local and semilocal features, but introduces a scale parameter, t, which measures the loss of locality (and the win in globality) of the extracted features, see [Koe], [Wit], [MaMo], ... We shall prove that, under a simple list of obvious axioms, the multiscale shape analysis is given by a curvature motion equation. In a second part of this conference, we shall show numerical experiments illustrating the advantages of such an analysis to discuss a theory of early vision due to B. Julesz : the texture preattentive discrimination theory. We shall see that this theory is formalized by a curvature motion and can be discussed in this framework. In particular, we shall see that the Julesz texton theory is ”too good” as a model for the human perception. Indeed, we prove that a formal multiscale analysis based on this model must be a curvature motion, and achieves ”hyperdiscrimination” of textures which should be undiscriminable for Julesz theory.

1Conference held at Trento, July 92. To appear in proceedings on ”Surface Tension and Movement by Mean Curvature”. De Gruyter Publishers (Berlin).

In section 2, we state the four principles of shape analysis. In section 3, we prove that they lead to a single special curvature motion. Section 4 is devoted to the generalization of the model for image processing. In section 5, we present the Julesz axiomatic of texture discrimination. Section 6 gives a translation of Julesz axioms into a mathematical theory and section 7 gives numerical results on Julesz textures.

3.2 The four principles of Multiscale Shape Analysis.

Define a shape or silhouette as a closed set X whose boundary is a Jordan curve of ℜ^{2}. We denote by T_{t}(X) the shape analysed at scale t. X is identified with its characteristic function X(x) = 1 if x ∈ X and 0 else. We call multiscale analysis any family operators (T_{t})_{t≥0} acting on shapes and we set X(t) = T_{t}(X).

We first consider a principle which is called in vision theory the pyramidal architecture property. According to this assumption, T_{t} can be computed from T_{s} for any s ≤ t and T_{0} is of course the identity. This is natural, since a coarser analysis of the original picture is likely to be deduced from a finer one without any dependence upon the original picture. Of course the finest picture analysis is the identity. In order to formalize the pyramidal architecture we include T_{t} = T_{t,0} into a family of transition operators T_{t,s} indexed by 0 ≤ s ≤ t and satisfying [Pyramidal architecture] T_{t+s} = T_{t+s,s}T_{s} for all 0 ≤ s, 0 ≤ t

The shape analysis process is thought of as a flow of the shape through different filters, each associated with a scale t. If we discretize such a process, it yields a serie of operators S_{n} = T_{(n+1)h,nh} applied successively.

Assume that X and Y are two silhouettes and that for some x ∈ ∂Y and some r > 0, one has X ∩ B(x, 2r) ⊂ Y ∩ B(x, 2r). Assume further that the inclusion is strict in the sense that ∂X and ∂Y only meet eventually at x. Then we shall say that the shape X is included in shape Y around x.

[Shape local inclusion] If X is included in Y around x, then for h small enough, T_{t+h,t}(X) ∩ B(x, r) ⊂ T_{t+h,t}(Y) ∩ B(x, r).

This last axiom implies that the value of T_{t+h,t}(X) for h small, at any point x, is determined by the behaviour of X near x. We are allowed to take r infinite. Therefore the shape local inclusion principle also implies that if a shape is globally contained in another, this order is preserved for every scale [MaMo]. Notice that the combination of both pyramidal and shape local inclusion principles is generally called in Image Processing causality [Koe], [Wit], [Kim], [SaTa], [KTZ]. Indeed, the causality assumption in Image Processing states that when the scale increases, no new structure appears in the image.

Both preceding principles allow, as we shall see, to localize the shape analysis in space and time and we therefore only need to state what the multiscale analysis makes of very simple shapes in order to specify it. So we add a basic ”principle” which will state what happens to disks. As we shall see, disks are somehow the ”basis” on which every shape can be decomposed.

[Basic principle] Let D = D(x, 1/r) be a disk with curvature 1/r and center x. Then T_{t+h,t}(D) is a disk with radius ρ(t, h, 1/r) and center x. Moreover, the function h → ρ(t, h, 1/r) is differentiable with respect to h at h = 0 and the differential is continuous with respect to 1/r.

The ”basic principle” implies that the multiscale analysis behaves in a smooth and isotropic way.

In the following, we set

g(t, 1/r) = \frac{∂ρ}{∂h}(t, 0, 1/r)

Note that g(t, s) is defined for t ≥ 0 and s ∈ ℜ, now to define g(t, 0) we must assume that lim g(t, 1/r) exists when r tends to +∞ or −∞ . The radius r may be positive or negative, according to the orientation of the normal n⃗(x). In the case where the curve is a Jordan curve enclosing a set X, we take as a convention that n⃗(x) is pointing outside X and the curvature is negative if X is convex at x, positive else. It may seem natural to assume therefore that g(t, κ) is odd with respect to κ. Indeed, this corresponds to the assumption that a black disk on grey background and a white disk behave in a same way.

Lemma 1.The function g(t, κ) = \frac{∂ρ}{∂h}(t, 0, κ) is nondecreasing with respect to κ.

Consider two disks D and D′ with respective radii |r| and |r′| > |r| and keep in mind the above convention about the orientation. So we have algebraically: r′ < r and κ < κ′. Let D be included in D′ and tangent to D′ at some point x. By the shape local inclusion principle T_{t+h,t}(D) is contained in T_{t+h,t}(D′). By the basic principle both sets are disks concentric to D and D′ respectively and with radii ρ(t, h, κ) and ρ(t, h, κ′). The inclusion of both disks therefore yields

ρ(t, h, κ) − ρ(t, 0, κ) ≤ ρ(t, h, κ′) − ρ(t, 0, κ′)

Dividing by h and passing to the limit when h tends to 0 yields the result.

The next optional shape preserving principle is the affine invariance of shape analysis. Set AX(x) = X(Ax) for any linear map A. The affine invariance can be stated as:

[Affine inv.] For any A and t ≥ 0, there exists a C^{1} function t′(t, A) ≥ 0 such that AT_{t′(t,A),t′(s,A)} = T_{t,s}A. Moreover, the function ϕ(t) = \frac{∂t′}{∂λ} (t, Id) is positive for t > 0.

This relation means that the result of the multiscale analysis T_{t} is independent of the size and position in space of the analyzed features : this is very important in the world where we live, since the same plane object X can be seen at very different distances and orientations and therefore with any affine change X(x) → X(Ax). A linear map corresponds to the anamorphosis of a plane image when it is presented to the eye at an arbitrary distance and with an arbitrary orientation in space.

Lemma 2. [AGLM](Normalization of scale) Assume that t → T_{t} is a one to one family of operators verifying [affine invariance]. Then the funtion t′(t, B) only depends on t and |detB|: t′(t, B) = t′(t, |detB|^{1/2}) and is increasing with respect to t. Moreover, there exists an increasing differentiable rescaling function h : [0, ∞] → [0, ∞], such that t′(t, B) = h^{−1}(h(t)|detB|^{1/2}) and if we set S_{t} = T_{h^{−1}(t)} we have t′(t, B) = t|detB|^{1/2} for the rescaled analysis.

Remarks.

a) Thus for any multiscale analysis, we can normalize the scale-space relation and give a spatial interpretation to t under a simple rescaling.

b) The assumption on t′, ϕ(t) = \frac{∂t′}{∂λ} (t, Id) > 0, can be interpreted by looking at the relation AT_{t′} = T_{t}A when the scale λ = |detB|^{1/2} increases, i.e when the image is shrinked before analysis by T_{t}. Then, the corresponding analysis time before shrinking is increased. In more informal terms we can say that the analysis time increases with the size of the picture.

The proof of the lemma 2 is given in [AGLM] and in appendix.

3.3 The fundamental equation of Shape Analysis.

When a point x belongs to an evolving curve, we denote by ẋ the time derivative of x, which is a vector of ℜ^{2}. By curv(x) we denote the curvature of a curve which is C^{2} at x. Recall that the curvature is defined as the inverse of the radius of the osculatory circle to the curve at x. The curvature is zero if the radius is infinite. By n⃗(x) we denote the vector normal to the curve. In the case where the curve is a Jordan curve enclosing a set X, we take as a convention that n⃗(x) is pointing outside X and the curvature is negative if X is convex at x, positive else.

Theorem 1. (i) Under the three principles (pyramidal, local shape inclusion, ”basic”), the multiscale analysis of shapes is governed by the curvature motion equation

ẋ = g(t, curv(x))n⃗(x)  (1)

where g is defined in lemma 1. (ii) If we assume [Affine inv.], then the equation of the multiscale analysis is, up to rescaling:

ẋ = γ(t.curv(x))n⃗(x)  (2)

where γ is defined as following: γ(x) = a.x^{1/3} if x ≥ 0 and γ(x) = b.x^{1/3} if x ≤ 0 and a, b are two nonnegative values. (iii) If we add that T_{t}(X^{c}) = T_{t}(X)^{c} [Reverse contrast invariance] then the function g in (i) is odd and a = b in (ii).

By the expression ”governed by”, we mean that the equation must be verified at any point (t, x) of X(t) where the silhouette is smooth enough to give a classical sense to both terms of the equation.

Remark: It has been proved [GaHa], [Gra], [Ang] that the equation (1) with g(t, s) = s has smooth solutions. This equation has been introduced in picture processing by [KTZ], [MaMo] and [AGLM] in different contexts. An early version of an algorithm leading to equation (1) in [KoDo].

The equation (2) has been introduced and axiomatically justified (with a more complicated axiomatic however) in [AGLM]. It has also been proposed in [SaTa]. Existence and regularity of the solution are proved in [SaTa].

(i) Let X be a silhouette and assume that T_{t}(X) has a boundary which is a C^{2} manifold in a neighbourhood of a point x of ∂X. Then it has a curvature κ at point x and we consider a subosculatory and a surosculatory disk, that is, a disk D with curvature κ − ϵ and a disk D′ with curvature κ + ϵ, both tangent to the silhouette at x. Applying the same two principles as in the lemma, we see that

T_{t+h,t}(D) ∩ B(x, r) ⊂ T_{t+h,t}(X) ∩ B(x, r) ⊂ T_{t+h,t}(D′) ∩ B(x, r)

Thus, denoting by x(t + h) the point of ∂T_{t+h}(X) such that x(t + h) − x(t) is parallel to n⃗(x), we obtain

ρ(t, h, κ − ϵ) − ρ(t, 0, κ − ϵ) ≤ (x(t + h) − x(t)).n⃗(x) ≤ ρ(t, h, κ + ϵ) − ρ(t, 0, κ + ϵ)

Dividing by h and passing to the limit when h tends to 0 yields

\frac{∂ρ}{∂h}(t, 0, κ − ϵ) ≤ lim inf \frac{x(t + h) − x(t)}{h}.n⃗(x)

lim sup \frac{x(t + h) − x(t)}{h}.n⃗(x) ≤ \frac{∂ρ}{∂h}(t, 0, κ + ϵ)

We obtain equation (1) by passing to the limit when ϵ tends to 0 and using the fact that κ → \frac{∂ρ}{∂h}(t, 0, κ) is continuous.

(ii) After renormalization, we can use the identity

T_{t+h,t}D_{λ} = D_{λ}T_{(t+h)λ,tλ}

(where D_{λ} = λId) and so, we can deduce that the function ρ of the basic principle must verifye ρ(t, h, \frac{λ}{r}) = \frac{1}{λ}ρ(λt, λh, 1/r) Therefore, we obtain the relation (after differentiation with respect to h at 0)

g(t, λs) = g(λt, s)

for any t > 0, λ > 0 and s ∈ ℜ Changing t in t/λ and taking λ = 1/t we get g(t, s) = g(1, ts) = β(ts) for the function β defined as β(x) = g(1, x).

(iii) After renormalization, we can use the identity T_{t+h,t}A = AT_{t+h,t}, where A is the linear transform whose determinant is one

(x, y) → (λx, (1/λ)y) λ > 0

Let us apply this identity to the unit disk Δ. Look at the point x_{0} = (1, 0) on the boundary of Δ. Then the velocity of x_{0} is β(−t), and this velocity is transformed into λβ(−t). Now, look at AΔ. Since AΔ is an ellipse with curvature −λ^{3} at point Ax_{0}, the velocity of Ax_{0} is β(−t.λ^{3}). Using the first identity, we obtain β(−t.λ^{3}) = λβ(−t). Taking t = 1, we get β(x) = b.x^{1/3} for x < 0 (b = β(−1)). Now, apply the same technique to Δ^{c} and we get the result β(x) = a.x^{1/3} for x > 0 (a = β(1)).

(iii) With the same technique as above we obtain that the function β is odd.

3.4 The fundamental equation of Image Processing.

Since the multiscale analysis satisfies the obvious inclusion principle that if A ⊂ B then T_{t+h,t}(A) ⊂ T_{t+h,t}(B), then we can, as well known in the ”mathematical morphology school” [Ser], [Mara],

associate with a picture u the set of its level sets

X_{a}u = {(x, y), u(x, y) ≥ a}

Then, assuming that each level set is a union of silhouettes, we can simply define T_{t}(u) as

[morphological principle] For any a, t, h and u: X_{a}T_{t+h,t}(u) = T_{t+h,t}(X_{a}u)

The corresponding multiscale analysis of images is governed by

\frac{∂u}{∂t} = g(t, curv(u))|Du|  (3)

We say that T_{t} is governed by equation (3) if for any u C^{2} and every value a such that u(x) = a implies Du(x) ≠ 0 , the equation is true. However, a generalized sense can be given by using the concept of ”viscosity solution” [AGLM]. Conversely, one can find in [CGG] a proof of the existence of viscosity solutions for this equation. We shall only prove here that the equation (3) is verified when the image u is smooth enough.

Let us prove (3) at any (x, a) such that u(t, x) = a and (u(t, y) = a implies Du(t, y) ≠ 0). The last condition implies by the implicit function theorem that x belongs to a Jordan curve Γ enclosing two regions Γ^{+} and Γ^{−}, with u(t, x) < a on Γ^{−} , u(t, x) > a on Γ^{+} and u(t, x) = a on Γ. Set u(t) = u(t, .). The relation X_{a}T_{t+h,t}u(t) = T_{t+h,t}X_{a}u(t) implies that T_{t+h,t}Γ^{+} is equal to a connected component of X_{a}u(t + h). Therefore, the point x(t+h) defined by equation (1) with initial value x(t) = x belongs to the boundary of X_{a}u(t+h). So we obtain u(t+h, x(t+h))−u(t, x(t)) = 0. Dividing by h and passing to the limit yields

\frac{∂u}{∂t}(t, x) + Du(t, x)ẋ(t) = 0

and using equation (1) and n⃗(x) = −\frac{Du}{|Du|} we obtain equation (3).

3.5 The ”Texton Density” Doctrine of Julesz.

According to Julesz theory [Jul1], texture discrimination is based on what he calls ”local features” (See figure 4 for an example of a pair of discriminable textures and figure 5 for an example of undiscriminable one). He adds that ”those local features that might yield texture discrimination should be very special ones. They should be regarded as the basic elements of pre-attentive perception, and I have named such hypothetical elements ”textons”. (...) It seems that elongated blobs of different width, length, and orientation, and perhaps of different first order statistics in small, local neighborhoods are the local conspicuous features yielding texture discrimination. (...) Besides colour, there are two additional texton classes: elongated blobs (of given orientation, width and aspect ratio) and their terminators. Pre-attentive texture discrimination is based either on the difference in the textons or the difference in the first-order statistics (or some parameter defined by the first-order statistic, such as density or standard deviation) of the textons.” (...) One can raise (...) the possibility that the pre-attentive textural system (...) can evaluate only first-order statistics of textons. (...) Whatever the final number of these texton classes will be, it is evident that for texture discrimination, the Gestaltist’s view has not been corroborated : the global percept of texture was found to be decomposable into more basic elements.”

According to Julesz [Jul2], the axioms of the texton theory are as follows :

Axiom 1: Human vision operates in two distinct modes : Preattentive vision. A system that is parallel (independent of number of elements), acting almost instantaneously, without scrutiny, and covering a large visual field, as in texture discrimination. Basically a preconscious system that knows where texton gradients occur, but is unaware of what these textons are. Attentive vision. Serial search (scrutiny) by a small aperture of focal attention in steps lasting 15-50 ms each (much faster than eye-movement scanning), as in form perception. Only in this attentional aperture is it consciously known what the textons are.

Axiom 2: Definition of textons (Restricted to line drawings) -Nonoverlapping line segments (black on white background, or white on black background, with specific length, orientation, width, velocity, binocular disparity, and flicker rate). -Terminators (ends-of-lines) of line segments. -Crossings of line segments. -Blobs, whose shape is unimportant, except for average length, width, and orientation.

Axiom 3: Texture formation, texton density, texton gradients -The area called the Δ-neighborhood, in which differences in texton densities (numbers) are determined and texton-gradients computed, is rather limited, and proportional to the size of the texture elements. Thus, textures are formed only if the textons are not to sparse, that is, if the adjacent textons lie within the Δ-neighborhood. -The Δ-neighborhood and the aperture of focal attention are scaled by the average extent of the texture elements. So the only global operation postulated by the texton theory is the determination of the average element size needed to scale the range Δ of texton interactions.

Axiom 4: The role of textons -Preattentive vision directs attentive vision to the location where texton-gradients (i.e. differences in textons, or in the densities (number of textons)) occur, but ignores the positional relationships between adjacent textons (i.e. the coupling (”glue”) between textons is missing).

To summarize, preattentive vision first determines the scale from an aggregate of similar elements, i.e. the average element size, and then computes the texton gradients at the boundary of textures, according to the axioms given above. Thus, in preattentive vision, precise positional relationships between ”textons” are not important; only densities matter. The figure proposed by Julesz [Jul1] illustrates this fact by showing two undiscriminable textures (fig. 5) wich, for Julesz, have the same density of textons but, however, made of different patterns.

3.6 How to define Texton Densities. An Axiomatic Approach.

The texton doctrine, which is based only on synthetic images, neglects to consider an important issue : how can such textons be reliably extracted from images of natural scenes ? Julesz [Jul2] left the question open by accepting a priori any plausible system proposed by biologists or vision researchers : ”Are line segments just thin, elongated blobs of white and black (or vice versa), or do they have the luminance profile of the ”Mexican-hat-function” or a Gabor-function? These questions can be decided only by doing the psychophysics.”

We can give a rigorous answer to this question by noting that the texton characteristics indicated by Julesz are shape elements based on curvature and orientation. Furthermore, these elements are multiscale because of the a priori unknown size of the textons. Therefore texton densities are nothing but multiscale curvature and orientation densities. Now, we have proved in theorem 1 that there is only one way to compute multiscale curvature and orientation. In order to completely formalize the Julesz theory, we can consider the following equivalences between formal and intuitive concepts:

”preattentive” ⟺ [pyramidal]

(because no time is given to a feedback in the perceptual system, which therefore acts as a serie of filters applied successively.)

”line drawings” ⟺ [morphological]

Indeed, the fact that Julesz admits that the performance of human vision can be characterized by its behavior on line drawings is equivalent to the assumption that the structure of a picture is contained in its level sets. Notice that the Julesz line drawings have some fitness and their boundaries are ”silhouettes” or finite union of silhouettes as discussed above.

”local features” ⟺ [Basic principle], [local inclusion of shape]

(because the analysis of shape is local for the filters mentionned above. Furthermore, the inclusion preservation is due to the fact that the eye is able to give the right order of sharpness (i.e. curvature) between shapes observed locally).

”Black on white or white on black background” ⟺ [Reverse contrast inv.]

”only texton densities matter” ⟺ local mean values of multiscale curvature and orientation are computed.

We add the principle [Affine invariance] because it is well known that the discriminability of plane textures is little affected by distance and change of position of the plane observer.

Notice that this formalization of Julesz theory simplifies a lot the kind of texton to be considered: we reduced the blobs, terminators,... to multiscale curvature and orientation. Now, this reduction is obtained from theorem 1 since no other local feature can be computed multiscale. We shall now see that, in any case, terminators and elongated blobs are characterized by their multiscale evolution.

3.7 Experimental verification of Julesz Theory.

The preceding equivalences lead to define directly texton densities as local means of multiscale curvatures and orientations. Thus, we do not need to really decompose a picture into ”blobs” of several sizes and orientations [MaPe]. However,we almost do it in the following sense: Grayson and Gage-Hamilton proved that any silhouette becomes convex before disappearing when its evolution is governed by equation (1) with g(t, s) = s. In the affine invariant case g(t, s) = s^{1/3}, a similar result has been announced by Sapiro and Tannenbaum: any silhouette becomes convex and behaves asymptotically like an ellipse. Thus the ”blobs” characterizing any shape in the theory of early vision can be quite well associated with those asymptotic ellipses. This gives a strong confirmation of intuitions contained in the book of D. Marr [Marr]. According to D. Marr, early vision obeys to what he calls the raw primal sketch and this sketch is essentially a decomposition of the image into a family of elliptical blobs of several shapes and orientations.

Formalization and computation of textons densities Let u_{0}(x) be an image and u(t, x) its multiscale analysis by equation (3). Then the curvature and orientation at scale t and location x are given respectively by curv(u)(x) and \frac{Du}{|Du|}(x) which are both well defined and computed by equation (3).

Notice that ”blobs” and textons can be white on black or white on black, so we define for any θ ∈ [0, π] and t > 0 the ”density on a ball B(x, Δ) of textons at scale t and orientation θ” by:

G_{B_{Δ}}(x) ∗ (1_{Θ_{θ,δθ}}curv^{+}(u)(t, x))
G_{B_{Δ}}(x) ∗ (1_{Θ_{θ,δθ}}curv^{−}(u)(t, x))

where Θ_{θ,δθ} is the set of points where \frac{Du}{|Du|} = (cosϕ, sinϕ) verifies |θ − ϕ| ≤ δθ. and G_{B_{Δ}}(x) is the isotropic gaussian with variance Δ. By a^{+} we mean max(a, 0) and by a^{−} max(−a, 0).

Then we get a formal definition of texton densities at each scale and orientation. Based on the obtained texton density channels, a segmentation of the image can be achieved. If the Julesz theory is correct, all the pairs of textures which we present as discriminable should be clearly discriminated by at least one of the texton densities.

Conversely, any pair of undiscriminable textures should be undiscriminated by all channels.

The basic facts that we discovered by experimentation are following. The textures proposed by Julesz, and most textures proposed in the litterature as discriminable are easily discriminated by texton densities even rougher than the announced above: Indeed, it is enough to compute global curvature densities, that is, we take δθ = π and do not take into account orientation.

Let us now show numerical experiments. First, we analyse basic patterns of Julesz and present their multiscale analysis obtained with the last equation of theorem 1. We store the results until scale 9 in figure 1. From this analysis we extract the multiscale positive curvature (figure 2) and the negative one (figure 3) at the same scales.

Second, we analyse in the same way two pairs of textures. In figure 4, we store the 256X256 original image, the analysed image at scale 4, the positive and negative curvature at scale 4 and the result of a segmentation algorithm applied on the 18 scalecurvature images (the 9 positive curvatures and the 9 negative). We also show the negative curvature at scale 1 to emphasize the need of the scale. In figure 5, we store the original 512X512 image, the corresponding images at scale 7 and we present the result of the segmentation algorithm. It is very surprising because in that case, we achieve hyperdiscrimination.

The discrimination is easily achieved and one can even compute the boundary by using the segmentation algorithm commented in [KLM].

N.B. The length unit is the pixel size. In order to make the scale understable in experiments, we normalized the curvature equation in such a way that T_{t} lets disappear a disk with radius r in time r.

3.8 Conclusion.

The axiomatisation of image processing allows an identification of a single invariant model and subsequent discussion of all psychophysical theories of vision. Moreover, the verification of these theories can face a priori two types of answers: a confirmation or a negation. Now, we have seen a third answer: rationalized vision can be better than the human eye! Psychophysics tends to idealize, as well as mechanics, the behavior of real beings. So the ”formalization” of vision goes beyond the experimental expectations and opens the possibility of hyperacuity properties.

A. Proof of the renormalization lemma.

(i) For linear transforms B and C and any t one has the semigroup property t′(t, BC) = t′(t′(t, B), C). Indeed, we have BCT_{t′(t,BC)} = T_{t}BC = BT_{t′(t,B)}C = BCT_{t′(t′(t,B),C)}. The map which associates T_{t} with t being one to one, this implies the stated relation.

(ii) t′(t, B) is increasing with respect to t.

Let us prove that t′(t, A) is one to one with respect to t for any A. Indeed, if not, there would be some A and some (s, t) such that t′(t, A) = t′(s, A). Thus T_{t}A = AT_{t′(t,A)} = AT_{t′(s,A)} = T_{s}A and therefore t = s because T_{t} is one to one. Since t′(t, A) is non negative (by definition), one to one, continuous with respect to t and t′(0, A) = 0 we can deduce that it is increasing with respect to t.

Moreover if B is orthogonal then t′(t, B) = t.

Let R be orthogonal. Then iterating the formula of (i) we have:

t′(t′(t′(t′(...t′(t, R)..., R), R), R) = t′(t, R^{n})

Remark that there is a subsequence of R^{n} tending to Id. Indeed, there is a subsequence R^{n_{k}} which converges to H orthogonal because the orthogonal group is compact. Therefore, the subsequence R^{n_{k+1}−n_{k}} converges to Id . Since there exists a subsequence of R^{n} tending to Id and since t′ is continuous we have for this subsequence lim t′(t, R^{n}) = t′(t, Id) = t. Assume by contradiction that t′(t, R) = t” with t” < t then t′(t′(t, R), R) = t′(t”, R) ≤ t′(t, R) = t” and by recursion, t′(t, R^{n}) = t′(t′(t′(t′(...t′(t, R)..., R), R), R) ≤ t” < t. This is a contradiction. Thus t′(t, R) ≥ t. We prove the converse inequality in the same way and we obtain t′(t, R) = t.

(iii) For any orthogonal R, we have proved that t′(t, R) = t. Now, we note that any linear transform B of ℜ^{2} can be obtained as a product of orthogonal transforms and of linear transforms of the kind A(λ): (x, y) → (λx, y) where λ is non negative. We only need to make a singular value decomposition of B: B = R_{1}DR_{2}, where R_{1} and R_{2} are both orthogonal transforms and D is a transform of the kind (x, y) → (λ_{1}x, λ_{2}y) where λ_{i} are non negative. Now, it is clear that D can be decomposed as D = A(λ_{1})RA(λ_{2})R^{−1} where R is the orthogonal transform: (x, y) → (−y, x). Using (i), t′(t, R_{i}) = t for any orthogonal transform R , the singular value decomposition and A(λ_{1})A(λ_{2}) = A(λ_{1}λ_{2}), we obtain

t′(t, B) = t′(t, λ_{1}λ_{2}) = t′(t, |detB|^{1/2})

(iv) Using (i) and (ii), we have t′(t, λμ) = t′(t′(t, μ)λ) for any positive λ and μ. Differentiating this relation with respect to μ at μ = 1 yields

(0) λ\frac{∂t′}{∂λ}(t, λ) = \frac{∂t′}{∂λ}(t, 1)\frac{∂t′}{∂t} (t, λ)

Choose h such that ϕh′ = h and set t′(t, λ) = G(t, h(t)λ), where

h(t) = exp(∫_{1}^{t} ds/ϕ(s))

Then the preceding relation (0) yields \frac{∂G}{∂x}(x, y) = 0. Thus G(x, y) = β(y) for some differentiable nondecreasing function β. We obtain that t′(t, λ) = β(h(t)λ). Returning to the definition of ϕ(t), we have ϕ(t) = \frac{∂t′}{∂λ} (t, 1) = \frac{∂β(h(t)λ)}{∂λ}(t, 1) and ϕ(t) = h(t)β′(λh(t)) = ϕ(t)h′(t)β′(h(t)λ) thus the derivative of β(h(t)) is 1 and integrating this last relation between 0 and t yields β(h(t)) = t + β(h(0)). Using the fact that t′(0, λ) = 0 (which derives from injectivity of the T_{t}), we obtain β(h(0)) = 0 and therefore t′(t, λ) = h^{−1}(λh(t)). To finish the proof, we set S_{t} = T_{h^{−1}(t)} and we prove that the affine invariance is true for S_{t} with t′(t, λ) = λt. S_{t}B = T_{h^{−1}(t)}B = BT_{t′(h^{−1}(t),λ)} = BT_{h^{−1}(λh(h^{−1}(t)))} = BT_{h^{−1}(λt)} = BS_{λt}.

Chapitre 3 71

fig1. Affine Invariant Multiscale Analysis. The original 240x240 image is analysed from scale one untill scale 9.

72 Chapitre 3 fig2. Multiscale positive curvature.

Chapitre 3 73 fig3. Multiscale negative curvature.

74 Chapitre 3

fig4. Multichannel segmentation using multiscale positive and negative curvatures (18 channels).

Chapitre 3 75 fig5. Multichannel segmentation.

76 Chapitre 3 fig6. Multichannel segmentation on four pairs of textures.

Chapitre 3 77 fig7. Multichannel segmentation.

78 Chapitre 3 fig8. Multichannel segmentation.

References :

[AGLM1] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomes et equations fondamentales du traitement d’images. Report 9216, CEREMADE. Université Paris Dauphine, 1992.

[AGLM2] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomatisation et nouveaux operateurs de la morphologie mathématique. C.R. Acad. Sci. Paris, t.315, Serie I, p. 265-268, 1992.

[AhSc] N. Ahuja and B. Schachter. Pattern Models, Wiley, New York 1983.

[ALM] L. Alvarez, P.L. Lions and J.M. Morel. Image selective smoothing and edge detection by nonlinear diffusion (II). SIAM Journal of numerical analysis 29, pp 845-866. 1992.

[Ang] S. Angenent. Parabolic equations for curves on surfaces I, II. University of Wisconsin-Madison Technical Summary Reports, 89-19, 89-24, 1989.

[AsBr] H. Asada and M.Brady. The curvature primal sketch. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 8, No. 1, January. 1986.

[Bar] G. Barles. Remarks on a flame propagation model. Technical Report No 464 INRIA Rapports de Recherche, December 1985.

[BaSo] G. Barles and P.E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equation, Asymp. Anal., to appear.

[BeJu] J. Bergen and B. Julesz. Rapid discrimination of visual patterns. IEEE Trans. on Systems, Man, and Cybernetics, 13 (5), 1983.

[BPR] J. Beck, K. Prazdny and A. Rosenfeld. A theory of textural segmentation, Human and Machine Vision, J.Beck, B. Hope and A. Rosenfeld eds., New York Academic Press, 1983, 1-38.

[Cae] T. Caelli. Three processing characteristics of visual texture segmentation. Spatial Vision, 1(1), 19-30, 1985.

[CGG] Y-G. Chen, Y. Giga and S. Goto. Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations, Preprint, Hokkaido University, 1989.

[CIL] M.G. Crandall, H. Ishii and P.L. Lions. User’s guide to viscosity solution of second order partial differential equation. Preprint, CEREMADE, 1990.

[Enn] J. Enns. Seeing textons in context. Perception and Psychophysics, 39 (2), 143-147, 1986.

[EvSp] L.C. Evans and J. Spruck. Motion of level sets by mean curvature, I., preprint.

[FMZ] D.Forsyth, J.L.Mundy, A.Zisserman, al. Invariant descriptors for 3-D object recognition and Pose. IEEE Transaction of pattern analysis and machine intelligence. Vol 13 N.10 1991.

[FRKV] L. Florack, B. ter Haar Romeny, J.J. Koenderink and M. Viergever. Scale and the differential structure of images. Image and vision computing. Vol 10, n. 6 July/August 1992.

[FRKV] L. Florack, B. ter Haar Romeny, J.J. Koenderink and M. Viergever. General Intensity transformations and Second Order Invariants. Procedings of the 7th Scandinavian Conference on Image Analysis. Aalborg, august 13-16, 1991.

[GaHa] M. Gage and R.S. Hamilton. the heat equation shrinking convex plane curves. J. Differential Geometry 23, pp. 69-96, 1986.

[GBS] M. Graham, J. Beck and A. Sutter. Two nonlinearities in texture segregation. Proceedings of the ARVO conference, p161, 1989.

[Gray] M. Grayson. The heat equation shrinks embedded plane curves to round points. J. Differential Geometry 26, pp. 285-314, 1987.

[JGV] B. Julesz, E.N. Gilbert and J. D. Victor. Visual discrimination of textures with identical third order statistics. Biological Cybernetics, 31, 137-140, 1978.

[JuBe] B. Julesz and J.R. Bergen. Textons, the fundamental elements in preattentive vision and perception of textures. Bell System Tech. J. 62(6) 1619-1645, 1983.

[JuKr] B. Julesz and B. Kroese. Features and spatial filters. Nature, 333, 302-303, May 1988.

[Jul1] B. Julesz. Textons, the elements of texture perception, and their interactions. Nature, 290, 12 March 1981.

[Jul2] B. Julesz. Texton gradients : the texton theory revisited. Biol. Cyber. 54, 245-251, 1986.

[Kim] B.B. Kimia. Toward a computational theory of shape. Ph. D. Dissertation, Department of Electrical Engineering, McGill University, Montreal, Canada, August 1990.

[KLM] G. Koepfler, C. Lopez, J.M. Morel. A multiscale algorithm for image segmentation by variational method . Accepted at SIAM journal of Numerical Analysis.

[KoDo1] J.J. Koenderink and A.J. van Doorn. Dynamic shape. Biol. Cyber. 53 : 383-396, 1986.

[KoDo2] J.J. Koenderink and A.J. van Doorn. Representation of Local Geometry in the Visual System. Biol. Cyb. 55 : 367-375, 1987.

[Koe1] J.J. Koenderink. Solid Shape. MIT Press, cambridge, MA 1990.

[Koe2] J.J. Koenderink. The brain a geometry engine. Psychol. Res., Vol 52 (1990) pp. 122-127.

[KTZ] B.B. Kimia, A; Tannenbaum, and S.W. Zucker. On the evolution of curves via a function of curvature, 1: the classical case. To appear in J. of Math. Analysis and Applications.

[LSW] Y.Lamdan, J.T.Schwartz and H.J.Wolfson. Object recognition by affine invariant matching. In Proc. CVPR 88, 1988.

[MaMo1] A. Mackworth and F.Mokhtarian. Scale-Based description and recognition of planar curves and two-dimensional shapes. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 8, No. 1, January. 1986.

[MaMo2] A. Mackworth and F.Mokhtarian. A theory of multiscale, curvature-based shape representation for planar curves. IEEE Trans. Pattern Anal. Machine Intell. 14, pp 789-805, 1992.

[MaPe] J. Malik and P. Perona. Preattentive texture discrimination with early vision mechanisms. Draft, Department of Electrical Engineering and Computer Science, University of California, Berkeley, CA 94720, June 1989.

[Mara1] P. Maragos : tutorial on advances in morphological image processing and analysis. optical engineering, july 1987, Vol. 26 No. 7.

[Marr1] D. Marr. Analysing natural images : a computational theory of texture vision. Cold Spring Harbor Symposium on Quantitative Biology, XL, 647-662, 1976.

[OsSe] S. Osher and J. Sethian. Fronts propagating with curvature dependent speed : algorithms based on the Hamilton-Jacobi formulation. J. Comp. Physics, 79, 12-49, 1988.

[SaTa1] G. Sapiro and A. Tannenbaum. On affine plane curve evolution. Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 821. February 1992, submitted.

[SaTa2] G. Sapiro and A. Tannenbaum. Affine shortening of non-convex plane curves . Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 845. February 1992, submitted.

[Ser] J. Serra Image analysis and mathematical morphology, Vol 1. Academic Press 1982.

[Tre] A. Treisman. Preattentive processing in vision. Comput. Vision, Graphics, Image Process., 31(2), 156-177, 1985.

[VoPo1] H. Voorhees and T. Poggio. Detecting textons and texture boundaries in natural images. Proc. Intern. Conf. Comput. Vision, IEEE, 250-258, 1987.

[Wei] I.Weiss. Projective invariants of shapes. In Proc. DARPA IU Workshop, 1989, pp 1125-1134.

[WeLi] D. Wermser and C.-E. Liedtke. Texture analysis using a model of the visual system. Proc. 6th Intern. Conf. Pattern Recog., 1078-1080, Munich 1982.

[Wit] Witkin, A. P. Scale-space filtering. Proc. of IJCAI, Karlsruhe 1983, 1019-1021.

[YuPo] A. Yuille and T. Poggio. Scaling theorems for zero crossings. IEEE Transactions on Pattern Analysis and Machine Intelligence, 8, January 1986.

Formalisation et vérification d’une théorie psychophysique: la théorie des textons.

Ce chapitre a été écrit en collaboration avec Jean-Michel Morel et a été soumis aux C.R. de l’Académie des Sciences de Paris.

Abstract. Bela Julesz proposed an (informal) axiomatic of the texture preattentive discrimination capability of human perception. We translate this axiomatic into four formal axioms and prove that the axioms lead to a single multiscale analysis, governed by a partial differential equation. Therefore, we can numerically check wether texton theory is true on discriminable and undiscriminable textures. We obtain a surprising result : texton theory achieves hyperdiscrimination!

Résumé. Bela Julesz a proposé une axiomatique informelle de la discrimination de texture en vision préattentive. Nous traduisons cette axiomatique par quatre axiomes formels et prouvons que ces axiomes conduisent à une seule analyse multiéchelle régie par une E.D.P. Par conséquent, nous pouvons valider numériquement la théorie des textons sur des paires de textures discriminables ou indiscriminables en vision préattentive. Nous obtenons un résultat surprenant : la théorie des textons réalise l’hyperdiscrimination!

4.1 Les axiomes de Julesz.

Selon la théorie de Julesz, la discrimination de textures est liée à ce qu’il appelle ”local features”. Il ajoute: ”those local features that might yield texture discrimination should be very special ones. They should be regarded as the basic elements of pre-attentive perception, and I have named such hypothetical elements ”textons”.

Les axiomes de la théorie des textons sont les suivants :

Axiom 1: Human vision operates in two distinct modes: Preattentive vision. A system that is parallel (independent of number of elements), acting almost instantaneously, without scrutiny, and covering a large visual field, as in texture discrimination. Attentive vision. (...)

Axiom 2: Definition of textons (Restricted to line drawings) -Nonoverlapping line segments (black on white background, or white on black background, with specific length, orientation, width, velocity, binocular disparity, and flicker rate). -Terminators (ends-of-lines) of line segments. -Crossings of line segments. -Blobs, whose shape is unimportant, except for average length, width, and orientation.

Axiom 3: Texture formation, texton density, texton gradients -The area called the Δ-neighborhood, in which differences in texton densities (numbers) are determined and texton-gradients computed, is rather limited, and proportional to the size of the texture elements. (...)

4.2 Formalisation des textons.

Les textons sont clairement définis par Julesz comme des formes. Définissons une forme ou une silhouette comme un ensemble fermé X dont le bord est une courbe de Jordan de ℜ^{2}. La recherche en reconnaissance des formes a défini comme paradigme le concept d’analyse multiéchelle . L’analyse multiéchelle vise à régulariser les formes. Cette régularisation permet d’extraire de manière fiable des caractéristiques locales mais introduit un paramètre d’échelle t, qui mesure la perte de localité et le gain en globalité de la caractéristique extraite (voir [KoDo], [Wit]). En réalité, il est implicite que les concepts d’orientation, de courbure, etc ... mentionnés par Julesz pour caractériser les textons sont des concepts multiéchelle : On parle de courbure et d’orientation en un point et à une échelle donnée. Nous notons T_{t}(X) la forme analysée à l’échelle t. X est identifié à sa fonction indicatrice X(x) = 1 si x ∈ X and 0 ailleurs. Nous appelons analyse multiéchelle une famille d’opérateurs (T_{t})_{t≥0} agissant sur les formes et nous posons X(t) = T_{t}(X).

Nous allons d’abord considérer un principe qui s’appelle en théorie de la vision l’architecture pyramidale. D’après cette hypothèse, T_{t} est déduit de T_{s} pour tous réels vérifiant s ≤ t et T_{0} est bien sûr l’identité. C’est une hypothèse nécessaire en théorie de la vision préattentive (irréversibilité de l’analyse). Afin de formaliser l’architecture pyramidale nous introduisons une famille d’ opérateurs de transition T_{t,s} indexés par 0 ≤ s ≤ t et satisfaisant

[architecture pyramidale] T_{t+s} = T_{t+s,s}T_{s} pour tout 0 ≤ s, 0 ≤ t Nous posons donc l’équivalence informelle/formelle suivante:

Axiom 1 ⟺ [Architecture pyramidale]

Supposons que X et Y soient deux silhouettes et que pour tout x ∈ ∂Y et tout r > 0, on ait X ∩ B(x, 2r) ⊂ Y ∩ B(x, 2r). Supposons, de plus, que l’inclusion soit stricte i.e., ∂X et ∂Y se rencontrent uniquement en x. Alors, nous dirons que la forme X est incluse dans la forme Y au voisinage de x. Le principe qui suit exprime que l’ordre entre courbures est préservé lors de l’évaluation multiéchelle.

[inclusion locale de forme] Si X est incluse dans Y au voisinage de x, alors pour h assez petit, T_{t+h,t}(X) ∩ B(x, r) ⊂ T_{t+h,t}(Y) ∩ B(x, r).

Le principe qui suit exprime qu’il n’y a pas d’orientation privilégiée: un disque reste un disque.

[principe de base] Soit D = D(x, κ) un disque de courbure κ and de centre x. Alors T_{t+h,t}(D) est un disque de rayon ρ(t, h, κ) et de centre x. En outre, la fonction h → ρ(t, h, κ) est régulière.

En bref, on a l’implication informelle/formelle suivante:

Axiom 2 ⟹ [principe de base], [inclusion locale de forme]

( Car l’analyse de forme est locale et isotrope. De plus, la conservation de l’inclusion est due au fait que l’oeil est capable de comparer les angles et les courbures. Ceci est implicite dans la description informelle des textons: Les ”terminators”, ”crossings” et ”blobs” sont des indicateurs de courbure.)

L’analyse multiéchelle vérifie donc le principe d’inclusion

si A ⊂ B alors T_{t+h,t}(A) ⊂ T_{t+h,t}(B). On peut alors, comme l’a fait l’école de ”morphologie mathématique” [Ser], associer à une image u ses ensembles de niveau

X_{a}u = {(x, y), u(x, y) ≥ a}

Alors, remarquant que chaque ensemble de niveau est une union de silhouettes, nous pouvons facilement associer à une analyse multiéchelle de forme une analyse multiéchelle d’image en définissant T_{t}(u) par

[principe morphologique] Pour tout a, t, h et u: X_{a}T_{t+h,t}(u) = T_{t+h,t}(X_{a}u)

Quand Julesz se contente de définir les textons sur des ”line drawings”, il admet implicitement le principe morphologique. ”line drawings” (Axiom 2) ⟹ [morphologique]

Le principe suivant de l’analyse de forme est l’invariance affine. Posons AX(x) = X(Ax) pour toute transformation linéaire A. Nous ajoutons le principe [invariance affine] car la discriminabilité de textures planes est peu affectée par la distance et le changement de position de l’observateur.

[invariance affine]. Pour tout A et t ≥ 0, il existe une fonction t′(t, A) ≥ 0 de classe C^{1} telle que AT_{t′(t,A),t′(s,A)} = T_{t,s}A. Il a été prouvé dans [AGLM] que l’on peut toujours supposer t′(t, A) = t.|detA|^{1/2}. Cette relation signifie que le résultat de l’analyse multiéchelle T_{t} est indépendant de la taille et de la position dans l’espace de la forme analysée. Une transformation linéaire correspond à une anamorphose d’une image plane lorsqu’elle est présentée à l’oeil à une distance arbitraire avec une orientation quelconque dans l’espace. Finalement notons que le texton peut être noir sur blanc ou blanc sur noir, et c’est équivalent :

[invariance par inversion du contraste] T_{t+h,t}(X^{c}) = T_{t+h,t}(X)^{c}

”Black on white or white on black background (Axiom 2)” ⟹ [invariance par inversion du contraste]

4.3 Les équations fondamentales de l’analyse de formes et d’images.

Lorsqu’un point x appartient à une courbe plane de classe C^{2} qui évolue, nous allons noter ẋ la dérivée temporelle de x, qui est un vecteur de ℜ^{2}. On note curv(x) la courbure de la courbe en x.

Théorème 1. (i) Sous les cinq principes (pyramidal, inclusion locale de forme, principe de base, invariance affine, invariance par inversion du contraste), l’analyse multiéchelle de forme est régie par l’équation

ẋ = (t.curv(x))^{1/3}n⃗(x)  (1)

(ii) Si nous ajoutons le principe morphologique, l’analyse d’images associée est régie par l’équation

\frac{∂u}{∂t} = (t.curv(u))^{1/3}.|∇u|  (2)

L’équation (2) a été introduite et justifiée axiomatiquement (d’une façon plus compliquée cependant) dans [AGLM]. L’ équation (1) a été proposée et étudiée dans [SaTa]. En fait, avec cette axiomatique, la démonstration est élémentaire:

(i) Soit X une silhouette. Supposons que T_{t}(X) ait un bord C^{2} qui soit une variété de classe C^{2} au voisinage d’un point x de ∂T_{t}(X). Alors elle admet une courbure κ au point x et nous pouvons considérer un disque sousosculateur et surosculateur, i.e., un disque D de courbure κ − ϵ et un disque D′ de courbure κ + ϵ, tangents à la silhouette en x. En écrivant le principe pyramidal et d’inclusion locale des formes, nous voyons que

T_{t+h,t}(D) ∩ B(x, r) ⊂ T_{t+h,t}(X) ∩ B(x, r) ⊂ T_{t+h,t}(D′) ∩ B(x, r)

Donc, en notant x(t + h) le point de ∂T_{t+h}(X) tel que x(t + h) − x(t) soit parallèle à n⃗(x), nous obtenons

ρ(t, h, κ − ϵ) − ρ(t, 0, κ − ϵ) ≤ (x(t + h) − x(t)).n⃗(x) ≤ ρ(t, h, κ + ϵ) − ρ(t, 0, κ + ϵ)

En divisant par h et en passant à la limite lorsque h tend vers 0 et ensuite en faisant tendre ϵ vers 0, nous obtenons grâce au principe de base l’ équation

ẋ = \frac{∂ρ}{∂h}(t, 0, κ)n⃗

Des calculs algébriques assez simples exploitant l’invariance affine puis l’invariance par inversion du contraste conduisent à l’équation (1) (voir [AGLM]). Ensuite, Il suffit de dériver une fois la relation u(t, x(t)) = λ par rapport à l’échelle t pour obtenir (2).

4.4 Vérification expérimentale de la théorie de Julesz.

Les équivalences précédentes et l’axiome 3 permettent de définir les densités de textons à partir de la courbure et de l’orientation multiéchelle . Soit u_{0}(x) une image et u(t, x) son analyse obtenue avec l’équation (2). Alors la courbure et l’orientation à l’échelle t et la position x sont données respectivement par curv(u)(x) et \frac{Du}{|Du|}(x) qui sont bien définies et calculées par l’équation (2). Si la théorie de Julesz est correcte, n’importe quelle paire de textures indiscriminables devrait avoir des courbures multiéchelles indiscriminables. L’expérimentation numérique nous a permis de constater les faits suivants. Les textures proposées par Julesz et celles proposées dans les publications comme discriminables préattentivement sont facilement discriminables en utilisant les densités de textons. Voyons maintenant une paire indiscriminable classique due à Julesz. Dans la figure 1, apparaissent l’image originale (512x512), l’image analysée à une certaine échelle (t = 7 pixels), les courbures négatives et positives de l’image analysée et le résultat d’une segmentation multiéchelle lorsqu’on décide de garder deux régions. Dans ce cas (et dans toutes les textures indiscriminables proposées dans la littérature), nous réalisons une hyperdiscrimination! ( Le dernier quadrant illustre la nécessité d’extraire une courbure à plusieurs échelles pour l’hyperdiscrimination car la courbure négative de l’image originale ne permet pas de conclure.)

La théorie est-elle donc ”trop bonne” pour rendre compte de la vision humaine? Nous pensons qu’il faut préciser l’axiome 3 en disant que la perception évalue l’échelle significative des textons, proportionnelle au ”Δ-voisinage”.

1) Si les textons sont trop isolés, ils sont discriminés comme des formes car personne ne doute qu’un ”5” et un ”10” ne soient deux formes distinctes!

2) Si les textons sont suffisamment rapprochés, alors on doit supposer que l’analyse multiéchelle va être inhibée [MaPe] dès qu’on dépasse une échelle petite devant le ”Δvoisinage”. Si on lève cette inhibition, ce qui ne pose aucun problème numérique, on obtient l’hyperdiscrimination.

References.

[AGLM] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomes et équations fondamentales du traitement d’images. Rapport 9216, CEREMADE. Université Paris Dauphine, 1992.

[Jul] B. Julesz. Texton gradients : the texton theory revisited. Biol. Cyber. 54, 245-251, 1986.

[KoDo] J.J. Koenderink and A.J. van Doorn. Dynamic shape. Biol. Cyber. 53, 383-396, 1986.

[KTZ] B.B. Kimia, A. Tannenbaum, and S.W. Zucker. On the evolution of curves via a function of curvature, 1: the classical case. To appear in J. of Math. Analysis and Applications.

[MaPe] J. Malik and P. Perona. Preattentive texture discrimination with early vision mechanisms. J. Optical Society of America, May 1990, pp. 923-932.

[SaTa] G. Sapiro and A. Tannenbaum. On affine plane curve evolution. Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 821. February 1992, submitted.

[Ser] J. Serra Image analysis and mathematical morphology, Vol 1. Academic Press 1982.

[Wit] Witkin, A. P. Scale-space filtering. Proc. of IJCAI, Karlsruhe 1983, 1019fig.1 Vérification de la théorie de Julesz et algorithmes hyperdiscriminants.

Auteurs: Lopez - Morel.

Julesz avait synthétisé des textures qui sont indiscriminables pour le système visuel humain et il pensait que cette indiscriminabilité était conforme à sa théorie des ”textons”. En réalité, l’équation fondamentale du traitement d’images, conforme à sa théorie, permet de discriminer ces textures! Premier quadrant (512x512): une paire de textures de Julesz indiscriminables en vision préattentive. Deuxième quadrant: image obtenue après évolution par l’équation à l’échelle de 7 pixels. Troisième et quatrième quadrants: courbure négative et positive à l’échelle de 7 pixels. Cinquième quadrant: segmentation obtenue en appliquant l’algorithme de croissance de régions aux canaux de courbures multiéchelle (Le dernier quadrant illustre la nécessité d’extraire une courbure à plusieurs échelles car la courbure négative de l’image originale ne permet pas de conclure.).

90 Chapitre 4

Numerical Analysis of the fundamental equation of image processing.

This chapter has been written in collaboration with T. Cohignac, F. Eve, F. Guichard, and J.M. Morel and has been submitted to the International Journal of Computer Vision.

Abstract. The fundamental equation of image processing

\frac{∂u}{∂t} = |Du|(curv(u))^{\frac{1}{3}}

has been discovered recently with existence proofs in the sense of viscosity solutions [AGLM-1-2-3]. Another approach, by modelling curve evolution has been introduced in [SaTa] and yields the intrinsic equation \frac{∂C}{∂t} = C_{ss} where s is the ”affine length” along the curve. The second equation describes the motion of the level curves of a solution u of the first equation; therefore both models are equivalent. We introduce a third reformulation for this equation, coming from the ”mathematical morphology” set analysis. This diversity of formulations leads to a rich discussion about numerical schemes: classical finite differences, curve evolution schemes, and set evolution schemes. We compare the numerical results and discuss the new kind of problems arising in their numerical analysis: Indeed, the fundamental equation has been created for achieving multiscale shape analysis. Therefore, geometric, affine and scale invariance properties should be very accurately satisfied. We introduce two new curve evolution schemes for the equation which achieve all invariance properties. We also state (and justify) several convergence conjectures for the new introduced schemes and discuss the experimental results obtained on several shapes. They lead to a crossvalidation of two very different schemes.

5.1 Introduction.

In the following, we consider an image u_{0}(x) defined on the plane and with real values representing the brightness. We also consider smooth Jordan curves C(s) where s denotes either the euclidean parameter, defined up to a shift by the relation |C_{s}| = 1, or the affine parameter, defined by |det(C_{s}, C_{ss})| = 1 ([SaTa]). We shall discuss some numerical schemes for five equations in IR^{2} related with image processing.

The most important equation to be considered is the fundamental equation of image processing (the Affine Morphological Scale Space (AMSS)),

\frac{∂u}{∂t} = |Du|(curv(u))^{\frac{1}{3}} ,  u(0) = u_{0}  (5.1)

It can be interpreted as the motion of the level curves of u by the intrinsic equation

\frac{∂C}{∂t} = C_{ss}  (5.2)

where s is the affine parameter. Indeed, this last equation can be rewritten [SaTa2] as

\frac{∂C}{∂t} = curv(C(s))^{1/3}n⃗(s)  (5.3)

where n⃗(s) is the normal at s to the curve C(s) and curv(C(s)) is the curvature of C at s.

In order to make the discussion more fruitful, we shall also include in the discussion the euclidian analogs of these equations, that is

\frac{∂u}{∂t} = |Du|curv(u),  u(0) = u_{0}  (5.4)

\frac{∂C}{∂t} = C_{ss} = curv(C(s))n⃗(s)  (5.5)

where s is the euclidian parameter. Indeed, at a point x = C(s) where u is smooth and Du ≠ 0, one has

curv(u).\frac{Du}{|Du|} = curv(C(s))n⃗(s) = C_{ss}

where C denotes the level curve of u passing by x. (What s is, length or affine parameter, will be clear in context). The relation between the equations (4), (5) has been discovered by Barles, Osher and Sethian ([Bar], [OsSe]).

Equations (2) and (5) have been solved, with existence, regularity and uniqueness results when the initial Jordan curve is smooth in [Grayson, Gage-Hamilton, Sapiro-Tannenbaum]. A general theory has been made for equations (1) to (5) by using the viscosity solution framework [AGLM-1-2-3]. It also yields existence and uniqueness results and implies the equivalence of both approaches ([EvSp], [CGG] ...). Therefore, (at least) two kinds of schemes are available for the applications in image processing :we can either identify a picture with the set of its level curves and then evolve them jointly by using the curve evolution framework. We also can evolve the image itself by using the PDE formulation. We shall introduce a third numerical approach, inspired from the Mathematical Morphology [Serra]. This approach consists in considering the level sets of the image X_{a} = {u(x) ≥ a} and to handle them by operations which are called ”openings”, ”closings”, ”median filters”, etc.... We shall use recent works by Evans-Spruck, Merriman-Bence-Osher, Barles-Georgelin, Evans, Lasry-Morel which prove that the iteration of these filters tends theoretically towards equation (1) or (4). Before starting with the discussion, let us mention that many ideas for the schemes which we discuss here come from the image processing litterature, but they were in many cases not correctly implemented or interpreted because the subjacent PDE models were not identified. We shall make this assertion more precise in the next sections.

In section 2, we shall start with an explanation of the use of the fundamental equation in image processing. This explanation is necessary, because it leads to ask for new kinds of (long time) invariance properties when defining numerical schemes. These invariance properties are a key to applications, e.g. in shape recognition. Therefore, we shall recall an axiomatic characterization of equation (2) which was developed in [LoMo], simplifying the founding approach of [AGLM1]. This axiomatization leads us to a list of invariance requirements for numerical schemes which will be the basis of our discussion through the paper.

In section 3, we present and discuss classical finite difference schemes, close to the Osher-Sethian formulation.

In section 4, we relate the fundamental equation to a new ”affine morphological” set theory and derive new schemes for these equations.

In section 5, we discuss curve evolution schemes and show why they can be understood as special implementations of morphological schemes.

Section 6 is devoted to the experimental comparison of the numerical schemes and to the precise statement and justification of several conjectures about their convergence.

Aknowledgements: The last author thanks Jean-Michel Lasry and Stanley Osher for valuable conversations. Part of this work was made while the last author was visiting Cognitech-Inc.

PLAN

1. Introduction.

2. The four principles of Multiscale Shape Analysis and the fundamental equation.

2.1 The fundamental equation of Shape analysis.

2.2 The fundamental equation of Image Processing.

2.3 Invariance and stability requirements for numerical schemes of the fundamental equation.

3. Quasilinear finite difference Schemes.

3.1 Luis Alvarez’s scheme.

3.2 Luis Alvarez, Frederic Guichard’s scheme.

4. Morphpological (set evolution) schemes.

4.1 A theoretical scheme.

4.2 Iterated median filters and curvature motion.

5. Purely morphological (curve evolution) schemes.

5.1 Introduction.

5.2 Quasilinear intrinsic diffusion schemes.

5.3 A geometrical curve evolution scheme.

6. Experimental results.

7. Conclusion and conjectures. Appendix

5.2 The four principles of Multiscale Shape Analysis and the fundamental equation.

Let us start with a list of four principles of shape analysis which lead to the fundamental equation (For more details, see [AGLM], [LoMo]). These principles will be essential for establishing requirements for numerical schemes and in order to justify the conjectures.

Define a shape or silhouette as a closed set X whose boundary is a Jordan curve of IR^{2}. We denote by T_{t}(X) the shape analysed at scale t. X is identified with its characteristic function X(x) = 1 if x ∈ X and 0 else. We call multiscale analysis any family operators (T_{t})_{t≥0} acting on shapes and we set X(t) = T_{t}(X).

We first consider a principle which is called in vision theory the pyramidal architecture property. According to this assumption, T_{t} can be computed from T_{s} for any s ≤ t and T_{0} is of course the identity. This is natural, since a coarser analysis of the original picture is likely to be deduced from a finer one without any dependence upon the original picture. Of course the finest picture analysis is the identity. In order to formalize the pyramidal architecture we include T_{t} = T_{t,0} into a family of transition operators T_{t,s} indexed by 0 ≤ s ≤ t and satisfying

[Pyramidal architecture] T_{t+s} = T_{t+s,s}T_{s} for all 0 ≤ s, 0 ≤ t

The shape analysis process is thought of as a flow of the shape through different filters, each associated with a scale t. If we discretize such a process, it yields a serie of operators S_{n} = T_{(n+1)h,nh} applied successively.

Assume that X and Y are two silhouettes and that for some x ∈ ∂Y and some r > 0, one has X ∩ B(x, 2r) ⊂ Y ∩ B(x, 2r). Assume further that the inclusion is strict in the sense that ∂X and ∂Y only meet eventually at x. Then we shall say that the shape X is included in shape Y around x.

[Shape local inclusion] If X is included in Y around x, then for h small enough, T_{t+h,t}(X) ∩ B(x, r) ⊂ T_{t+h,t}(Y) ∩ B(x, r).

This last axiom implies that the value of T_{t+h,t}(X) for h small, at any point x, is determined by the behaviour of X near x. We are allowed to take r infinite. Therefore the shape local inclusion principle also implies that if a shape is globally contained in another, this order is preserved for every scale [MaMo]. Notice that the combination of both pyramidal and shape local inclusion principles is generally called in Image Processing causality [Koe], [Wit], [Kim], [SaTa], [KTZ]. Indeed, the causality assumption in Image Processing states that when the scale increases, no new structure appears in the image. Hummel [Hum] noticed that the causality assumption can be reduced to maximum principle.

Both preceding principles allow, as we shall see, to localize the shape analysis in space and time and we therefore only need to state what the multiscale analysis makes of very simple shapes in order to specify it. So we add a basic ”principle” which will state what happens to disks. As we shall see, disks are somehow the ”basis” on which every shape can be decomposed.

[Basic principle] Let D = D(x, 1/r) be a disk with curvature 1/r and center x. Then T_{t+h,t}(D) is a disk with radius ρ(t, h, 1/r) and center x. Moreover, the function h → ρ(t, h, 1/r) is differentiable with respect to h at h = 0 and the differential is continuous with respect to 1/r.

The ”basic principle” implies that the multiscale analysis behaves in a smooth and isotropic way.

In the following, we set

g(t, 1/r) = \frac{∂ρ}{∂h}(t, 0, 1/r)

Note that g(t, s) is defined for t ≥ 0 and s ∈ IR, now to define g(t, 0) we must assume that lim g(t, 1/r) exists when r tends to +∞ or −∞ . The radius r may be positive or negative, according to the orientation of the normal n(x). In the case where the curve is a Jordan curve enclosing a set X, we take as a convention that n(x) is pointing outside X and the curvature is negative if X is convex at x, positive else. It may seem natural to assume therefore that g(t, κ) is odd with respect to κ. Indeed, this corresponds to the assumption that a black disk on grey background and a white disk behave in a same way.

By the shape local inclusion principle, it is easily seen that the function g(t, κ) = \frac{∂ρ}{∂h}(t, 0, κ) is nondecreasing with respect to κ.

The last shape preserving principle is the affine invariance of shape analysis. Set AX(x) = X(Ax) for any linear map A. The affine invariance can be stated as:

[Affine inv]. There exists a C^{1} function t′(t, A) ≥ 0 defined for any A and t ≥ 0, such that AT_{t′(t,A),t′(s,A)} = T_{t,s}A.

It was proved in [AGLM] that we can always take t′(t, A) = t.|detA|^{1/2} up to a rescaling (See Appendix).

This relation means that the result of the multiscale analysis T_{t} is independent of the size and position in space of the analyzed features. This is very important in the world where we live, since the same plane object X can be seen at very different distances and orientations and therefore with any affine change X(x) → X(Ax). A linear map corresponds to the anamorphosis of a plane image when it is presented to the eye at an arbitrary large distance and with an arbitrary orientation in space.

5.2.1 The fundamental equation of Shape Analysis.

When a point x belongs to an evolving curve, we denote by ẋ the time derivative of x, which is a vector of IR^{2}. By curv(x) we denote the curvature of a curve which is C^{2} at x. Recall that the curvature is defined as the inverse of the radius of the osculatory circle to the curve at x. The curvature is zero if the radius is infinite. By n⃗(x) we denote the vector normal to the curve and pointing towards the concavity of the curve (These definitions for curves are compatible with the ones given for sets.)

Theorem 1. [LoMo](i) Under the three principles pyramidal, local shape inclusion and ”basic”, the multiscale analysis of shapes is governed by a curvature motion equation

ẋ = g(t, curv(x))n(x)  (5.6)

where g is a nondecreasing function. (ii) If we assume [Affine inv.], then the equation of the multiscale analysis is, up to a rescaling :

ẋ = γ(t.curv(x))n(x)  (5.7)

where γ is defined as following: γ(x) = a.x^{1/3} if x ≥ 0 and γ(x) = b.x^{1/3} if x ≤ 0 and a, b are two nonnegative values. (iii) If we add that T_{t}(X^{c}) = T_{t}(X)^{c} [Reverse contrast invariance] then the function g in (i) is odd and a = b in (ii).

By the expression ”governed by”, we mean that the equation must be verified at any point (t, x) of X(t) where the silhouette is smooth enough to give a classical sense to both terms of the equation.

Remark: It has been proved [GaHa], [Gra], [Ang] that the equation (5.6) with g(t, s) = s has smooth solutions. This equation has been introduced in image processing by [KTZ], [MaMo] and [AGLM1,2,3] and [Yui] in different contexts. An early version of an algorithm leading to equation (5.6) is in [KoDo].

The equation (5.7) has been first introduced and axiomatically justified in [AGLM1,2,3]. It has been independently proposed in [SaTa1]. Existence and regularity of the solution are proved in [SaTa2].

5.2.2 The fundamental equation of Image Processing.

Since the multiscale analysis verifies the obvious inclusion principle that if A ⊂ B then T_{t+h,t}(A) ⊂ T_{t+h,t}(B), then we can, as well known in the ”mathematical morphology” school [Ser], [Mara], associate with a picture u the set of its level sets

X_{a}u = {x, u(x) ≥ a}

Then, assuming that each level set is a union of silhouettes, we can simply define T_{t}(u) as [morphological principle] For any a, t, h and u: X_{a}T_{t+h,t}(u) = T_{t+h,t}(X_{a}u)

It is easily seen [LoMo] that if we define a multiscale analysis of images from the preceding analysis of shapes by using [morphological principle], it is governed by

\frac{∂u}{∂t} = g(t, curv(u))|Du|  (5.8)

Conversely, this equation has a unique solution in the sense of ”viscosity solutions” [AGLM]. So a consequence of Theorem 1 above is that there is a single image multiscale analysis, satisfying [pyramidal], [local shape inclusion], [basic principle], [affine invariance] and the [morphological principle]. Therefore, we call the associated equation the fundamental equation of image processing

\frac{∂u}{∂t} = |Du|(t.curv(u))^{\frac{1}{3}} , u(0) = u_{0}

Remark : It is easy to obtain the equations (1) to (5) by rescaling the last equations (t → ϕ(t)).

5.2.3 Invariance and stability requirements for numerical schemes of the fundamental equation.

In order to compare schemes, we shall fix in this section a list of formal requirements which derive from the above developed axiomatic analysis. We have first the properties which imply convergence of the schemes toward the equations, and second the invariance properties of the equations themselves. All schemes will be defined as the iteration of a nonlinear discrete filter, which we call T. We set u^{n} = T^{n}u_{0}, where u_{0} is the (discretized) initial image. All schemes will depend on a scale step Δt and a space step Δx.

We shall consider three different discrete representations of u^{n} :

1. The classical discrete representation on a grid u_{i,j}^{n} = u^{n}(i, j), with 1 ≤ i ≤ N 1 ≤ j ≤ N .

2. ”The morphological representation” of u by a bunch of level sets X_{a} = {u(x) ≥ a}, where a admits, for technological reasons, 255 values.

3. The ”purely morphological representation” of u by its level curves, a sequence C_{k} of Jordan curves discretized by affine finite elements. (Each curve is given by a sequence of vertices, with real components).

Let H_{λ}, R, and A the operators defined by H_{λ}u(x) = u(λx), Ru(x) = u(Rx) where R is an isometry of IR^{2}, Au(x) = u(Ax) where A is a linear application of IR^{2}. Since the schemes to be considered will be very different in structure, we cannot state all following properties in a completely formal way. What they mean in each particular case will be clear in context.

• [Consistency] : We shall say that a scheme is consistent if the discrete operator T approximating the second member of (1) or (2) tends to the differential one, when the steps Δx and Δt tend to 0 in a suitable way.

• [Convergence] : ∀t T^{n}u → T_{t}u a.e. when Δt and Δx tend to 0 in a suitable way, and nΔt → t.

• [L^{∞} stability] : Let c,d be real numbers. If c≤ u^{n}(x) ≤d then c≤ u^{n+1}(x) ≤d.

• [Order Preserving] : If u^{n}(x) ≥ v^{n}(x) for all x in IR^{2}, then u^{n+1}(x) ≥ v^{n+1}(x).

• [Morphology] : For all nondecreasing functions h : T h(u) = h(T u). This means that the level sets are handled independently.

• [Pure Morphology] : There is no interaction between the connected components of level sets. (This is a more restrictive property than [Morphology] and means that the level curves are handled independently).

• [Isom.  invariance] : T R_{θ} = R_{θ}T.

• [Scale invariance] : ∀l, n ∈ IN ∃n′(n, l), T^{n}H_{l} = H_{l}T^{n′}.

• [Affine invariance] : ∀n ∈ IN ∃n′(n, A), T^{n}A = AT^{n′} (for any matrix A with entire coefficients (only for equation 5.1.))

• [Quasilinear] : Let T be a discrete operator approximating a differential operator of u on a grid. (For instance, the second derivative in a direction). T is quasilinear if it can be written as a linear combination of the u_{ij}^{n}, where i,j belong to a fixed stencil. The coefficients of this linear combination may vary.

We add this (not very formal) notion because quasilinear schemes are easy to implement.

In the next sections, we define schemes and discuss which of the above properties is satisfied.

5.3 Quasilinear finite difference schemes.

These schemes are based on a diffusive interpretation of the equation (4). Indeed, if ξ is a unit vector, such that ∇u.ξ = 0, and u_{ξξ} the second derivative of u in the direction ξ, we have :

(4) ⇔ \frac{∂u}{∂t} = u_{ξξ}

This formulation yields an anisotropic diffusion equation. Anisotropic, because we diffuse only in one direction ξ, depending on the gradient. (compare with the isotropic diffusion (\frac{∂u}{∂t} = Δu) i.e. the heat equation).

Equation (1) can be rewritten as \frac{∂u}{∂t} = (∇u)^{\frac{2}{3}} (u_{ξξ})^{\frac{1}{3}} , where (∇u)^{\frac{2}{3}} is understood as the speed of the diffusion.

The schemes which we wish to discuss compute u_{ξξ} by a quasilinear formula in the sense defined above.

5.3.1 Luis Alvarez’s scheme.

The numerical scheme of Alvarez uses the main directions of the grid, given by the vectors

(1, 0), (2, 1), (1, 1), (1, 2), (−1, 0), (−2, 1), (−1, 1), (−1, 2)

At every point and every time there will be diffusion in only one of these directions. So, the algorithm is:

Step 1 : We search the nearest main direction η of ξ.

Step 2 : We assume that u_{ξξ} is equal to u_{ηη}. And we discretize u_{ηη} by

u_{ηη}^{n+1} = \frac{u_{i−1,j}^{n+1} − 2u_{i,j}^{n+1} + u_{i+1,j}^{n+1}}{(Δx)^{2}} , if η = i⃗

u_{ηη}^{n+1} = \frac{u_{i−1,j−1}^{n+1} − 2u_{i,j}^{n+1} + u_{i+1,j+1}^{n+1}}{2(Δx)^{2}} , if η = i⃗ + j⃗

and so on... Then define an implicit scheme by

u^{n+1} − u^{n} = Δtu_{ηη}^{n+1}

Step 3 : Solve the implicit scheme and return to step 1.

This scheme satisfies [Isom. invariance] for rotations of radius \frac{π}{4}, and symmetries with respect to main axis. [L^{∞} stability] is also true. But it does not satisfy [Consistency] nor [Morphology] ! We remark experimentally that the scheme is robust ([L^{∞} stability] ). But, the code has not enough geometrical invariants. For instance, a disk becomes after some time a polygon... This can be checked on the analysis of the spiral (fig. Spiral.img) So, it cannot be used for sharp shape analysis because it creates corners which do not exist in the original shape. We need more rotational invariance and this will be achieved by the next scheme.

5.3.2 Luis Alvarez, Frederic Guichard’s scheme.

In the last subsection, we tried to use an adaptative diffusion (in order to keep [L^{∞} stability] ), now we seek for a good estimate of the curvature and gradient direction. We want to have more isometrical invariants. Let us begin by computing a discretization of curv(u) as rotation invariant as possible. We write curv(u) as a linear combination of the values of u in the 9 points.

(u_{ξξ})_{i,j} = ∑_{−1≤l≤1} ∑_{−1≤k≤1} α_{l,k}^{i,j}(θ_{i,j})u_{i,j} + o(Δx)^{2}

where θ_{i,j} is the angle of the vector ξ. At each point (i,j),we would like to define the weights α_{l,k}^{i,j}, in order to have :

(i) the preceding equality.

(ii) the obvious geometrical invariants of the grid.

(iii) A relation which implies [L^{∞} stability]

∑_{−1≤l≤1} ∑_{−1≤k≤1} |α_{l,k}^{i,j}(θ_{i,j})| − 2 ∗ α_{0,0}^{i,j}(θ_{i,j}) ≤ 0

Now, if we have (i) and (ii) , we cannot have (iii) for all angles θ. Thus, we shall not have [L^{∞} stability] . In order to have so much stability as possible, we change therefore (iii) into

∑_{−1≤l≤1} ∑_{−1≤k≤1} |α_{l,k}^{i,j}(θ_{i,j})| − 2 ∗ α_{0,0}^{i,j}(θ_{i,j}) as small as possible.

At last the functions α_{l,k} must be as smooth as possible for stability.

Properties (i), and (ii) imply that

\{ \begin{matrix}α_{1,0} = α_{−1,0} = λ_{1}(θ) \\ α_{0,1} = α_{0,−1} = λ_{2}(θ) \\ α_{1,1} = α_{−1,−1} = λ_{3}(θ) \\ α_{−1,1} = α_{1,−1} = λ_{4}(θ) \\ α_{0,0} = −4λ_{0}(θ) \end{matrix}

where  \{ \begin{matrix}λ_{1}(θ) = 2λ_{0}(θ) − sin^{2} θ \\ λ_{2}(θ) = 2λ_{0}(θ) − cos^{2} θ \\ λ_{3}(θ) = −λ_{0}(θ) + 0.5(sin θ cos θ + 1) \\ λ_{4}(θ) = −λ_{0}(θ) + 0.5(− sin θ cos θ + 1) \end{matrix}

and

λ_{0}(θ + \frac{π}{2}) = λ_{0}(θ) λ_{0}(θ + \frac{π}{4}) = λ_{0}(\frac{π}{4} − θ) λ_{0}(0) = 0.5 λ_{0}(π/4) = 0.25.

We choose for λ_{0} the simplest trigonometric polynomial which verifies these equalities, and property (iii) as much as possible (see Appendix). So we have λ_{0}(θ) = 0.5 − cos^{2}(θ) + cos^{4}(θ) and we deduce the value of λ_{1}(θ), λ_{2}(θ), ...

We can show that it is equivalent to say that the computation of u_{ξξ} is exact for a set of polynomials of degree as high as possible. We can include all polynomials of degree 3 plus the polynomial : (x^{2} + y^{2})(x sin θ − y cos θ)^{2}. (See Appendix).

Our algorithm is following:

• We estimate the gradient direction θ, and ξ the orthogonal direction using the 9 points of the stencil 3x3.

• Computation of curv(u) as above.

• We proceed, in case of equation

\begin{matrix}(5.4) :  u^{n+1}(x) = u^{n}(x) + Δtu_{ξξ}^{n}(x) \\ (5.1) :  u^{n+1}(x) = u^{n}(x) + Δt(Du)^{2/3}(u_{ξξ}^{n}(x))^{1/3} \end{matrix}

In the case of (5.4), we must choose Δt ≤ 0.5(Δx)^{2} in order to have some stability. Whereas for (5.1), the maximum value of Δt cannot be imposed theoretically. This scheme verifies [Consistency] , because it is exact at order 2. All the geometrical

properties ([Isom.  invariance] , [Scale invariance] , ...) are true for all the

polynomials mentionned above. For other functions, we have these properties only for the affine transforms which leave the grid invariant.

\underline{But, the scheme is not morphological.} And, it does not satisfy [Order Preserving] nor [L^{∞} stability] . However, it is experimentally quite stable. In fact, we do not expect to have [Order Preserving] and [Consistency] in such a quasilinear scheme with a fixed stencil. If you impose both, as noticed in [ALM] and proved by [CrLi], you do not limit the size of the stencil.

5.4 Morphological (set evolution) schemes.

5.4.1 A theoretical scheme.

In this subsection, we discuss theoretical schemes inspired from the set-theory of the Mathematical Morphology School. Denote by C a set of convex sets which is stable either by isometries or by linear maps. For example, we can take for C the set of disks, or the set of ellipses, or the set of triangles,... And, for all t ≥ 0 , we set

C_{t}(x) = {B ∈ C,  area(B) = t^{3/2},  and x is the barycenter of B}

Then, we define two operators

IS_{t}(u)(x) = inf_{B∈C_{t}(x)}(sup_{y∈B}(u(y))

SI_{t}(u)(x) = sup_{B∈C_{t}(x)}(inf_{y∈B}(u(y))

For example, if we take for C the set of the circles, then IS_{t} is the ”closing” of radius t, and SI_{t} is the ”opening” of radius t which are wellknown operators of the mathematical morphology [Serra]. If we impose that C is stable by all linear maps, we obtain ”special” opening and closing which verify [Affine invariance] and which are clearly not idempotent. Consider the schemes

(i) u^{n+1} = IS_{Δt}(u^{n})
(ii) u^{n+1} = SI_{Δt}(u^{n})

It is proved in [LaMo] that the scheme (i) converges to the equation

\frac{∂u}{∂t} = |Du|(curv^{+}(u))^{1/3}

and that the scheme (ii) converges to the equation

\frac{∂u}{∂t} = |Du|(curv^{−}(u))^{1/3}.

In order to approximate equation (1), we alternate affine closing and affine opening and we are led to the scheme:

(iii)  \{ \begin{matrix}u^{n+ \frac{1}{2}} (x) & = & inf_{B∈C_{t}(x)}(sup_{y∈B}(u^{n}(y))) \\ u^{n+1}(x) & = & sup_{B∈C_{t}(x)}(inf_{y∈B}(u^{n+ \frac{1}{2}} (y))) \end{matrix}

The mathematical morphology framework also leads naturally to the scheme

(iv) u^{n+1}(x) = med_{B∈C_{t}(x)}(med_{y∈B}(u^{n}(y)))

If A is a finite set of real numbers, we call med(A) = a any real number such that Card{r ∈ A, r ≥ a} =Card{r ∈ A, r ≤ a}. Of course, the possible a’s make an interval and if we want to specify a, we take for a the middle point of this interval.

All of these schemes satisfy:

[L^{∞} stability] , [Order Preserving] and [Morphology] .

Lasry-Morel [LaMo] have introduced the schemes (i), (ii) and (iii) and proved their convergence. And we conjecture for (iii) and (iv) that convergence is true for equation (5.1), if the set C is stable by affinity for convenient values of Δt and n. Evans [Eva] and Barles [Bar] proved that convergence is true for equation (5.4) and the scheme (iv), when the set C is reduced to one disk with radius t^{1/2} (see more details below). Among the schemes discussed in this paper, the morphological schemes (i)-(iv) discretized on a fixed grid are the only ones to be both morphological and orderpreserving. (Whatever discretization of sup, inf, med on the grid is.)

But, they have not all geometrical properties. They satisfy only the geometrical invariance of the grid and do not satisfy [Scale invariance] . Of course, their numerical implementation ”à la lettre” is problematic. It could be tried however by using a fast Legendre transform recently proposed by Y. Brenier [Bre].

5.4.2 Iterated median filters and curvature motion.

Variants of (4) have been proposed by Koenderink - Van Doorn [ KoDo ] and Merriman-Bence-Osher [MBO]. They define a weighted medium filter which weighs the contribution of points y close to a point x according to a gaussian law of distance. So the algorithm can be rewritten ( when applied to a characteristic function ũ^{0} of a level set)

1. We define the image ũ^{0} by :

\{ \begin{matrix}ũ^{0}(x) = 1 if x belongs to the level set E \\ ũ^{0}(x) = 0 else. \end{matrix}

2. We solve the heat equation, with initial datum ũ^{n}, for a small time Δt, by a convolution with a gauss function. We obtain a new function v^{n}.

3. We set (median filter !)

\{ \begin{matrix}ũ^{n+1}(x) = 1 if v^{n}(x) ≥ \frac{1}{2} \\ ũ^{n+1}(x) = 0 else \end{matrix}

4. We turn back to step 2.

Evans [Eva] and independently, Barles-Georgelin [BaGe], proved that this theoretical scheme is consistent and convergent. Earlier works of Yuille [Yui] and Mascarenhas [Masc] proved the consistency. However, we shall see that its implementation on a fixed grid is problematic. Because of the truncature, and the fixed step Δx of the grid, for little curvature, curv(u) will be implicitely wrong calculated. Indeed, in a grid, if the curvature is small, in order to calculate it, you should diffuse a long time. But this is impossible. We cannot diffuse too much at step 1, if we do not want to miss the evolution of high curvature points. Then, the small curvatures will be implicitly computed as equal to 0, and so we can have u^{n+1} = u^{n}, which stops the evolution. In other words, when Δt tends to 0, the scheme tends to an equation

\frac{∂u}{∂t} = |Du|F (curv(u)) where F is defined by :

F (s) = \{ \begin{matrix}s if s ≥ C \\ 0 if - C ≤ s ≤ C \\ s if s ≤ -C \end{matrix}

and C is a constant depending on the grid step. So consistency is in practice not true; there is a thresholding effect on the curvature. We refer to section 6 for the numerical confirmation.

5.5 Purely morphological (curve evolution) schemes.

5.5.1 Introduction.

For any discretized function u, define its level set of value a as the set of points

X_{a} = {x ∈ IR^{2}; u(x) ≥ a}

We divide the datum image u into a set of level sets S. Purely morphological schemes will handle directly the boundaries of the elements of S which are themselves unions of Jordan curves.

5.5.2 Quasilinear intrinsic diffusion schemes.

We shall first discuss a finite element scheme for equations (2) and (5) based on Mokhtarian, Mackworth’s theoretical scheme. Let us recall the idea of these authors for analyzing a curve. We extract a level curve of the image, and it can be represented by a parametric curve C(s) = (x(s), y(s)) where s is either the parametrization by arc length or the parametrization by affine length.

Equation (2) and (5) are equivalent to

\frac{∂x}{∂t} = \frac{∂^{2}x}{∂s^{2}}  (5.9)

\frac{∂y}{∂t} = \frac{∂^{2}y}{∂s^{2}}  (5.10)

Thus, the equation can be interpreted as an oriented diffusion and we are now ready to display a scheme, inspired from [MaMo] and [BSS].

Step 0 We extract the level curves of the picture and parametrize them by length. We obtain (x(s), y(s)).

Step 1 We solve the equations (5.9, 5.10) for a little step Δt, by convolving the functions x(s) and y(s) with a gaussian function.

After processing, we reconstruct the image by superposing the level sets at scale t.

\underline{So, we have here a theoretical scheme which satisfies} \underline{all invariants except} [Order

Preserving] .

Let us discretize this theoretical scheme in the case of equation (2). We start with a finite element representation of the curve to analyze, C̃^{n}(i), i = 1,...,N. The affine parametrization is given by having the ”affine length” m_{i} of the interval (C̃^{n}(i), C̃^{n}(i + 1)), defined by

m_{i} = |area[C̃^{n}(i − 1), C̃^{n}(i), C̃^{n}(i + 1), C̃^{n}(i + 2)]|^{1/3}

where the area is computed algebraically by

area(a, b, c, d) = det(b − a, c − a) + det(c − a, d − a)

We obtain a reparametrization of C̃^{n} by affine length as

C^{n}(∑_{j=1}^{i} m_{j}) = C̃^{n}(i)

C^{n} is defined on [0, ∑_{j=1}^{i} m_{j}] by affine interpolation. Then, as in the theoretical scheme, we define the C̃^{n+1} from C^{n} by a convolution with a positive kernel (this corresponds to the computation of some barycenter). For simplicity, we can take as kernel the weighted characteristic function \frac{1}{2ϵ}1_{[−ϵ,ϵ]} of the interval [−ϵ, ϵ]. Indeed C̃^{n+1} appears as the set of barycenters of ”chord sets” of C^{n}. We call chord set a set enclosed on the one side by a chord of the curve and on the other by the curve. Roughly speaking, the preceding scheme computes a serie of chord sets with constant area Δt and takes their barycenters as a new polygon. It is intuitively true that the inf (resp the sup) is attained in the Lasry-Morel alternate scheme (iii) when the convex set in consideration is nearby a chord set (of course, if Δt is small enough). Let us now discuss the properties of our scheme. The first one to take in consideration is the fact that it is fully affine invariant, because if |detA| = 1 then the definition of m_{i} and of the barycenter is not affected by a change of C into A(C).

Now, we are facing a dilemma when we want to achieve full scale invariance (i.e. when |detA| ≠ 1). Indeed, we can achieve it by assuming that the affine parametrization of any curve runs between 0 an 1 (Normalization of affine length) and convolving it with a fixed kernel. In this case, we get a full affine invariance; This opens the possibility of successful multiscale shape recognition. Now, this renormalization of affine length changes the equation: an easy calculation shows that in this case, T^{nΔt} should converge, when nΔt tends to t, towards the solution of

\frac{∂C}{∂t} = |Length(C(t))|^{2}C_{ss}

Thus, assuming we want to analyze simultaneously several curves (which is the case of image processing! ) we face a problem of resynchronization. This can be done by returning to an universal time τ, defined by:

\frac{dτ}{dt} = (Length(C(t))^{2}

However, discretization of this last equation may in practice lead to errors in the resynchronization (then, scale invariance would be violated because each curve is correctly analysed, no matter its scale, but we loose accuracy on its time-scale). The other face of the dilemma is to give up perfect scale invariance: we convolve with a fixed kernel the affine curve C parametrized on the interval [0, Length(C)]. Therefore, we keep affine invariance, but loose scale invariance because, for instance, the analysis of the curve C/2 will be rougher than the analysis of C. This must be taken into account in shape recognition tasks. We conjecture, in this last case, that the scheme converges strongly to the solution C(t) of the equation (2). In any case we shall see that the ”desynchronized scheme” is excellent in experiments.

5.5.3 A geometrical curve evolution scheme.

Curve sampling problems.

We have seen that the preceding scheme changes time when the affine length of the curve changes. This problem is enhanced by the unstability of the affine length of a polygon when the discretization changes. Take as an example the case where the polygon is a square. The affine length of a square (and in general, of a polygon) is zero. If we discretize the square as discussed above, with constant euclidian step ϵ, it is easily seen (see figure below) that the affine length of the square is proportional to ϵ^{2/3} ( ϵ^{2} corresponds to the surface of the triangles allowing the computation of the area yielding the affine length). Thus, the affine length of the square can be any value between 0 and its area to the power 2/3. So, there is a high dependence on the number of points involved in the discretization. (Of course, the dependence is also in the position of the points). In conclusion, all invariance properties are true if the sampling has been defined in an unambiguous way. But this is impossible in practice.

The same problem arises in a more acute way in the scheme defined by Bruckstein, Sapiro, Shaked [BSS]. They define the affine length of a polygon as the number of vertices and iterate a convolution kernel on that. Therefore, it is a Mackworth-Mockhtarian scheme.

A geometrical scheme.

In order to avoid the above discussed problems of sampling and synchronization, we shall define an easy to implement version of the alternate scheme (iii) of section 4.1.

Let C^{0} be a Jordan curve, X^{0} the set enclosed by this Jordan curve, and x a point of C^{0} = ∂X^{0}. Let K_{t} be an affine invariant and translation invariant set of convex sets (e.g. all ellipses,...) with area equal to t^{3/2}. Then, we define the t-affine opening O_{t}X^{0} as the set of the barycenters of the elements of K_{t} contained in X^{0}. In the same way, we define the closing F_{t}X^{0} as the complementary of the opening of X^{c}. And we consider the alternate scheme

X^{n+1/2} = O_{t}X^{n}
X^{n+1} = F_{t}X^{n+1/2}

This is the set evolution scheme associated with the convergent image evolution scheme 4.1.(iii).

We shall now define a practical scheme associated with this theoretical scheme. For that, we notice that in order to define the opening, we can restrict the sets in K_{t} to be chord sets, that is, sets bounded on the one side by the boundary and on the other by a chord. This leads to a computationally easier variant. X is a set enclosed by a Jordan curve C(s), let (C(s_{0}), C(s_{1})) be a “chord”, that is, a segment joining two points of C with s_{0} < s_{1}.

We define the Δt-opening as follows. For every s_{0}, we consider the chords (C(s_{0}), C(s)) and we call A(s) the area of X between the chord and the curve C(σ) for s_{0} ≤ σ ≤ s. The function A(s) is nondecreasing and we take s_{1} = sup{s, A(s) ≤ Δt^{3/2}}. We call CH(s_{0}) the chord set associated with the chord (C(s_{0}), C(s_{1})). It has an area less or equal than Δt^{3/2}.

Then the Δt-opening of X is defined by

O_{Δt}X = X \ (∪_{s}CH(s))

In the same way, we can define the closing F_{Δt}X, by F_{Δt}X = (O_{Δt}X^{c})^{c}. Of course, O_{Δt}X may be no Jordan set. Now, if we discretize X and C, we can ensure that O_{Δt}X remains a finite union of Jordan sets. We now define a practical algorithm based on the same principle.

1) C(i) i = 1 ⋯ N is a polygon curve enclosing X.

2) For every i, we find the last j such that the area of the chord set CH(C(i), C(j)) is less than Δt^{3/2}. Then, we define a new vertex C̃(i), contained in the vertex (C(j), C(j+ 1)). C̃(i) is chosen so that either C̃(i) = C(j) or the area of the chord set associated with the chord (C(i), C̃(i)) is Δt^{3/2}. Denote CH(C(i), C̃(i)) this chord set.

The new polygon O_{Δt}C is defined as

O_{Δt}X = X \ (∪_{i}CH(C(i), C̃(i))

It is a polygon or a union of polygons. Of course, the algorithm for computing O_{Δt}X is complex if Δt is large, but very simple if Δt is of the same order as the distance between two consecutive vertices.

We notice that such a method is “selfsampling” since the number of vertices cannot increase. The final algorithm consists in alternating O_{Δt} and F_{Δt} as explained at the beginning of this section.

5.6 Experimental results.

In order to decide which schemes are invariant and robust enough for shape analysis and which is best, we shall simply run them on the same test images (with both simple invariant and complex shapes) and check all invariance properties askable for shape analysis. We aim at very reliable algorithms and cross-validation of algorithms with very different principles is a key to be completely sure : so we shall match the best invariants algorithms on complex shapes.

Our plan follows the plan of the preceding sections. We shall test in order two finite difference schemes (Alvarez, Alvarez-Guichard), a set evolution scheme (Merriman-Bence-Osher) and finally the curve evolution schemes which make a synthesis of Sapiro-Tannenbaum, Mackworth-Mocktarian and Lasry-Morel ideas.

5.6.1 The Alvarez anisotropic diffusion scheme.

We first discussed a nonconsistent stable finite difference scheme, based on the diffusive interpretation of equation (3). We shall get rid immediately of the code, because, while it is the only one of the finite difference schemes to satisfy [L^{∞} stability] its lack of rotational invariance tends to transform a circular shape into a polygon, as it is clearly seen in figure Spiral.img. Furthermore, long range diffusion effects can be clearly seen on the same image and can be related to the nonconsistency of the scheme and has a poor behaviour when fronts are present in the image. Therefore, this code, which is quite useful for image denoising (cf. [ALM]) is not suitable for shape analysis : it creates shape elements which did not exist, and is therefore not ”causal” in the sense discussed in section 2.

5.6.2 The Alvarez-Guichard rotationally invariant finite difference scheme.

Let us now test the almost fully rotationally invariant code of Alvarez-Guichard. To summarize its theoretical properties, this code is consistent but not order-preserving nor L^{∞} stable. We shall check one by one all desirable properties on test images.

• A disk remains a disk (basic principle) : see figure (Disque.img), where a strong diffusion effect is also manifest.

• Affine invariance. We test it by generating a shape and its images by linear maps B with respective eigenvalues (1/2, 2), (1/3, 3) and (1/4, 4). The determinants of these maps being 1, commutation between analysis and B should be verified at any scale. We display in the figure helice.img the corresponding experimental results. First quadrant of (helice.img) : the above mentionned anamorphosis applied to the last shape on the right. Quadrant 2 : these shapes are analysed at some medium scale. Quadrants 3 and 4 : the inverse anamorphosis is applied to the original shapes as well as to the analyzed shapes. Theoretically, all original shapes should be identical (and they are close to it, except for some perceptible grid effect). The analysed shapes should also be identical and they are similar in size and main features, but some (important) inflexion points are lost. This is due to the fact that the shapes are small and thin, and therefore the diffusive effect implies that some information is lost. In any case, we are tempted to consider this experiment as positive.

• Local shape inclusion. It is tested in the figure (Cygne.img), where two shapes are enclosed in each other, the interior one being oscillating. One sees that the inclusion of shapes is quite well verified at any scale until the smaller one disappears.

• Pure morphology. In principle, connected components of level sets should not interact. In this figure, however, imitated from Kanisza’s book La grammatica del vedere [Ka], we see how three of the black squares, which are close to each other, get into contact by a diffusion effect and create at large scale a ”grey bean” which is clearly an artefact. This bad phenomenon, however, implies more a limit in accuracy of the algorithm as a true invalidation.

5.6.3 ”Morphologization” of the Alvarez-Guichard scheme.

We know that the [morphology] property can be restored by using the ”morphological school” idea of running the algorithm separately on each level set of the image and then reconstructing it. Moreover, a very efficient numerical idea of Osher and Sethian [OsSe] is to run the algorithm not on a set, but on some Lipschitz function having this set as a level set, therefore allowing a good behaviour of numerical schemes by avoiding explicit front tracking. Let us now make some tests on simple sets. We can first check again the affine invariance by noting (fig. ellipsemorpho.img) that an ellipse remains an ellipse with the same excentricity, when we run the algorithm on the distance function, and a circle remains a circle (figure disque.img). See also fig. Muller.img.

Morphology and scale invariance. Let us make another experiment to check wether the evolution of a set is not altered when we run it with two different Lipschitz functions. In figure (Cspire. img), we run on a function whose value is 0 on the spiral and 10 outside. This function has roughly speaking a gradient less or equal than 10. In figure (spire.img), we run on a function whose value is 0 on the spiral and 255 elsewhere. Thus the Lipschitz constant is 255. Both above mentionned figures display the evolution of the spiral at the same number of iterations, extracted as the 5 level set of the first function and as the 128 level set of the second. One sees that the evolution of the shape and its final position short before disappearance are close to identical!

5.6.4 Set morphology on a grid.

Here, theory and experimentation diverge : as explained above, Merriman, Bence and Osher conjectured that a particular iterated median filter converges towards the mean curvature motion while computer vision researchers know quite well that a median filter, iterated on a grid, remains steady after some iterations. We display in figure (planchetriangle.img) the result of these iterations of the Merriman-Bence-Osher scheme on both a noisy binary image (above) and a grey level noisy image (below). The good denoising result is certainly imputable to the mean curvature effect (elimination of high curvature). Now, medium and low curvature simply remain steady in the right figures. In these figures, one can look at the smoothed corners of the triangle and rectangle : they correspond to pieces of disks with a radius of five pixels ; roughly speaking the size of the support of the used Gauss function... If we wanted to smooth further the shapes, we should enlarge the gaussian until its support exceeds by much the size of the picture. But then, the analysis would of course loose any local character... We think that no good simple solution can be given to this problem if we keep the grid discretization of the analysed sets. So we pass to curve evolution schemes, which allow an arbitrary accuracy on the position of the curves and therefore of their curvature and speed.

5.6.5 Purely morphological (curve evolution) scheme.

Here again, we check all invariance properties. As a matter of fact, both above defined curve evolution schemes behave in exactly the same way and we shall comment them jointly.

• Shape inclusion principle : Figure (planchecygne). Notice how close this evolution is to the Alvarez-Guichard scheme.

• Local shape inclusion principle : Figure (plancherosace). This principle is not theoretically true for the algorithm. In pratice, it seems to be quite well true.

• Affine invariance : Figure (plancheellipse). An ellipse remains an ellipse with the same excentricity. Figures (plancheaffine) and (planchereconstruite) : we repeat the same experiment as in the case of the Alvarez-Guichard code. The affine invariance result is almost perfect.

• Scale invariance. The algorithm (in its ”desynchronized” version), has been written to be perfectly scale invariant. If we dilate a polygonal shape (that is, an array of points with floating point coordinates) by any floating point matrix, then the discretized affine parameter, the convolution, etc... will be computed in exactly the same way with the same samplings, etc... So the algorithm yields experimentally and in practice exactly the same results. The only problem could arise from ”resynchronization” between different curves. We can look again at (planchereconstruite) to check how well this resynchronization works. If we look at the ”blobs” of the last quadrant, we see that three of them are identical and one is slightly bigger. This difference comes from a grid effect, the anamorphosis having increased the area of the initial fourth shape. This was not perceptible at small or medium scale, but is at larger scale because the multiscale analysis tends to sharpen the size differences of shapes. In any case, we conclude to an almost perfect scale invariance after resynchronization.

Crossvalidation of the ”best schemes”. This crossvalidation can be made in a very simple way. Indeed, the theoretical times in both algorihms are proportional and we can make them in practice very close by ”calibrating” the time so that a circle with radius t disappears in time t. Then, crossvalidation consists in running both algorithms at exactly the same times on a complex shape. The spiral is well adapted, because it has a mix of high and very small curvature and lives a long time before disappearing. The figure (planchespirale) displays the evolution of the spiral at exactly the same times as in the figures (spire.img) and (Cspire.img). Crossvalidation is clear.

5.7 Conclusion and conjectures.

We have defined and discussed in this paper several numerical schemes for the fundamental equation of image processing. Since the equation is new, the schemes are also new but we took the occasion to match for a first time numerical ideas with very different origins (finite difference schemes, mathematical morphology, multiscale shape analysis theory...). We also checked sharp consistency criteria, based on invariance.

Among the discussed schemes, two appear as very reliable because of their crossvalidation : a consistent classical finite differences scheme and an adaptive finite element curve evolution scheme, with nonlocal aspects. The second is, however, almost perfect from the shape analysis viewpoint : it preserves accurately the invariance properties of the fundamental equation on very long times. Of course, its use for complex image analysis will lead to run it on an huge number of curves, but this seems feasible with the actual technology. For many tasks in image analysis, the first scheme is reliable enough.

Let us finish by stating how the preceding schemes could be mathematically proved to converge. For the finite difference scheme, there is not much hope because it is not order preserving and therefore a formal convergence proof towards the viscosity solution is problematic. We conjecture that the curve evolution schemes converge (e.g. for the Hausdorff distance) towards the exact evolution of any initial Jordan curve with finite length. Because of the lack of a shape inclusion principle for the fully discretized scheme, it seems however advisable to restrict this conjecture to the semidiscretized schemes (i), to (iv), because those schemes verify all four invariance and order preserving principles stated in section 2. Since we know that those principles can only lead to the fundamental equation (and its ”affine opening” and ”affine closing” variants), the conjecture seems to be more than sound.

A. Appendix.

Denoting by ξ the direction orthogonal to the gradient of u, |∇u|curv(u) is equal to u_{ξξ}. Suppose now that the direction ξ is known :

ξ = (− sin θ, cos θ)

Then, we can write

u_{ξξ} = sin^{2}(θ)u_{xx} − 2 sin(θ) cos(θ)u_{xy} + cos^{2}(θ)u_{yy}  (5.11)

We want to write u_{ξξ} as a linear combination of the values of u on the fixed stencil 3x3. Of course, the coefficients of the linear combination may depend on ξ.

Because the direction of the gradient (and then ξ) is defined modulo π, by symmetry we must assume that the coefficients of points symmetrical with respect to the central point are the same.

In order to have consistency, we must find λ_{0}, λ_{1}, λ_{2}, λ_{3}, λ_{4}, such that :

u_{ξξ} = \frac{1}{Δx^{2}} (−4λ_{0}u_{i,j} + λ_{1}(u_{i+1,j} + u_{i−1,j}) + λ_{2}(u_{i,j+1} + u_{i,j−1})  (5.12)
+λ_{3}(u_{i−1,j−1} + u_{i+1,j+1}) + λ_{4}(u_{i−1,j+1} + u_{i+1,j−1})) + o(1)

We write

u_{i+1,j} = u_{i,j} + h(u_{x})_{i,j} + \frac{h^{2}}{2} (u_{xx})_{i,j} + o(Δx^{2}),

and the same relation for the other points of the stencil. By feeding (12) with these relations and by using relation (11), we obtain four relations between our five coefficients

λ_{1}(θ) = 2λ_{0}(θ) − sin^{2} θ
λ_{2}(θ) = 2λ_{0}(θ) − cos^{2} θ
λ_{3}(θ) = −λ_{0}(θ) + 0.5(sin θ cos θ + 1)
λ_{4}(θ) = −λ_{0}(θ) + 0.5(− sin θ cos θ + 1)

It remains only one degree of freedom for our coefficients : for example the choice of λ_{0}.

Now, in order to have the [L^{∞} stability] , we must have

4λ_{0} ≥ 2(|λ_{1}| + |λ_{2}| + |λ_{3}| + |λ_{4}|)

Now, since 2λ_{1} + 2λ_{2} + 2λ_{3} + 2λ_{4} − 4λ_{0} = 0, λ_{1}, λ_{2}, λ_{3}, λ_{4} must be nonnegative.

But, because of the relations between our coefficients, it is impossible to assume such a thing, except for particular values of θ(0, \frac{π}{4}, \frac{π}{2}, ...).

Then, for the choice of λ_{0}, we impose the [L^{∞} stability] , for the particular values of θ, when we can have it. We have

λ_{0}(0) = 0.5

λ_{0}(\frac{π}{4} ) = 0.25

Then, we impose some obvious geometric invariance properties.

Symmetry with respect to the central point,

(i) λ_{0}(θ + π) = λ_{0}(θ)

Invariance by rotation of angle \frac{π}{4},

(ii) λ_{0}(θ + \frac{π}{2} ) = λ_{0}(θ)

Symmetry with respect to the axis i+j,

(iii) λ_{0}(\frac{π}{2} − θ) = λ_{0}(θ)

Of course, we seek for a smooth function for λ_{0}(θ). We shall take it C^{1}. Then,

(i) ⇒ λ′_{0}(0) = 0  and  (iii) ⇒ λ′_{0}(\frac{π}{4} ) = 0

Moreover, (i) and (ii) imply that λ_{0} is \frac{π}{2}-periodic, and symmetric with respect to \frac{π}{4}. We shall now study it on [0, \frac{π}{4}].

Let us summarize : we are looking for a function λ_{0} on the interval [0, \frac{π}{4}], such that

\{ \begin{matrix}λ_{0}(0) = 0.5 \\ λ_{0}(\frac{π}{4}) = 0.25 \\ λ′_{0}(0) = 0 \\ λ′_{0}(\frac{π}{4}) = 0 \\ λ_{1}, λ_{2}, λ_{3}, λ_{4} as “positive” as possible \end{matrix}  (5.13)

We remark that for θ in [0, \frac{π}{4}],

λ_{1} ≥ λ_{2} and λ_{3} ≥ λ_{4}

0.25 0.3 0.35 0.4 0.45 0.5 0  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 λ_{0}(θ)  θ Search of the optimal λ_{0} 0.5 − cos^{2}(θ) + cos^{4}(θ) cos^{2}(θ)/2 (1 − sin(θ) cos(θ))/2

But

λ_{2}(θ) ≥ 0 ⇒ λ_{0}(θ) ≥ \frac{cos^{2}(θ)}{2}

λ_{4}(θ) ≥ 0 ⇒ λ_{0}(θ) ≤ \frac{1 − sin(θ) cos(θ)}{2}

So, we cannot find λ_{0}(θ) satisfying both inegalities, since

\frac{cos^{2}(θ)}{2} ≥ \frac{1 − sin(θ) cos(θ)}{2}

Then, if we choose λ_{0}(θ) ≥ \frac{cos^{2}(θ)}{2} we have λ_{4}(θ) very negative. If we take λ_{0}(θ) ≤ \frac{1−sin(θ) cos(θ)}{2} we obtain λ_{2}(θ) very negative.

We prefer to choose λ_{0} between both functions, and then to have λ_{2} and λ_{4} negative, but slightly. (see figure above.)

The periodicity of λ_{0} leads us to choose it as a trigonometric polynomial, of degree as small as possible for numerical simplicity. The simplest polynomial verifing (5.13), and between both functions is

0.5 − cos^{2}(θ) + cos^{4}(θ)

We deduce the others λ’s:

−4λ_{0} = −2 + 4 cos^{2}(θ) sin^{2}(θ)
λ_{1} = cos^{2}(θ)(cos^{2}(θ) − sin^{2}(θ))
λ_{2} = sin^{2}(θ)(sin^{2}(θ) − cos^{2}(θ))
λ_{3} = cos^{2}(θ) sin^{2}(θ) + 0.5 sin(θ) cos(θ)
λ_{4} = cos^{2}(θ) sin^{2}(θ) − 0.5 sin(θ) cos(θ)

We shall find the above coefficients by a faster method.

A polynomial basis approach.

We still consider the stencil 3x3. And, we want now to have an estimate of the second derivative of u in the direction ξ, exact for a set of polynomials. Of course, we want the degree of polynomials as high as possible.

In order to have [Consistency] , we want at least the polynomials of degree two. that gives us six independant relations on our α_{ij}.

We remark that on our stencil x^{3} = x and y^{3} = y and that (x)_{ξξ} = (x^{3})_{ξξ} = (y)_{ξξ} = (y^{3})_{ξξ} = 0. Then, our calcul will be exact for the two more polynomials : x^{3} and y^{3}.

The polynomial x^{2}y and y^{2}x give us two other independent relations.

It remains us only one degree of freedom, we can choose one more polynomial. We choose, in order two have the coefficients defined before, the polynomial

(x^{2} + y^{2})(x sin θ − y cos θ)^{2}

Then our estimation of u_{ξξ}, is exact for all the linear combinations of the polynomials of degree 0,1,2 and 3, plus the polynomial :(x^{2} +y^{2})(x sin θ −y cos θ)^{2}. Of course, the calculus can be false for the other polynomials.

118 Chapitre 5 fig.1 Spiral.img Chapitre 5 119 fig.2 Disque.img 120 Chapitre 5 fig.3 helice.img Chapitre 5 121 fig.4 cygne.img 122 Chapitre 5 fig.5 illusion.img Chapitre 5 123 fig.6 ellipsemorpho.img 124 Chapitre 5 fig.7 disque.img Chapitre 5 125 fig.8 Cspire.img 126 Chapitre 5 fig.9 spire.img Chapitre 5 127 fig.10 Muller.img 128 Chapitre 5 fig.11 planchetriangle Chapitre 5 129 fig.12 planchecygne 130 Chapitre 5 fig.13 plancherosace Chapitre 5 131 fig.14 plancheellipse 132 Chapitre 5 fig.15 plancheaffine Chapitre 5 133 fig.16 planchereconstruite 134 Chapitre 5 fig.17 planchespirale References

[AGLM1] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomes et équations fondamentales du traitement d’images. Report 9216, CEREMADE. Université Paris Dauphine, 1992. Accepted at Arch. for Rat. Mech., July 1992.

[AGLM2] L. Alvarez, F. Guichard, P.L. Lions and J.M. Morel. Axiomatisation et nouveaux opérateurs de la morphologie mathématique. C.R. Acad. Sci. Paris, t.315, Serie I, p. 265-268, 1992.

[AhSc] N. Ahuja and B. Schachter. Pattern Models, Wiley, New York 1983.

[AlGu] L. Alvarez, F. Guichard. Software MEGAWAVE, user’s guide of the mean curvature motion package.

[ALM] L. Alvarez, P.L. Lions and J.M. Morel. Image selective smoothing and edge detection by nonlinear diffusion (II). SIAM Journal of numerical analysis 29, pp 845-866. 1992.

[Ang] S. Angenent. Parabolic equations for curves on surfaces I, II. University of Wisconsin-Madison. Technical Summary Reports, 89-19, 89-24, 1989.

[AsBr] H. Asada and M.Brady. The curvature primal sketch. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 8, No. 1, January. 1986.

[Bar] G. Barles. Remarks on a flame propagation model. Technical Report No 464 INRIA Rapports de Recherche, December 1985.

[BaGe] G. Barles and C. Georgelin. A simple proof of convergence for an approximation scheme for computing motions by mean curvature. Preprint 1992.

[BaSo] G. Barles and P.E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equation, Asymp. Anal., to appear.

[BPH] E.B.Barret, P.M.Payton, N.N.Haag and M.H.Brill. General methods for determining projective invariants in imagery. Computer Vision Graphics Image Proc. Jan. 1991.

[Bre] H. Brezis. Analyse Fonctionnelle, Théorie et Applications. Masson, Paris, 1987.

[BSS] A.M. Bruckstein, G. Sapiro and D. Shaked. Affine-invariant evolutions of planar polygons. Preprint June 1992.

[CCLM] F. Catte, T. Coll, P.L. Lions and J.M. Morel. Image selective smoothing and edge detection by nonlinear diffusion. Preprint, CEREMADE 1990. To appear in SIAM Journal of Numerical Analysis.

[CGG] Y-G. Chen, Y. Giga and S. Goto. Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations, Preprint, Hokkaido University, 1989.

[CGHH] Y-G. Chen, Y. Giga, T. Hitaka and M.Honma : Numerical analysis of motion of a surface by its mean curvature. Preprint, Hokkaido University.

[CIL] M.G. Crandall, H. Ishii and P.L. Lions. User’s guide to viscosity solution of second order partial differential equation. Preprint, CEREMADE, 1990.

[CrLi] M.G. Crandall and P.L. Lions. Personal communication.

[Dia] J.I. Diaz. A nonlinear parabolic equation arising in image processing. Extracta Matematicae, Universidad de Extremadura 1990.

[DuTs] G. Dubek and John Tsotsos. Recognising planar curves using curvaturetuned smoothing. Third Conference on computer vision 1990. IEEE computer Society Press.

[Eva] G. Evans. Convergence of an algorithm for the mean curvature. Preprint 1992.

[EvSp] L.C. Evans and J. Spruck. Motion of level sets by mean curvature, I., preprint.

[FMZ] D.Forsyth, J.L.Mundy, A.Zisserman, C.Coelho, A.Heller and C.Rothwell. Invariant descriptors for 3D object recognition and Pose. IEEE Transaction of pattern analysis and machine intelligence. Vol 13 N.10 1991.

[Gab] D. Gabor. Theory of communication. J. Inst. Elect. Engr. 93, 429-457, 1946.

[GaHa] M. Gage and R.S. Hamilton. The heat equation shrinking convex plane curves. J. Differential Geometry 23, pp. 69-96, 1986.

[GGIS] Y. Giga, S. Goto, H. Ishii and M.-H. Sato. Comparison principle and convexity preserving properties for singular degenerate parabolic equation on unbounded domains. Preprint, Hokkaido University, 1990.

[GiGo] Y. Giga and S. Goto. Motion of hypersurfaces and geometric equations. J. Math. Soc. Japan, Vol. 44, No. 1, 1992.

[GKW] G.H. Granlund, H. Knutsson and R. Wilson. Anisotropic non stationary image and its applications. Picture Processing Laboratory, Linkoeping University S581 83 Linkoeping, Sweden.

[Gray] M. Grayson. The heat equation shrinks embedded plane curves to round points. J. Differential Geometry 26, pp. 285-314, 1987.

[Hart] R.Hartshorne. Foundations of Projective Geometry. Benjamin, 1967.

[HoNo] K. Hollig and J.A. Nohel. A diffusion equation with a nonmonotone constitutive function. Proc. on Systems of Nonlinear Partial Differential Equations, John Ball ed., pp 409-422, Reidel 1983.

[Hum] R. Hummel. Representations based on zero-crossings in scale-space. Proc. IEEE Computer Vision and Pattern Recognition Conference, pp. 204-209, 1986.

[Ka] G. Kanisza. La grammatica del Vedere, Il Mulino, Bologna 1980.

[Kim] B.B. Kimia. Toward a computational theory of shape. Ph. D. Dissertation, Department of Electrical Engineering, McGill University, Montreal, Canada, August 1990.

[KoDo] J.J. Koenderink and A.J. van Doorn. Dynamic shape. Biol. Cyber.,53:383-396, 1986.

[Koe1] J.J. Koenderink. The structure of images, Biol. Cybern. 50, 363-370, 1984.

[Koe2] J.J. Koenderink. Solid Shape. MIT Press, Cambridge, MA 1990.

[Koeh] W. Koehler. Psychologie de la forme. Idées nrf-Gallimard, 1964.

[KTZ] B.B. Kimia, A. Tannenbaum, and S.W. Zucker. On the evolution of curves via a function of curvature, 1: the classical case. J. of Math. Analysis and Applications. Vol. 163, No 2, Jan. 15 1992.

[LaMo] J.M. Lasry and J.M. Morel. A monotone consistent theoretical scheme for the fundamental equation of image processing. preprint.

[Lio] P.L. Lions. Generalized solutions of Hamilton-Jacobi equations, Research notes in Mathematics, 69, Pitman Boston 1982.

[LiSr] H-C. Liu and M. D. Srinath Partial Shape Classification Using Contour Matching in Distance Transformation. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 12, No. 11, Nov. 1990.

[LoMo] C. Lopez, J.M. Morel. Axiomatisation of shape analysis and application to texture hyperdiscrimination. Proceedings of the Trento Conference on Surface Tension and Movement by Mean Curvature. De Gruyter Publishers. Berlin 1992.

[LSU] O.A. Ladyzhenskaja, V.A. Solonnikov and N.N. Ural’tseva. Linear and quasilinear equation of parabolic type. American Mathematical Society, Providence, R.I. 1968.

[LSW] Y.Lamdan, J.T.Schwartz and H.J.Wolfson. Object recognition by affine invariant matching. In Proc. CVPR 88, 1988.

[MaHi] D. Marr and E. Hildreth, Theory of edge detection. Proc. Roy. Soc. Lond. B207 1980, 187-217.

[MaMo1] A. Mackworth and F.Mokhtarian. Scale-Based description and recognition of planar curves and two-dimensional shapes. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 8, No. 1, January. 1986.

[MaMo2] A. Mackworth and F.Mokhtarian. A theory of multiscale, curvature-based shape representation for planar curves. IEEE Trans. Pattern Anal. Machine Intell. 14, pp 789-805, 1992.

[MaPe] J.Malik and P. Perona. A computational model of texture perception. Preprint UCB/CSD 89/491, Dept of Electrical Engineering and Computer Sciences, Univ. of California Berkeley 1989.

[Mara1] P. Maragos : tutorial on advances in morphological image processing and analysis. optical engineering, july 1987, Vol. 26 No. 7.

[Mara2] P. Maragos. Pattern spectrum and multiscale shape representation. PAMI, 11(7):701-716, july 1989.

[Marr1] D. Marr Vision. Freeman and Co. 1982.

[Masc] P. Mascarenhas. Diffusion generated motion by mean curvature. Preprint 1992.

[MBO] B. Merriman, J. Bence, S. Osher. Diffusion Generated motion by Mean Curvature. April 1992. CAM Report 92-18. Depart. of Mathematics. University of California, Los Angeles CA. 90024.1555.

[NiSh] M. Nitzberg and Takahiro Shiota. Nonlinear Image Smoothing with Edge and Corner Enhancement. Technical Report n¡ 90-2, Division of Applied Sciences, Harvard University, Cambridge, MA 01238, 1990.

[Nor] K. N. Nordstrom. Biased anisotropic diffusion- A unified Approach to Edge Detection. Preprint. Dept. of Electrical Engineering and Computer Sciences, University of California, Berkeley, CA 94720. (1989).

[OsRu] S. Osher and L. Rudin, Feature-oriented image enhancement using shock filters. SIAM J. on Numerical Analysis, to appear.

[OsSe] S. Osher and J. Sethian. Fronts propagating with curvature dependent speed : algorithms based on the Hamilton-Jacobi formulation. J. Comp. Physics, 79, 12-49, 1988.

[PeMa] P. Perona and J. Malik. Scale space and edge detection using anisotropic diffusion. Proc. IEEE Computer Soc. Workshop on Computer Vision, 1987.

[RaCh] A. Rattarangsi and R. T. Chin. Scale-Based Detection of Corners of Planar Curves. Third Conference on computer vision 1990. IEEE computer Society Press.

[RoTh] A. Rosenfeld and M. Thurston. Edge and curve detection for visual scene analysis. IEEE Trans. on Computers, C-20, 562-569, May 1971.

[SaTa1] G. Sapiro and A. Tannenbaum. On affine plane curve evolution. Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 821. February 1992, submitted.

[SaTa2] G. Sapiro and A. Tannenbaum. Affine shortening of non-convex plane curves . Department of electrical engineering. Technion Israel Institute of Technology. Haifa. Israel. EE Pub 845. February 1992, submitted.

[Ser] J. Serra Image analysis and mathematical morphology, Vol 1. Academic Press 1982.

[SiGi] D. Sinha and C.R.Giardina. Discrete black and white object Recognition via morphological functions. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 12, No. 3, March 1990.

[SMCM] P. Saint-Marc, J.-S. Chen, and G. Medioni Adaptative Smoothing: A General Tool for Early Vision. IEEE Transactions on Pattern analysis and machine intelligence, Vol. 13, No. 6, June 1991.

[Son] M. Soner. Motion of a set by the curvature of its mean boundary. Preprint.

[Wei] I.Weiss. Projective invariants of shapes. In Proc. DARPA IU Workshop, 1989, pp 1125-1134.

[Wer] M. Wertheimer. Untersuchungen zur Lehre von der Gestalt II, Psychologische Forschung, 4, 301-350, 1923.

[Wit] Witkin, A. P. Scale-space filtering. Proc. of IJCAI, Karlsruhe 1983, 1019-1021.

[YuPo] A. Yuille and T. Poggio. Scaling theorems for zero crossings. IEEE Tran-140 Chapitre 5 sactions on Pattern Analysis and Machine Intelligence, 8, January 1986.

[Yui] A. Yuille. The creation of structure in dynamic shape. Proc. of the second international conference on computer vision, Tampa 1988, pp. 685-689.

Vu, le Président  Vu, les suffrageants
Monsieur Y. Meyer  Monsieur G. Aubert

Monsieur M. Barlaud Monsieur J.P. Boujot Monsieur P.L. Lions Monsieur J.M. Morel Monsieur P. Perona Monsieur C. Puech Monsieur M. Rascle

Vu et permis d’imprimer, Le Président de l’Université PARIS IX DAUPHINE

Résumé : Qu’est ce qu’une texture? Il est impossible de répondre à cette question de manière intensive, mais il est tout à fait possible de donner des exemples: un champ de blé, la mer, la peau d’un reptile... Le but de cette thèse est de formaliser mathématiquement la discrimination de textures en prenant comme source la Théorie (psychophysique) des Textons de Julesz.

Nous formalisons les deux aspects essentiels de cette théorie :

D’une part, la définition des caractéristiques géométriques locales, ou textons. Nous montrons en suivant l’axiomatique de Julesz qu’elle ne peut être faite que par l’intermédiaire d’une analyse multiéchelle associée à une équation aux dérivées partielles parabolique et géométriquement invariante.

D’autre part, il faut réaliser le “grouping”: comment segmenter une image à partir des multiples images données par les densités de textons. Nous proposons une méthode variationnelle qui permet d’effectuer cette segmentation.

Enfin, nous faisons l’analyse numérique des deux outils de base (segmentation multicanaux, EDP géométrique) de la théorie de Julesz. Il est alors possible de la vérifier et de la discuter numériquement.