fig 1a- Trois couleurs primaires (lumières)
fig 1b- Trois couleurs primaires (pigments)
Segmenter une image ? Le but de la segmentation est l'extraction d'attributs caractérisant des entités telles que des surfaces homogènes ou texturées, distinguables de leur environnement à la manière des tableaux cubistes. Il ne faut pas confondre ce travail avec celui qui consiste à discerner les formes ou les objets, qui semble pourtant si intuitif, si primitif à notre conscience. C'est l'étape de localisation informelle qui synthétise le perçu visuel avant le "traitement" qui lui donne une consistence cognitive.
On raconte qu'un enseignant avait donné à un élève un travail d'été qui consistait à faire "voir" un ordinateur. C'était il y a un demi siècle, et tous les scientifiques étaient tombés amoureux de leur nouveau calculateur programmable (ordinateur), à tel point qu'ils y voyaient en toute confiance la possibilité de simuler l'intelligence ! A ce prix, la sensation ne devait être qu'une piètre besogne à concocter. L'étudiant s'appelait Gerald Sussman, l'enseignant, Marvin Minsky. Sussman fut écoeuré, Minsky qui a reçu le plus gros budget de l'époque n'a pas réussi à synthétiser son observateur-descripteur. Et on y travaille encore.[1]
La vision informatique est depuis devenu un domaine de recherche indépendant,
et des plus prolifiques.
David Maar a distingué trois niveaux de traitement de l'information visuelle
[2] :
i- extraction de caractéristiques 2D (primal sketch);
ii- mise en relation de ces caractéristiques entre elles en incluant
le point de vue de l'observateur (ébauche 2,5D);
iii- description de la scène en terme d'objets et de relations entre
objets, indépendamment du point de vue (information 3D).
La segmentation représente l'ambition de l'étape (i). La diversité des images et la difficulté du problème ont conduit à l'introduction d'une multitude d'algorithmes que nous n'aurons pas ici l'ambition de résumer (se reporter par exemple à [3]). Les résultats obtenus, évalués souvent de façon très empirique, ainsi que l'évolution des capacités de traitement des ordinateurs, ont souvent amené des scientifiques de profils variés à traiter le problème de manières très différentes.
Historiquement, deux approches duales se sont concurrencées : la première, appelée détection de contours, s'est appliquée à rechercher les discontinuités locales, à discerner leur pertinence, à recomposer leurs connexions. La deuxième, la segmentation à proprement parler (approche par régions), extrait des régions en considérant leur homogénéité vis-à-vis de caractéristiques considérées comme pertinentes, et leur proximité. Si l'objectif est le même, les résultats n'ont pas convergé.
Dans cet article, nous vous proposons de revenir sur le travail de Ohlander,
Price et Reddy [4], qui, en 1977 à Carnegie Mellon University (PA, USA),
ont proposé une des premières méthodes de segmentation
par classification un peu évoluée. Avant d'en détailler
le contenu, nous vous proposons quelques rappels sur les notions de couleurs
et de segmentation.
NB : par commodité, "O-P-R"
fera référence à l'article [4].
[1] CREVIER, A la recherche de l'IA, Champs-Flammarion
[2] HORAUD, MONGA, Vision par ordinateur : outils fondamentaux, Hermès
[3] MAITRE, Le traitement des images, Lavoisier-Hermès
[4] OHLANDER, PRICE, REDDY, in Computer Graphics and Image Processing 8,
Academic Press
La lumière est un support propagatif multidimensionnel particulièrement riche en informations : le vue stéréo photométrique (environ douze niveaux de gris) nous permet à la fois de comprendre les distances, et de savoir les ordonner ; la vue chromatique nous permet d'aller beaucoup plus loin dans la discrimination des éléments de notre environnement (plusieurs centaines de couleurs et d'ombres sont détectées : cf modèle Retinex) : nous pouvons par ce biais ressentir leur nature microscopique (où les photons interagissent avec la matière) !
En analyse d'images, l'exploitation de la donnée chromatique permet
de simplifier effectivement l'identification des objets et leur extraction d'une
scène. Jusqu'aux années 80, la plupart des images "colorées"
étaient en fait des images en pseudo-couleurs (une couleur étant
appliquée à une intensité monochromatique particulière)
: ce n'est pas le cas des images "couleur" dont nous disposons à
l'heure actuelle et dont nous allons étudier le codage. Nous partirons
cependant de l'hypothèse que l'article de O-P-R utilisait
les mêmes considérations de colorimétrie que celles définies
par la CIE (Comission Internationale de l'Eclairage) qui avait posé dès
les années 30 les longueurs d'onde des trois couleurs "primitives".
Un nombre incroyable d'applications ont été développées
autour de l'exploitation systématique des données spectrales.
Nous citerons la télédétection et son application à
l'agriculture de précision (prévision des rendements) ou à
la géothermie (climatologie) (ex: SPOT
du CNES).
Les physiciens utilisent trois caractéristiques pour décrire
les qualités d'une source de lumière polychromatique :
- la radiance (l'énergie émise, en W)
- la luminance (l'énergie qu'un observateur humain perçoit, en
lm)
- la brillance (traduit l'intensité chromatique perçue, subjective).
C'est finalement la grandeur qui nous intéresse.
En tenant compte des propriétés "normales" des yeux humains, il apparait que le spectre perceptible (d'environ 400 nm à 700 nm) est décompasable en trois couleurs primaires : rouge, bleu, et vert. La différence entre les couleurs primaires des lumières et celles des pigments est importante : pour les premières, elles sont additives (superposition) ; pour les secondes, elles sont soustractives (absorption d'une couleur primaire de la lumière, réflection des deux autres), de sorte que l'on considère que les couleurs primaires des pigments sont plutôt le magenta, le cyan, et le jaune (qui sont les couleurs secondaires, c'est-à-dire mélange des primaires deux à deux, de la lumière). (voir figures 1a et 1b)
fig 1a- Trois couleurs primaires (lumières) |
fig 1b- Trois couleurs primaires (pigments) |
Les trois couleurs "primaires" sont :
Rs = 700 nm
Gs = 546,1 nm
Bs = 435,8 nm
Toutes les couleurs C du spectre visuel (ou presque)
peuvent être exprimées dans la base RGB :
C = X.Rs + Y.Gs + Z.Bs
L'application classique de cette décomposition est la télévision.(voir
figures 2a et 2b)
Remarque : une jolie démo à ce sujet existe sur
: http://pages.infinit.net/niuton/ntv/lumino.html
 fig 2a - Coupe d'un écran |
 fig 2b - Luminophores |
Pour une brillance donnée, on obtient le cône de la figure
4.
Les coefficients trichromatiques standarts x, y
et z sont tels que :
; x + y + z = 1
soit
Le cône standardisé donne le diagramme de chromaticité (la fameuse "langue", "horse shoe" en anglais) _ figure 3 : les axes correspondent à x (rouge) et à y (vert).
|
fig 3 - Diagramme de chromaticité (CIE) |
|
|
Un modèle de couleur vise à généraliser et à standardiser la représentation des couleurs. Comme vu précédemment, on se ramène pour le visible à un espace borné 3D : chaque couleur y est représentée par un seul point (x,y,z).
Il existe plusieurs types de modèles ; certains sont privilégiés selon que l'on a une approche "hardware" (écran -> RGB, imprimante -> CMY) ou une approche "manipulation de palettes" (graphisme -> HSV).
- Le modèle RGB (rouge, vert, bleu) explicité plus haut (cf figures 2) est un codage qui sert essentiellement aux tubes lumineux à balayage.
- Nous avons vu que les primaires des pigments étaient le magenta, le
cyan, et le jaune. Pour les imprimantes, on utilisera donc préférentiellement
le modèle CMY.
On a évidemment :
[C, M, Y] = [1, 1, 1] - [R, G, B]
- Le modèle YIQ est le standard pour la télétransmission des signaux TV. De la même manière qu'il s'est agit de trouver une astuce pour rendre compatible les récepteurs radiostéréos avec les récepteurs radiomonos plus anciens, il a fallu conserver le signal d'intensité pour les postes noir et blanc. "Y" correspond donc à la luminance, et les deux autres termes aux composantes chromatiques ("Inphase" et "Quadrature").
Nous partons ici de la grandeur "brillance" ( le V
) : elle traduit indépendament de la fréquence l'intensité
de notre flux lumineux perçu.
Le terme suivant, la saturation ( le S ), explicite l'éloignement
vis-à-vis de la lumière blanche : on dira qu'elle traduit la pureté
relative de notre grandeur.
Le dernier terme, la teinte ou hue, ( le H ) est un
attribu associé à la longeur d'onde dominante dans la distribution
polychromatique perçue. Pour un H correspondant à 700 nm et une
saturation maximale, on aura donc une couleur purement rouge.
NB : la teinte se code de manière angulaire, et est donc périodique.
La saturation et la teinte caractérisent la chromaticité de notre signal : cela est très proche du système perceptif de notre corps. La brillance est une grandeur séparée, ce qui permet de faire des considérations sur les dégradés d'intensité indépendemmant des considérations de couleurs.
Pour bien comprendre l'avantage de ce modèle par rapport au codage RGB par exemple, prenons l'exemple du rehaussement d'un visage humain partiellement ombragé. L'outil idéal pour faire ce genre de translation dans le contraste est le relèvement isotonique d'histogramme : on égalise dans les zones souhaitées les trois histogrammes de R, G, B. Cependant, ces modifications ne se feront pas de manière identique pour chaque couleur, et la resuperposition va donc entraîner un décalage dans le spectre très surprenant, la couleur de la chair prenant par exemple des teintes incohérentes. Dans le modèle HSV, une égalisation de V, voire éventuellement de S, suffisent à résoudre correctement le problème.
Les formules de passage entre RGB et HSV étant complexes, on se rapportera par exemple à [5] pour plus d'informations.
Exemple : |
 Luminosité plus faible |
 Désaturation |
|
 Rotation de teinte (hue) |
[5] GONZALEZ, WOODS, Digital Image Processing, Addison-Wesley Publishing Company
La communauté des chercheurs s'est entendue pour définir une segmentation d'image comme étant un partition de cette image (les parties, disjointes et connexes, étant nommées des régions). Relativement à un prédicat d'homogénéité P donné (en général l'écart min-max, la variance ou la moyenne), une segmentation sera dîte maximale si pour tout (i,j) tels que Ri et Rj connexes, P( Ri U Rj ) => ( i = j )
Ce problème a été correctement posé dans son approche par minisation d'énergie [Mumford et Shah, 1989] : il n'existe pas de solution analytique à ce problème. Les autres approches ne sont donc pas disqualifiées.
Il existe deux approches principales au problème de la
segmentation :
- méthodes de partition de l'espace des luminances (dites "de classification")
: on y exploite ensuite la notion de connexité des régions.
- méthodes (dites "de croissance de régions") exploitant
les relations spatiales des pixels et les associant sur la base de leur connexité
de luminance.
La méthode de O-P-R faisant partie de la première catégorie, nous allons ici résumer les travaux qui ont précédé celui décrit par les auteurs, qui a fait rupture à l'époque. Il s'agit de la technique de segmentation par seuillage.
L'espace des luminances peut être naturellement partitionné en utilisant les niveaux de gris (n.d.g.) des composantes (par exemple R, G, B) de l'image. On associe alors à chaque pixel la classe de niveaux de gris à laquelle il appartient : les ensembles maximaux connexes appartenant à la même classe constituent les régions.
L'illustration classique de ce point de vue est l'image de la mariée à la sortie de l'église, blanc sur fond noir. On pense alors à seuiller l'histogramme entre les tons clairs et les tons foncés. De manière plus générale, le seuillage signifie extraire de l'image I les pixels (k,l) correspondant à un intervalle de n.d.g. : T1 < I ( k, l ) < T2
Il apparait clairement que ce type de discrimination va donc dépendre très précisement de la possibilité de partitionner notre histogramme. Or, il faut se représenter l'histogramme comme la projection sur l'axe des luminances de l'ensembre 3D { k, l, I(k,l) }. L'idée qui sous-tend le seuillage est que l'information de localisation (k, l) n'apporte rien au problème de la segmentation.
Nous vous présentons une série d'exemples partant d'un cas simple qui illustre la difficulté d'appliquer l'idée intuitive du seuillage pour segmenter une image.
|
|
|||||||||||||||||||||
|
||||||||||||||||||||||
Avant de commenter l'ensemble de ces couples d'image/histogramme, nous remarquerons la proximité "visuelle" des six premiers cas au moins. La discrimination visuelle des divers ensembles ne parait pas nous poser de problème, et l'effort cognitif accompli pour repérer les régions nous semble équivalent.
On notera que la proximité "visuelle" des cas 3, 4 et 5 illustre
clairement la difficulté inhérente à attaquer le problème
de la segmentation par l'angle des luminances : nous savons par exemple que
nous pouvons donner à peu près n'importe quelle forme voulue à
l'histogramme de n'importe quelle image par translation isotonique. Il ne faut
donc pas trop compter sur quelques modifications du contraste pour clarifier
notre problème. Pas d'espoir de ce côté là !
Toutefois, comme dit l'expression populaire, “il ne faut pas jeter bébé
avec l'eau du bain”. L'idée de partitionner l'histogramme est une
"idée forte" que nous allons retrouver dans tout ce qui suit
: qu'il s'agisse dans un premier temps de faire une statistique unidimensionnelle
sur une localisation donnée (cas O-P-R décrit
plus loin), ou qu'il s'agisse comme nous le verrons ensuite de faire une statistique
multidimensionnelle qui améliore la manière de considérer
l'interdépendance des informations connues ("clustering" détaillé
encore plus loin).
Comme souligné plus tôt (cas 3), il nous faut réintroduire
la localisation. Cela veut dire définir des seuils "locaux"
qui tiennent non seulement compte de la valeur de la luminance de notre pixel
mais aussi de celles des pixels voisins : T ( k ,l , I(k,l), V(k,l)
), V(k,l) correspondant à une propriété locale
des points du voisinage de (k,l) _ par exemple la moyenne de leur n.d.g.
L'image seuillée S sera donc telle que : S(k,l) = I(k,l) x Test_T(k,l)
où Test_T est la fonction caractéristique qui vérifie si
I(k,l) est inclu dans l'intervalle T(k ,l , I(k,l), V(k,l) ).
Face à une formulation aussi complexe, la première idée
est donc de masquer la notion de localité en créant une nouvelle
image pour laquelle chaque point voit sa luminance pondérée avec
celles de ses proches voisins : c'est typiquement du filtrage par correlation.
On pourra faire par exemple un moyennage : c'est l'équivalent d'un filtre
passe-bas.
Sur l'exemple ci-dessous, nous voyons se profiler trois pics (correspondant
au trois classes) selon la taille du carré de moyennage. Ce n'est évidemment
pas un filtrage optimal. Il est toutefois évident que sans connaissance
préalable du bruit, on ne saura pas précisemment choisir le filtre
passe-bas : or, un filtrage trop abrupt risque d'écraser des détails
déterminants.
Une autre manière plus fine de réintroduire la notion de localité tout en utilisant l'idée d'un seuillage "global" est introduite par O-P-R : nous aurons donc l'occasion de repréciser cela ultérieurement.
On peut modéliser une image (matrice de luminances) comme le produit
de deux caractéristiques :
- l'illumination : quantité de lumière incidente sur la scène
regardée
- la réflectance : quantité de lumière réfléchie
par les objets de la scène.
Pour une image I, on écrira donc : I ( k, l ) = i (k, l) . r
(k, l)
Les intervalles théoriques pour chacune de ces composantes sont : 
Le n.d.g [0,255] correspond donc à 
Il n'est pas pratiquement possible de distinguer i et r partant de I ; cependant, nous allons nous servir de ces deux composantes pour expliciter le filtrage homomorphique.
Soit z(k,l) = ln I(k,l) = ln i(k,l) + ln r(k,l)
alors TF(z) = TF(ln i) + TF(ln r)
On écrira : Z(u,v) = Î(u,v) + R(u,v)
Le filtrage de Z par H donne Zf(u,v) = H(u,v).Z(u,v)
Par TF inverse : zf (k,l) = ITF ( H(u,v).Î(u,v) ) + ITF ( H(u,v).R(u,v)
) = if (k,l) + rf (k,l)
L'image traitée obtenue est donc :
It (k,l) = exp ( zf (k,l) ) = it(k,l) . rt(k,l)
H(u,v) est appelé le filtre homomorphique : il est construit pour traiter
la composante i et la composante r séparément. En effet, l'illumination
est caractérisée par des variations spatiales relativement lentes,
alors que la réflectance tend à varier de façon abrupte,
en particulier au niveau des contours des objets.
On associera donc les basses fréquences de la transformée de Fourier
du logarithme de l'image avec l'illumination.
H(u,v) est construit à l'aide d'un FIR : on fabrique la fenêtre
de Hamming optimale pour laquelle les basses fréquences inférieures
au cinquième de la fréquence de Nyquist sont très aténuées
; H est obtenu par symétrie cylindrique.
On applique ensuite le traitement décrit par le shéma bloc du
filtrage homomorphique : l'histogramme du résultat
devient seuillable. On notera toutefois que les hautes fréquences étant
favorisées, le bruit HF l'est aussi. On aura donc tendance à filtrer
passe-bas avant et après le filtrage homomorphique. Enfin, du "ringing"
apparait au niveau des contours, et cela risque d'induire des artefacts dans
notre processus de segmentation.
Ce prétraitement de l'illumination se justifie en particulier par des
considérations de perception : l'oeil est sensible en premier lieu aux
transitions marquées de contraste (qui correspondent en général
aux bords des objets, par cause du phénomène d'occlusion). Il
se comporte de ce point de vue comme un filtre passe-haut. Cependant, il s'agit
aussi de noter que nous observons des scènes en focalisant rapidement
notre attention sur certaines de leurs parties, démarquées des
autres : l'échelle de considération d'une scène renverrait
alors à une adaptation du seuil du filtre en fonction de la localisation
contrairement à l'équilibrage global que nous opérons.
Il serait plus adéquat d'appliquer un filtrage homomorphique dans les
sous-parties majeures. Il faut ajouter à cela que ce traitement global
présente deux défauts :
- amplification du bruit HF (création de textures indésirables)
- déformation des transitions de contraste ("ringing").
Nous nous plaçons dans le cas (... idéal ? facile ? gentil ?
) où l'histogramme étudié permet de retrouver les classes
des régions de l'image.
On cherche ici à décrire un traitement de l'histogramme permettant
de le partitionner de façon pertinente (déterminer les seuils).
Soit une image qui ne contient de deux régions principales de brillance. L'hitogramme d'une telle image est l'estimée de la fonction densité de probablilité de la brillance : p(z). Cette densité est la superposition de deux densités unimodales correspondant aux deux régions : p(z) = P1 p1(z) + P2 p2(z). On assimilera p1 et p2 à deux densités de probabilité gaussiennes, centrées resp. en z1 et z2 (avec z1<z2), et d'écart type resp. sigma1 et sigma2. (on a evidemment P1 + P2 = 1).
Un seuil T peut être choisi de manière à séparer
les deux modes. L'erreur de classification faite sera alors : (E1
étant la probalilité de se tromper en classant 2 dans 1, et vice
versa)
La minimisation de cette erreur est atteinte lorsque : P1 p1(T) = P2 p2(T)
Le résultat de cette équation passe par la résolution d'un
polynôme du second degré.
Le résultat classique pour une même variance sigma est :
En fait, l'interpolation d'une distribution avec son équivalent gaussien
est encore un problème statistique ouvert (on se reportera par exemple
à la thèse de Guillaume Saint Pierre, Université Paul
Sabatier, Toulouse, 2002).
L'idée première (non optimale) est d'approcher par minimisation
quadratique notre histogramme réel par une superposition de n densités
gaussiennes p(z) : on adoucie l'allure de l'histogramme h(z), on compte le nombre
de maxima ce qui détermine le nombre n de modes recherchés ; enfin
on minimise 
.
Dans notre cas, la résolution de ce problème complexe n'apporte pas grand chose. On place en général le seuil au milieu des pics voisins (cas où la variance est la même). On se souviendra qu'il s'agit d'une approximation très large.
Une image couleur (codée en tri-bande) contient énormement plus d'informations que son équivalent monochromatique (le cube !). L'idée de O-P-R est donc d'exploiter ces informations en les examinant dans un nombre important de bases possibles (c'est-à-dire avec divers points de vue).
Prenons l'exemple de deux nuages de points (clusters) dans l'espace des RGB : l'exemple ci-dessous est construit pour illustrer le fait qu'observés du point de vue des intensités de rouge (c'est-à-dire sans même avoir construit leur statistique), les deux clusters ne sont plus différentiables. Ils le sont parfaitement observés du point de vue des intensités de vert.
A l'époque, comme l'expliquent O-P-R, l'analyse des clusters multi-dimensionnels est encore trop coûteuse d'un point de vue calculatoire. Et même s'ils lui reconnaissent des qualités, O-P-R affirment : " l'analyse d'histogrammes unidimensionnels semble être plus efficace dans le cas général que les variations multi-dimensionnelles, et est adaptée pour beaucoup d'images réelles. Cette méthode est aussi très utile pour la segmentation guidée par l'interprétation où des connaissances extérieures sont disponibles pour sélectionner les régions désirées. Il y a cependant des cas où cette analyse ne fonctionne pas aussi bien que voulue, et il est possible qu'une de ses variations soit mieux adaptée".
O-P-R définissent donc un ensemble de caractéristiques pertinentes, qui sont autant de points de vue différents sur nos données. Ce sont l'analyse des composantes primaires (RGB, YIQ) et leur représentation dans le modèle perceptif (HSV). Neuf histogrammes sont donc construits pour chaque zone traitée (voir les explications sur la notion de "zone traîtée" plus loin). Les histogrammes sont éventuellement légèrement adoucis pour retirer les petits pics.
Les neufs panoramas obtenus sont donc très différents.
L'étape cruciale de la méthode O-P-R est le choix de la composante essentielle (correspondant à la classe dominante) sur l'ensemble des histogrammes de manière à appliquer un seuillage qui va permettre d'extraire les régions correspondant à la classe déterminée. "Usuellement, le meilleur pic est un pic bien séparé d'avec les autres (l'écart entre le sommet et ses vallées est bien marqué)". En fait, O-P-R définissent un ensemble de critères hierarchisés complexes permettant de choisir cette composante essentielle : c'est le résultat de la thèse d'Olhander[6].
Nous avons essayé de décrypter la signification exacte du tableau de description des priorités des caractéristiques discriminantes pour la segmentation. En premier lieu en s'appuyant sur l'article, mais celui-ci renvoie à la thèse (or, nous ne disposions pas d'une copie). En deuxième lieu, en s'appuyant sur la bribe de programme fournie : il est en "SAIL", un vieux langage utilisé en IA à Stanford qui a des fonctions non standards... Leur manière de programmer est confuse, mais il est clair qu'il ne permet pas d'expliciter les priorités. Difficile dès lors de discuter de la pertinence de leur méthode. D'autant plus qu'il apparait que leur quantification est purement empirique. Ceci justifie le fait que nous nous en sommes tenus à dégager la particularité de leur méthode, et à en reconstituer l'esprit avec les outils nouveaux disponibles étant donnée l'évolution des ordinateurs en matière de rapidité et de mémoire.
L'idée avouée de O-P-R était qu'un système de segmentation devrait être capable d'extraire les régions dominantes en premier lieu, et d'utiliser la représentation obtenue pour l'affiner selon les besoins du problème. C'est avec cet objectif que les auteurs ont donc construit une méthode récursive. Cherchant à obtenir une précision réglable, ils ont aussi permis d'aborder la question de l'optimalité de la différentiation statistique des classes de manière tout à fait neuve.
Nous avons vu que dans la plupart des cas, l'étude complète des
modes de l'hitogramme ne nous permettait pas de nous débrouiller pour
déterminer les diverses classes. Il s'agit alors de réintroduire
une notion de localisation dans l'étude même des luminances. Pour
cela, O-P-R introduisent une classification récursive :
1- les pixels déjà attachés à une classe donnée
(niveau dominant de l'histogramme) sont soustraits de l'étude et définissent
des sous-régions dont on étudie la connexité. Les pixels
restants sont de nouveau traités : leurs caractéristiques ne sont
alors "plus" confondues avec celles de la classe dominante maintenant
masquées par sa soustraction.
2- chacune des sous-régions constitue une nouvelle image qui est étudiée
de la même manière.
3- à chaque fois, un seuil de cohérence est testé, à
partir duquel on arrête la récursivité : le réglage
de ce seuil est le point critique de cette méthode.
4- lorsqu'il n'existe plus de zone à segmenter, on étiquette les
régions selon le résultat des masquages successifs.
Le schéma général de l'algorithme de segmentation proposé est donc :
O-P-R devaient, à la fin des années 70, faire face à une plaie qui n'est plus (?) d'actualité : la capacité de calcul des ordinateurs. En ce sens, un des soucis majeurs de l'époque était d'obtenir un résultat rapidement, quand le traitement d'image est souvent très "gourmand". Disons qu'aujourd'hui, nous pouvons remettre à plus tard cette question, voire laisser à un réel spécialiste le soin d'optimiser notre algorithme.
Un des points importants sur lequel insistent O-P-R est l'augmentation de la
rapidité de leur algorithme. Ils appelent cela une "planification".
Il s'agit de ne considérer que l'image réduite de taille (échantillonnée)
de manière à obtenir sans le détail l'aspect général
de la segmentation : de sorte que bien moins de pixels sont à traiter.
Ils exploitent en effet le bénéfice de la récursivité
de la division de régions, qui autorise à préciser le schéma
initial de la segmentation.
Pour être plus précis, ce n'est pas le mascage binaire qui est
élargi car cela n'est pas stable si on échantillonne trop l'image,
mais c'est le seuil qui est conservé et appliqué à l'image
de "taille réelle". Quelques opérations postérieures
sont éventuellement ajoutées.
[6] OHLANDER, Thesis : Analysis of natural scenes, CMU, Pittsburgh (PA,USA).
Ne pouvant se débrouiller de la méthode empirique de détermination du "meilleur pic" dans l'ensemble des histogrammes, nous allons utiliser une variante de la méthode de O-P-R. Plutôt que de considérer des histogrammes mono-dimensionnels, nous allons étudier la proximité des points dans l'histogramme 3D.
L'histogramme 3D n'est pas une statistique : il existe une bijection entre tout point de l'image et tout point de l'histogramme. Cependant, la position du pixel dans l'image n'apparait plus du tout. Comme la plupart des images sont codées en RGB, l'histogramme 3D est ici l'occupation de l'espace RGB par l'ensemble des points. L'occupation de l'espace 3D est à considérer de la même manière qu'un pic dans l'histogramme 1D. De sorte que ce que nous recherchons est la plus grande densité de points dans un voisinage de l'espace donné. Cette approche "physique" par densité, ou par concentration, n'est pas évidente à mettre en oeuvre. On cherchera donc à optimiser une mesure de l'écart entre les points d'un même nuage (un cluster) vis-à-vis de leur barycentre : il s'agit de l'analyse des clusters, aussi appelée l'analyse de segmentation ou analyse de taxonomie.
La méthode traditionnelle pour traiter ce problème est le calcul des "K-means".
L'algorithme K-MEANS est une méthode de partition des données
en K clusters, ou classes. Contrairement à d'autres méthodes dites
hiérarchiques (voir "arbre de clusters"
ci-dessous), l'algorithme des K-MEANS ne crée pas de structure d'arbre
pour décrire les groupements, mais crée plutôt un seul niveau
de clusters.
K-MEANS renvoie une partition des données dans laquelle les objets à
l'intérieur de chaque cluster sont aussi proches que possible les uns
des autres, et aussi loins que possible des objets des autres clusters. Selon
le type de donnée, certaines distances sont plus ou moins adaptées.
Chaque cluster de la partition est défini par ses objets et son centre
(ou centroïde). Le centre d'un cluster est le point qui minimise la somme
des distances aux objets du cluster.
K-MEANS est un algorithme itératif qui minimise la somme des distances
de chaque objet au centre de son cluster, pour chaque cluster. K-MEANS change
les objets de cluster jusqu'à ce que la somme ne puisse plus diminuer.
Le résultat est un ensemble de clusters compacts et clairement séparés,
sous réserve qu'on ait choisi la bonne valeur K du nombre de clusters
que l'on cherche.
Le programme Matlab (pour charger les images, taper : image = imread('nom_image'); puis imshow(km(image,Nombre_de_classes)) )
|
|
 |
 |
|
|
 |
 |
Les resultats sont satisfaisants dans le cas ou les zones sont uniformes : cela marche comme nous le souhaitons sur l'exemple 2, mais les vaches de l'exemple 1 ne sont pas uniformes sur l'image, et la segmentation n'est pas proche à notre avis de celle que donnerait un humain. Un autre inconvenient de la méthode est qu'il faut spécifier àl'avance le nombre de classes...
NB : Une étude plus détaillée des paramètres de cette méthode est consultable sur : http://www.tsi.enst.fr/tsi/enseignement/ressources/mti/kmeans/index.htm
A noter encore que rien ne guarantit la convergence vers le minimum global, et qu'il n'est pas non plus guaranti que l'algorithme produise k clusters, à moins que l'on ne modifie la phase d'allocation pour assurer que chacun des clusters a un nombre de points non nul.
Remarque : à explorer davantage, le Fuzzy C-means ...
Il se trouve que chercher à déterminer exactement le nombre K de classes pertinentes dans notre histogramme 3D revient à résoudre le problème du seuillage optimal présenté plus haut. Nous avons vu que cela n'est pas nécessairement une bonne idée pour les cas où du bruit vient perturber notre signal.
En appliquant l'idée de O-P-R, nous allons donc construire une méthode itérative en ne sélectionnant à chaque fois que les points inclus dans un "meilleur cluster" (à la place de "meilleur pic").
Nous utilisons pour cela l'outil fourni dans MATLAB (statistics toolbox) :
Hierarchical Clustering
(cf http://www.mathworks.com/access/helpdesk/help/toolbox/stats/multiv15.shtml
pour une description complète)
Hierarchical clustering is a way to investigate grouping in your data, simultaneously over a variety of scales, by creating a cluster tree. The tree is not a single set of clusters, but rather a multi-level hierarchy, where clusters at one level are joined as clusters at the next higher level. This allows you to decide what level or scale of clustering is most appropriate in your application.
To perform hierarchical cluster analysis on a data set using the Statistics Toolbox functions, follow this procedure:
- Find the similarity or dissimilarity between every pair of objects in the data set. In this step, you calculate the distance between objects using the pdist function. The pdist function supports many different ways to compute this measurement (default : Euclidean distance).
-> comme pour les K-means, l'idée est donc de représenter l'ensemble de nos points comme un "vecteur" de trois composantes (r,g,b), et de calculer leur distance (par exemple euclidienne même si cela n'est pas naturel pour un maillage discret ; on poura prendre 'cityblock' qui est la distance qui "suit" la quantification) : distance entre 1-2, puis 1-3, ..., puis 2-3, 2-4, etc...- Group the objects into a binary, hierarchical cluster tree. In this step, you link together pairs of objects that are in close proximity using the linkage function. The linkage function uses the distance information generated in step 1 to determine the proximity of objects to each other. As objects are paired into binary clusters, the newly formed clusters are grouped into larger clusters until a hierarchical tree is formed.
-> l'idée est de classer par ordre croissant de proximité nos couples de points, sachant qu'un couple classé constitue une nouvelle entité : la plus proche distance à l'un des points du couple sera attribuée comme la distance à cette entité. L'idée est de propager la recherche du point suivant le plus proche : de là vient l'idée de "hiérarchie".
L'illustration fournie dans l'aide de Matlab est explicite :
Il est possible d'afficher cette hiérarchisation sous forme d'arbre de clusters avec la commande dendrogram. Cette représentation graphique des distances relatives entre chaque entité hiérarchique ne sert absolument à rien du point de vue calculatoire : elle permet seulement de se convaincre qu'il est possible d'établir avec cette approche une manière rigoureuse de distinguer les clusters ; en effet, s'il apparait une longue branche relativement à l'espacement des ramifications inférieures, c'est qu'il devient raisonnable de considérer que le cluster obtenu par l'ensemble des points hiérarchiquement inférieurs est distingable à peu près des autres points. Cette relation entre les longeurs de branches est traduite au travers d'un coefficient de consistence relatif à chaque noeud de l'arbre de clusters. (The inconsistent function compares each link in the cluster hierarchy with adjacent links two levels below it in the cluster hierarchy).
Exemple :
- Determine where to divide the hierarchical tree into clusters. In this step, you divide the objects in the hierarchical tree into clusters using the cluster function. The cluster function can create clusters by detecting natural groupings in the hierarchical tree or by cutting off the hierarchical tree at an arbitrary point.
-> cela permet d'étiqueter directement nos points pour une classe donnée, qui sera soit choisie en fonction de la profondeur souhaitée dans l'arbre de clusters, soit en fonction d'un seuil de sélection de la consistence. Nous utiliserons la première option pour sélectionner seulement le départ de la branche la plus "porteuse" (équivalent du "meilleur pic"). Cependant, la notion de consistance nous sera très utile pour régler la "profondeur de segmentation", c'est à dire le niveau à partir duquel on considérera qu'une région est optimalement segmentée.
Lorsque l'on seuille l'espace des luminances, rien ne guarantit la connexité
des points sélectionnés : on obtient la plupart du temps des nuages
de points dont il s'agit d'étudier la densité spatiale de manière
à supprimer les éléments les plus petits. Si l'on se contente
de compter la nombre de points connexes (surface), on risque de détruire
la majeure partie des régions. La solution, adoptée par O-P-R,
est l'astuce classique qui consiste à réduire-étirer :
pour un masque binaire, on sous-échantillonne au plus proche (moyenne
seuillée à la moitié), puis on reforme l'image en l'agrandissant
à la taille originale.
Sur l'exemple ci dessous, l'image binaire en nuages de points est traitée
avec un facteur 2 puis 4, puis 8.
 |
 |
 |
 |
On relève aisément les principaux défauts de cette méthode
:
- si on obtient bien une reconnexion des proches voisins, on n'arrive pas à
isoler une zone à peu près compacte (quasi-convexe)
- les contours des zones sont rapidement très déformés
: pour se convaincre de l'instabilité du précédé,
on remarquera la différence entre les deux "grosses tâches",
pourtant symétriques originalement.
Le terme de morphologie renvoie plus couramment à une branche de la biologie qui traite des formes et des structures des animaux et des plantes. Dans le contexte de la morphologie mathématique, il s'agit d'extraire des images des composantes utiles pour la représentation et la description de la forme d'une de ses régions, telles que ses limites, son squelette ou sa coque convexe. Cet outil est basé sur la théorie des ensembles : dans le contexte du traitement des images, un ensemble conrrespond à une région. D'abord développé en binaire (ensemble de points appartenant à un intervalle de Z²), cet outil a été étendu pour les ensembles en niveaux de gris ( dans un intervalle de Z²x[0,255] ), plus rarement pour les espaces plus grands (couleurs, temps).
Nous rappelons juste ici que les opérateurs fondamentaux sont la dilatation
et l'érosion. On utilise fréquemment leur combinaison (l'ouverture
étant une érosion suivie d'une dilatation ; la fermeture une dilatation
suivie d'une érosion). On peut ainsi fabriquer un filtre morphologique
qui a tendance à éliminer les petites régions, et effacer
les petits défauts : pour cela, on enchaînera une ouverture puis
une fermeture (voir schéma ci-dessous).
Cette simple opération donne déjà de bien meilleur résultat
que le réduire-étirer (l'exemple suivant est obtenu avec pour
élément structurant un cercle de rayon 1, puis 3, puis 5) : on
notera encore que la symétrie est conservée.
La taille de l'élément structurant correspond à la densité
de points acceptable pour sélectionner un ensemble pertinent : le lien
n'est cependant pas directement évident.
|
 |
 |
 |
On pourra comparer le résultat d'une dilatation itérative conditionelle (remplissage de région), et l'extraction de la coque convexe ; cela risque de ne pas donner des résultats satisfaisants étant donné la répartition aléatoire des nuages de points.
Les meilleurs résulats seront cependant obtenus à l'aide d'une méthode que nous n'avons pas ici pris le temps de mettre en oeuvre : il s'agit de la méthode dite de "ligne de partage des eaux" (LPE) avec marqueurs. En effet, cette méthode permet de sélectionner les régions disposant chacune d'un minimum distinct. Les marqueurs permettent de ne pas considérer les minima qui ne nous intéressent pas. On peut prendre, par exemple, le résultat d'un réduire-étirer que l'on érode pour qu'il ne déborde pas, et l'appliquer comme marqueur à notre algorithme LPE.
NB : Une étude plus détaillée des LPE
est consultable sur : http://www.tsi.enst.fr/tsi/enseignement/ressources/mti/filtrage_lpe/mti1.htm
Pour un exemple d'utilisation, voir : http://www.tsi.enst.fr/tsi/enseignement/ressources/mti/bulles/i3m.html
Nous résumons ici les principales étapes de l'algo proposé
dans la démo
NB : vous retrouverez les explications à l'intérieur même
du code en commentaire.
Son schéma séquentiel est évidemment celui proposé par O-P-R
Etape 1 : on sélectionne l'image à traiter. Nous vous déconseillons de choisir une image de plus de 50 pixels de côté (voir pourquoi au niveau du clustering).
Etape 2 :
La consistance est le paramètre qui sert à sélectionner
la tolérance dans l'écart des luminances pour le choix des clusters.
Le nombre de pixels maximaux correspond à la surface maximale des détails
qui nous intéressent.
Etape 3 : le programme affiche le résultat de la segmentation à
chaque fois qu'une sous-région ne sera plus traitée (sub-divisée).
Cette partie est la boucle à proprement parler du schéma
O-P-R.
Revenons sur ses principaux enjeux :
L'idée première pour analyser l'histogramme 3D en terme de densité est de le corréler avec une boule dont le rayon serait relatif à la tolérance choisie. Cependant, il faut parcourir les 256*256*256 couleurs à chaque fois ! C'est irréaliste.
Nous utilisons donc le "hierarchical clustering". Nous mettons sous
forme de vecteur notre image pour en faire une série de points dans l'histogramme
3D. Cette étape n'est déjà pas vectorisable (cela veut
dire, en langage Matlab, optimiser son programme pour qu'il tourne plus rapidement
en exploitant les "méthodes" de Matlab : typiquement une programmation
de type C sous forme de boucle ira 100 fois plus vite sous forme vectorisée).
Mais la principale limitation en terme de calcul arrive à l'étape
d'après, pour le calcul de l'arbre des clusters :
Y = pdist(X);
Cette étape peut être très longue à calculer. En
effet, pour une image m par n, le vecteur de points traité est de taille
m fois n, et la table Y des distances réciproques est de taille (m.n)x(m.n-1)
: on a donc un calcul en N puissance 4 !! Sur notre machine, la mémoire
n'était pas suffisante pour traiter une image telle que N>50.
Z = linkage(Y,'centroid');
L'arbre de clusters est le rapprochement deux à deux des entités
les plus proches (les entités initiales étant les points, un rapprochement
étant à son tour une entité).
Ici la question du "plus proche" n'est plus une évidence. Avec
Matlab, nous avons différentes manière de choisir ce "plus
proche" : distance minimale ou maximale, distance au centre, distance moyenne
(celle que nous avons choisi parce que cela nous parlait davantage, mais c'est
à examiner de plus près).
Premier test de sélection : la consistance sera ici la distance entre les deux principaux clusters démarqués. Si elle est faible, c'est que la zone est relativement homogène. Si ce test n'est pas validé, on doit examiner les éléments du cluster démarqué en détail.
Ici se pose un problème essentiel qui est la représentativité
du plus petit cluster sélectionné. En effet, un cluster tout petit
peut par exemple se démarquer très nettement dans l'espace des
luminances : il sera donc sélectionné ; or le filtre va le rejeter
(du fait de sa petitesse), donc l'analyse de toute la région sera considérée
comme terminée. Nous touchons ici du doigt une des limitations principales
de l'approche par arbre de clusters : nous avons accès à la structure
de nos points en terme de distance, et non pas en terme de densité locale.
Illustration pour mieux nous faire comprendre : imaginons une image test où
on a deux plans constants séparés par une ligne de luminance très
différente. Il peut être bon de sélectionner cette bande
(par exemple, si on segmente une image aérienne de champs, la mince séparation
qui les distingue nous intéresse), ou au contraire, il peut être
préférable de l'ignorer comme du bruit (cas des lignes blanches
dans les éblouissements des CCD où le drain a saturé).
Nous arrivons donc ici à une limitation de notre algorithme : il faut
avoir un a priori sur l'image à traiter. Si nous considérons que
la ligne nous intéresse (on la garde en tant que sous-région),
cela exclue les images où le bruit ne suit pas une distribution gaussienne
autour de la valeur "normale" (cas pour lequel les spots de "bruit"
seraient conservés). Cependant, nous pourrions utiliser une astuce pour
dépasser cette limitation : celle-ci consisterait à masquer les
petits clusters bien séparés (en les soustrayant du masque), sans
les valider en tant que sous-région, et à recommencer le traitement
puisque le filtrage va combler les "trous" occasionnés (sauf
si la répétition de ce masquage vient à élargir
dramatiquement le trou : voir l'exemple du cône de luminance dans l'étude
du filtrage-> l'astuce sera alors de créer un masque de masque : nous
nous épargnerons ce niveau supplémentaire, mais nous nous souviendrons
que cela peut être la cause d'erreurs).
Le filtrage morphologique (décrit plus haut) que nous utilisons sert à réintroduire la notion de densité. La démonstration en est faite avec l'image bitmap suivante, que nous filtrons avec un élément structurant circulaire de taille 2 (cela va très vite sous Matlab).
 |
 |
Le choix de l'élément est crucial : c'est lui qui définit la forme et l'orientation de l'alignement de densité observé. Si nous cherchons à sélectionner les limites d'un champ dans une vue aérienne, un cercle va les écraser, alors qu'une croix permettra aux alignements verticaux et horizontaux de survivre (c'est une solution, ce n'est sûrement pas la meilleure).
La correlation taille de l'élément structurant / densité n'est pas si simple toutefois. Pousser la taille de l'élément structurant permet de retirer correctement les petites régions, mais le filtre ne prend plus en compte la densité initiale. En fait, l'enjeu du filtrage est de supprimer les problèmes liés aux forts gradients : dès qu'on a une forte variation dans l'espace des luminances, on n'a plus de proximité dans notre histogramme 3D. Si on cherche à retirer le bruit, il faut évidemment écraser les effets de contour, mais les éléments pertinents de faible épaisseur sont perdus... On retombe sur une problématique du filtre de Canny (en rapport avec la largeur de la gaussienne), récurrent à l'époque : adoucir les forts gradients ou conserver le bruit...
Pour cela, il faut favoriser l'élément le plus petit, et examiner à la fois sa consistance et sa taille, de manière à détecter les pièges qui consistent à sélectionner une seule forme, à la faire disparaitre par filtrage, et à renvoyer donc la meme sous-région à traiter (boucle infinie).
Illustration : prenons l'image d'un cône de luminance.
La pointe va être correctement sélectionnée. Pour les cercles
juste inférieurs, le filtrage morphologique va retirer le vide intérieur
: on va aussi les sélectionner correctement. Pour les cercles suivants,
facilement détectés dans l'espace des luminances, si la pente
du cône est trop forte, "on a perdu" car le filtre les fait
disparaitre...
Pour résumer, l'algorithme est le suivant :
numéro := 0;
pour tous les pixels P de gauche à droite et de haut en bas :
si P est différent du fond
alors si seulement un des deux voisins sup et gauche a une étiquette
e
alors marquer P avec e ;
sinon si ces deux voisins ont la même étiquette e
alors marquer P avec e ;
sinon si ces deux voisins ont deux étiquettes différentes e1 et
e2
alors marquer P avec min(e1, e2) ;
mettre (e1, e2) dans une table d'équivalence ;
sinon numéro := numéro+1;
marquer P avec numéro;
fin si
fin si
fin pour
prendre le plus petit élément comme représentant de chaque
classe;
Donnons l'exemple suivant :
Le premier balayage donne un résultat de la forme (ici deux images pour
faire ressortir les premières et les dernières étiquettes)
:
Il s'agit donc réordonner selon la table d'équivalence : la première
donne la deuxième (sur fond artificiellement dégradé pour
la visualisation).
Sur la base de cette décomposition, nous testons la taille de la sous-région obtnue ; passé un certain seuil, un nouveau masque est défini qui est renvoyé dans la pile des sous-régions à traiter.
La première version de notre algorithme (sous forme C) prenait plus de 120 secondes pour traiter 120000 pixels. Sa version vectorisée (forme Matlab) en prend moins d'une ! Le facteur limitant de ce programme est bel et bien la construction de l'arbre de cluters !
Compte-tenu de la limitation de notre programme, il nous apparait inutile de multiplier les exemples.
Nous présenterons le suivant :
Les étapes successives de la segmentation donne (consistance=0.2 ; surface_min=16)
:
Le résultat final donne très correctement :
Il s'agit donc de rechercher une méthode plus pragmatique de clustering
que celle proposée : son mécanisme est désuet, nous recherchons
en effet davantage un outil qui produise des régions utiles pour une
application non supervisée (choix automatique de la consistance ?). Ici,
le choix du seuil de "consistance" (qui correspond à l'homogénéité)
se fait par l'observation empirique des résultats obtenus pour diverses
valeurs. Rien ne guarantit qu'un seuil valable pour une image soit correct pour
une autre : c'est à explorer.
Nous remettons ce travail à vos bons soins.
Une autre amélioration notoire sera d'utiliser l'astuce de "PLANNING" (O-P-R) : nous vous proposons son corrolaire sous Matlab, qui consiste à utiliser un traitement prélimianire par k-means (plus rapide, en N puissance 2), de manière à obtenir des sous-masques de taille suffisament réduite pour être calculables dans les limitations de mémoire de votre machine.
Nous signalons ici qu'une définition optimale du point de vue perceptif de ce qu'est un cluster est donnée par Jean-Michel Morel, Agnès Desolneux et Lionel Moisan dans leur "Théorie de la perception et géométrie stochastique". Cette formulation mathématique rigoureuse ne nous a pas donné de solution d'implémentation simple du clustering.
On notera enfin que la segmentation des images en couleur s'avère être bien plus pertinente que le modèle de base en N&B dans la mesure où un grand nombre d'informations sont apportées par la chromaticité. Le pas suivant à franchir est l'utilisation d'une série temporelle d'images (un film) permettant de détecter les changements dans une scène donnée. Si vous voyez cette image pour la première fois ...
 |
||
 |
 |
 |
 |
||
...vous mettrez sans doute un certain temps avant d'arriver à distinguer
le dalmatien. La même image animée nous vous aurez posé aucun
problème. Les applications les plus courantes sont la soustraction d'arrière-plan,
et la détection automatique des limites d'une séquence.
Plus poussées encore, certaines méthodes exploitent les textures.
Daniel CREVIER, A la recherche de l'IA, Champs-Flammarion
Rafael GONZALEZ, Richard WOODS, Digital Image Processing, Addison-Wesley Publishing Company
Radu HORAUD, Olivier MONGA, Vision par ordinateur : outils fondamentaux, Hermès
Ron OHLANDER, Keith PRICE, Raj REDDY, Picture segmentation using a recursive region splitting method, Computer Graphics and Image Pocessing 8, 313-333 (1978)
Henri MAITRE, Le traitement des images, Lavoisier-Hermès
NB : tous les exemples et les codes MATLAB sont disponibles
en format WINZIP en cliquant sur l'icône
L'image intégrale et originale de cet exposé est disponible en
format WINZIP en cliquant sur l'icône 
:
l'utilisation en est libre à condition de siter la source. Merci.
Les auteurs : |
Projet réalisé dans le cadre du DEA MVA 2003/2004 http://www.cmla.ens-cachan.fr/Cmla/DeaMVA/ |
|
|
 |
|
