Partie obligatoire
TME 8: MMC et structure 2D
Dans ce TME, nous allons chercher à analyser automatiquement des images. La base est issue de lien, mais nous allons travailler sur un pickle pré-formaté pour gagner du temps.
Pour ceux qui le souhaite (en dehors de la séance de TME), le script de mise en forme des données est disponible ici.
Chargement des données
Après avoir récupéré le pickle : lien (pour ceux qui le souhaite, un second jeu de données au même format est disponible ici pour un test sans biais), charger le fichier:
data = pkl.load(file("fichierData.pkl","rb"))
Dans les buts de préserver votre quota d'espace disque et d'accélérer le chargement des données, n'hésitez pas à stoker le pickle dézippé dans le /tmp de la machine.
Il suffit de d'ajouter le chemin d'accès et le tour et joué.
Dans ce fichier se trouvent plusieurs champs accessibles avec la syntaxe suivante:
# chaque champ contient les caractéristiques de 100 images:
# dans chaque champ, on peut acceder à l'image i en faisant par exemple data['im'][i]
images = data['im'] # 100 im : 256x256 pixels x RGB
label_pixel = data['lab'] # 100 im : 256x256 pixels
# pour le TME, nous travaillerons sur des segments:
segments = data['seg'] # 100 im : 256x256 pixels => 16x16 segments
label_segments = data['segLab'] # 100 im : 256 segments x labels
coord_segments = data['coord'] # 100 im : 256 segments x 4 coordonnées dans l'image
description_segments = data['feat'] # 100 im: 256 segments x 9 descripteurs RGB
# pour plus de détail sur les descripteurs, voir les questions suivantes sur le modèle vectoriel
graphes = data['graph'] # 100 im: 256 segments x connexions g,d,h,b
labels = data['labMeaning'] # signification des étiquettes
#[u'window', u'boat', u'bridge', u'building', u'mountain', u'person', u'plant', u'river', u'road', u'rock', u'sand', u'sea', u'sky', u'sun', u'tree']
Ce qui correspond du point de vue graphique à:
Et au graphe:
Pour récupérer les fonctions d'affichage: aff
Modèle vectoriel standard
Cette partie n'étant pas le but premier du TME, la plupart du code est fournie. Nous nous reposerons sur la bibliothèque python sklearn (SciKitLearn). Normalement l'installation de la PPTI fonctionne bien. En cas de problème (certaines installation repose sur une version un peu ancienne de la bibliothèque):
pip install scikit-learn --user --proxy=proxy.ufr-info-p6.jussieu.fr:3128 --upgrade
Pour chaque segment, vous disposez de 9 caractéristiques (assez pauvres) correspondant à :
- histogramme de la distribution des valeurs de Red sur 3 niveaux
- histogramme de la distribution des valeurs de Green sur 3 niveaux
- histogramme de la distribution des valeurs de Blue sur 3 niveaux
Focus et discussion sur les caractéristiques:
Nous considérons les pixels d'une région de l'image, codés en RGB. Les valeurs des couleurs sont données entre 0 et 255: nous avons créé 3 intervalles de 0 à 255/3 (niveau faible), de 255/3 à 2*255/3 (niveau intermédiaire), et la fin de plage de définition pour les valeurs fortes. Chaque histogramme somme à 1 (car il modélise une distribution).
Ainsi par exemple, le vecteur de caractéristiques: {$0.98,0.02,0. ,0. ,1. , 0. , 0. , 0. ,1.$} désigne une région où le rouge est virtuellement absent, le vert toujours à un niveau intermédiaire et le bleu toujours au max. La région est donc majoritairement bleu-turquoise.
Discussion: ces descripteurs très simples ne sont pas les plus performants, il ne prennent pas en compte la texture. Même du point de vue colorimétrique, il est souvent plus intéressant de passer de RGB en HSV: n'hésitez pas à faire des tests en dehors du TME. Dans le cadre de la séance, nous cherchons seulement à mettre en évidence la distinction entre un classifieur de base et un classifieur exploitant la structure de l'image.
Pour accéder à la représentation d'un segment s de l'image i:
s=2
features = data['feat'][i][s] # => 9 valeurs
# array([ 0.98046875,0.01953125,0. ,0. ,1. , 0. , 0. , 0. ,1.], dtype=float16)
# pas de rouge, vert à un niveau moyen et bcp de bleu => turquoise
Etape 1: construction d'un jeu de données standard X,Y
X = np.array([ x for i in range(len(data['feat'])) for x in data['feat'][i]])
# Y : l'étiquette correspondante
Y = np.array([ x for i in range(len(data['segLab'])) for x in data['segLab'][i]])
Etape 2: apprentissage d'un classifieur
Nous allons utiliser un algorithme de regression logistique existant (pas besoin de récupérer votre TME 5):
import sklearn.naive_bayes as nb
classifier = lin.LogisticRegression();
classifier2 = nb.MultinomialNB(); # pas propre sur des données continues... Mais efficace
# apprentissage
classifier.fit(X,Y) # automatiquement multi-classes un-contre-tous
classifier2.fit(X,Y)
# inférence sur 1 ou plusieurs vecteurs
print classifier.predict(X[0]) # rend une classe
print classifier.predict_proba(X[0]) # rend le vecteur de proba d'appartenance aux classes pour X[0]
print classifier.predict_log_proba(X[0]) # log proba
Pour la documentation de sklearn: lien
Afin de comprendre le fonctionnement des modèles linéaires, il est amusant d'aller récupérer leur coefficients et de les analyser. La fonction suivante permet de faire ça, elle vous apprendra que la mer est bleue mais le ciel pas toujours !
def introspectionModel(classif, data, Y, filename=None):
# introspection: qu'est ce qui est associé à chaque classe de données
plt.figure()
plt.imshow(classif.coef_, interpolation='nearest')
localLabs = [data['labMeaning'][y-1] for y in np.unique(Y)]
plt.yticks(range(len(localLabs)),localLabs)
rgb = [color+' '+strength for color in ['red', 'green', 'blue'] for strength in ['(low)', '(med)', '(high)']]
plt.xticks(np.arange(len(rgb))-0.5,rgb, rotation=45)
plt.vlines(np.array([2.5,5.5]), -1, len(localLabs),linestyles='dashed')
if filename != None:
plt.savefig(filename)
Introspection d'un modèle Naive Bayes:
Etape 3: évaluation en validation croisée
L'avantage d'utiliser une bibliothèque de machine learning comme sklearn est de bénéficier des outils annexes comme les stratégies d'évaluation. La séparation en apprentissage/test et même la validation croisée deviennent très faciles à utiliser:
# définition d'un objet Validation Croisée:
cvk = cv.StratifiedKFold(Y, n_folds=5) # 5 sous-ensembles de données
classifier = lin.LogisticRegression();
k=0
for train_index, test_index in cvk: # parcours des 5 sous-ensembles
classifier.fit(X[train_index],Y[train_index])
ypredL = classifier.predict(X[train_index])
ypredT = classifier.predict(X[test_index])
print "(RL) iteration ",k," pc good (Learn) ",np.where(ypredL == Y[train_index],1.,0.).sum()/len(train_index)
print "(RL) iteration ",k," pc good (Test) ",np.where(ypredT == Y[test_index],1.,0.).sum()/len(test_index)
k+=1
L'ensemble du protocole ci-dessus vise à vous donner une performance de référence que l'on va tenter de battre en prenant en compte la structure de l'image.
Prise en compte de la structure de l'image
Le but est maintenant de passer à un apprentissage structuré de la forme MRF (Markov Random Field):
- Les états correspondent aux classes présentent dans les images
- Les observations sont multi-variées: ce sont les 9 descripteurs associés à une région de l'image
Construction de la matrice de transition
La matrice de transition complète correspond à {$nbClasses^5$} coefficients: nous ne disposons pas d'assez de données pour envisager un tel modèle. Nous allons donc faire l'hypothèse que les 4 directions sont indépendantes et nous allons construire 4 matrices de transitions {$Ag, Ad, Ah, Ab$} de {$nbClasses^2$} paramètres.
Donner le code de la méthode A = MRFLearn(data). A contiendra 4 couches correspondant à {$Ag, Ad, Ah, Ab$}.
Notes:
- Dans
data['graph'][i]se trouvent les connexions des 16x16=256 noeuds de l'imagei - Les connexions de chaque noeud sont données dans l'ordre:
[g,d,h,b] - Les connexions à -1 signifient que le noeud du graphe est au bord de l'image : ça vaut le coup de créer une classe bord de l'image pour améliorer les résultats.
Classification par Gibbs Sampling
La résolution de Viterbi dans un contexte cyclique est délicate, nous optons pour une alternative plus simple: une approche par échantillonnage. Même si ce TME intervient avant le cours dédié, la technique est suffisamment intuitive pour être comprise.
Nous cherchons à tirer des possibilités crédibles d'étiquetage de l'image, c'est à dire des échantillons {$s_1,\ldots,s_{256}$} tirés selon {$P(S_1, \ldots, S_{256}| X_1, \ldots, X_{256})$} en notant {$s$} les états et {$x$} les observations. Ensuite, nous étudierons quelles sont les étiquettes les plus probables pour chacun des états. Le tirage selon cette loi n'est pas aisé, mais l'algorithme de Gibbs va nous aider: comme le montre l'illustration ci-dessous, cet algorithme permet de tirer des des échantillon selon une loi jointe complexe en tirant itérativement une par une des variables aléatoire selon la loi conditionnelle.
Nous faisons toujours l'hypothèse d'un champs d'ordre 1: l'étiquette d'un noeud dépend des noeuds voisins et de l'observation. Nous ferons l'hypothèse que : {$$ p(S_i = s | S_1, \ldots, S_{i-1}, S_{i+1}, S_{256}, X_1, \ldots, X_{256} = p(S_i = s | vois(i), X_i )$$} Comme il n'est toujours pas aisé de tirer selon cette loi de probabilité, nous jouons sur la multiplication des tirages et nous tirons en 2 temps:
- choisir {$Y\in\{X_i, vois_{gauche}, vois_{droit}, vois_{haut}, vois_{bas} \} $}
- tirer {$S_i$} selon {$p(S_i = s | Y )$}
- L'algorithme demande un état initial: nous partirons des étiquettes données par la regression logistique
- Une fois que vous avez effectué suffisamment de tirages, vous pouvez déterminer la classe majoritaire associée à chaque segment. Le code suivant peut être utile:
# labeling = nbIterationsGibbs x nbNoeuds
labelEnd = np.array([Counter(labeling[:,j]).most_common(1)[0][0] for j in range(labeling.shape[1])])
- Afficher le résultat, calculer les segments pour lesquels la décision est ambiguë
- Optionnellement, on peut étudier l'introduction d'une période de brun: c'est à dire ne pas prendre en compte les premiers échantillons générés en considérant qu'ils ne sont pas pertinents.
- Optionnellement, calculer les performances en TOP3, c'est à dire compter le pourcentage de fois ou la bonne réponse se trouve parmi les 3 premières réponses.
Alternative plus simple pour la description du contexte :
Une autre stratégie toute simple permet d'améliorer les résultats: la prise en compte de la position du segment dans l'image.
Les tests sont très faciles à effectuer: il suffit d'ajouter les coordonnées des segments (présentes dans data['coord']) aux caractéristiques initiales des segments.
Les coordonnées sont décrites de manières redondantes dans la structure de données fournie:
r = 0 # region dans l'image
coord_image_i = data['coord'][i] # 4 coordonnées x 256 regions
# coord_image_i[r] = coordonnées en x (en pc de la largeur de l'image), 1 - coord_x, coordonnées en y (en pc de la hauteur de l'image), 1 - coord_y
En utilisant la commande np.hstack, on peut fusionner horizontalement 2 matrices: il suffit ensuite de construire un classifieur sur cette représentation enrichie.
Etudier la variation de performances entre les deux approches.
Partie facultative
La partie facultative consiste à implémenter l'algorithme de Viterbi pour les arbres binaires que nous avons vu en TD. Dans un premier temps, plusieurs fonctions de calcul de l'arbre vous sont données et expliquées. Dans un deuxième temps, vous devez implémenter l'algorithme étudié en TD.
L'ensemble de l'approche est décrite dans la référence: lien
Construction d'un arbre binaire (aléatoire) dans l'image
Par défaut, l'image segmentée est plutôt un graphe. La transformation en arbre binaire que nous proposons est arbitraire: il s'agit de lier chaque région soit à celle du dessus, soit à celle de gauche. Même si l'arbre est aléatoire, nous espérons que les transitions vont permettre de capter une sémantique (e.g. le ciel est généralement au dessus des autres concepts...).
Dans une représentation en arbre (telle que représentée ci-dessous), chaque noeud n'a qu'un seul parent: la génération aléatoire consiste donc pour chaque noeud à déterminer si le parent est au-dessus ou à gauche. Le code d'implémentation est donc réduit.
A partir d'un graphe contenant une liste de quadruplets g,d,h,b de liens:
def computeTree(graphe): # randomTree
return np.array([-1 if ((np.random.rand()<0.5 or h == -1) and not g == -1) else 1 for g,d,h,b in graphe])
 | => |  |
