Tuto: EM algorithm in mixture model¶

Stéphane Derrode, École Centrale de Lyon

Ce tuto montre comment appliquer un modèle de mélange gaussien pour la segmentation d’une image ; il s’agit d’affecter à chacun des pixels de l’image une numéro de classe, parmi [1, …, K]. Nous travaillerons sur une image de cible à deux classes (K=2). Pour être en mesure de mesurer le taux d’erreur, nous bruiterons nous-même l’image au préalable. Ensuite, en faisant mine de ne pas connaître les paramètres du bruit, nous essaierons de les estimer à l’aide de l’algorithme itératif EM (Expectation-Maximization), tel qu’il a été vu en cours. Enfin, sur la base des paramètres estimés, nous appliquerons la décision bayésienne et finirons par comparer l’image originale avec l’image segmentée pour estimer le taux d’erreur de classification.

Le tuto consiste à intégrer progressivement certaines fonctions d’un programme écrit en Python. Il s’agira ensuite d’examiner les courbes de convergence de l’algorithme EM en fonction des paramètres de bruitage de l’image.

Préambule¶

Téléchargez et décompressez le fichier zip disponible en suivant ce lien.

Le répertoire contient

• Le répertoire sources contient des images noires et blanches.
• Le répertoire results où seront stockés les images résultantes de l’exécution des programmes.
• Un programme pour ajouter du bruit à l’image de cible appelé BruitageImage.py.
• Le squelette du programme principal intitulé LabSession_MM_EM.py, ainsi que le fichier func.py contenant des fonctions partiellement codées, à compléter durant le tuto.

Bruitage d’une image¶

• Exécutez le programme BruitageImage.py et observez les images générées par le programme dans le répertoire results.
• Analysez le code source pour en comprendre les étapes. Jouez avec les paramètres du bruit et observez l’impact sur les images générées.
• Question : Dans le mélange, comment sont fixées les probabilités a priori ? Modifier le programme pour en afficher les valeurs.

Segmentation d’images : implémentation des algorithmes¶

Analysez le squelette de programme du fichier LabSession_MM_EM.py. Celui-ci fait appel à des fonctions qui sont toutes décrites dans le fichier func.py. Dans ce dernier fichier, plusieurs fonctions sont vides : à vous de les compléter en suivant les consignes ci-dessous !

Truc :

• Travaillez avec la petite image (cible_64.png), sinon vos essais vont prendre trop de temps !
• Réduisez à 3 ou 4 le nombre d’itérations de EM, le temps de réaliser vos tests. Sachant qu’il faut entre 50 et 100 itérations pour obtenir une bonne estimation des paramètres du modèle.

Travail à effectuer

Il vous suffit ici de comprendre les ajouts de code proposés (à intégrer au fichier func.py) et de faire les tests pour vérifier que vous avez bien compris (référez-vous au slides pour les explications).

1. Examinez la structure du programme. Vous y trouverez 5 étapes :

   1. La lecture des images.
   2. L’initialisation des paramètres du mélange.
   3. L’estimation des paramètres du mélange par l’algorithme EM.
   4. La classification de l’image selon la décision bayésienne.
   5. Le dessin des courbes.

   À chaque itération, les paramètres et l’erreur de classification sont sauvegardés, en vue du dessin des courbes de la dernière phase.

2. Complétez la fonction d’initialisation des paramètres avec l’algorithme suivant

def InitParam(K, imageNoisy):

    meanTab = np.zeros(shape=(K))
    varTab  = np.zeros(shape=(K))
    piTab   = np.zeros(shape=(K))

    # compute the min, max, mean and var of the image
    minB  = int(np.min(imageNoisy))
    maxB  = int(np.max(imageNoisy))
    meanB = np.mean(imageNoisy) # comment: not necessary
    varB  = np.var(imageNoisy)
    # print(minB, maxB, meanB, varB)
    # input('pause')

    # set the initial values for all the parameters
    for k in range(K):
        # Init Pi
        piTab[k] = 1./K
        # Init Gaussian
        varTab[k]  = varB / 2.
        meanTab[k] = minB + (maxB-minB)/(2.*K) + k * (maxB-minB)/K

    return meanTab, varTab, piTab

3. Complétez la fonction qui itère EM avec les formules de re-estimation vues en cours.

def EM_Iter(iteration, imageNoisy, meanTabIter, varTabIter, piTabIter):

    L, C = np.shape(imageNoisy)
    K    = np.size(meanTabIter[0, :])
    N = L*C

    # calculating gamma_n(k)
    gamma = np.zeros(shape=(N, K))
    for l in range(L):
        for c in range (C):
            n = l*C+c
            for k in range(K):
                gamma[n, k] = piTabIter[iteration-1, k] * norm.pdf(imageNoisy[l, c], loc=meanTabIter[iteration-1, k], scale=np.sqrt(varTabIter[iteration-1, k]))
            gamma[n, :] /= np.sum(gamma[n, :])

    # Updating pi a priori proba and Gaussian mean
    for k in range(K):
        temp  = 0.
        temp2 = 0.
        for l in range(L):
            for c in range (C):
                n = l*C+c
                temp  += gamma[n, k]
                temp2 += gamma[n, k] * imageNoisy[l, c]
        piTabIter  [iteration, k] = temp / N
        meanTabIter[iteration, k] = temp2 / temp

    # Updating Gaussian variance
    for k in range(K):
        temp  = 0.
        temp2 = 0.
        for l in range(L):
            for c in range (C):
                n = l*C+c
                temp  += gamma[n, k]
                temp2 += gamma[n, k] * (imageNoisy[l, c] - meanTabIter[iteration, k])**2
        varTabIter[iteration, k] = temp2 / temp

4. Ajoutez la fonction qui calcule la vraisemblance et la valeur de la fonction auxiliaire Q à toutes les itérations de EM.

def likelihood_Q_computing(imageNoisy, meanTab, varTab, piTab):

    L, C = np.shape(imageNoisy)
    K    = np.size(meanTab)

    Likelihood = 0.
    AuxiliaryQ = 0.

    # computing gamma_n(k)
    gamma = np.zeros(shape=(K))
    for l in range(L):
        for c in range(C):

            for k in range(K):
                gamma[k] = piTab[k] * norm.pdf(imageNoisy[l, c], loc=meanTab[k], scale=np.sqrt(varTab[k]))
            sumg = np.sum(gamma)

            temp1 = 0.
            for k in range(K):
                temp1 += np.log(gamma[k]) * gamma[k] / sumg

            AuxiliaryQ += temp1
            Likelihood += np.log(sumg)

    return Likelihood, AuxiliaryQ

---

Voici les allures des courbes que vous devriez obtenir, après 70 itérations de EM sur le mélange définit dans BruitageImage.py :

final results

Mean error rate by class = [0.26617251 0.02335375]
Global mean error rate   = 0.111328125
Pi estimated [0.33412083 0.66587917]
Mean estimated [ 98.85694127 109.32681671]
Std estimated [6.0567873  3.12997051]