Le code pour dessiner les graphes et les arbres est disponible lien
Les outils dont nous avons besoin pour manipuler les arbres sont les suivants:
- récupération des feuilles (pour vérification du code, inutile dans l'algorithme)
- détermination des noeuds parents pour chaque noeud
- détermination des noeuds enfants (gauche et droit)
- parcours en largeur de l'arbre = mise en ordre de l'indice des noeuds pour qu'un enfant soit toujours vu après ses parents (ou l'inverse si on parcours la liste à l'envers)
def computeTreeFather(graphe, tree):
return np.array([g if tr == -1 else h for (g,d,h,b),tr in zip(graphe,tree)])
# retourne une matrice de taille (2,nbnoeuds) donnant la référence des 2 noeuds fils pour chaque parent
def computeChild(graphe, tree):
# pour chaque noeud, le successeur en bas et à droite
return np.array([[b if (b != -1 and tree[b] == 1) else -1, d if (d != -1 and tree[d] == -1) else -1] for (g,d,h,b) in graphe])
# parcours en largeur de l'arbre pour avoir les noeuds dans le bon ordre lors du calcul de Viterbi
def orderedTree(childs):
root = 0 # definition
order = [root]
toTreat = childs[root].tolist()
while len(toTreat) > 0:
i = toTreat.pop(0)
if i==-1:
continue
order.append(i)
toTreat = toTreat + childs[i].tolist()
return order
def leafs(graphe,tree):
return np.setdiff1d(np.arange(len(graphe)), computeTreeFather(graphe, tree))