Python : détection de particules de forme circulaire sur une image par technique du watershed

Temps de lecture : 8 minutes

Voici ce que je te propose de réaliser dans cet article :

Détecter et classifier (grandes ou petites) des particules de forme circulaire sur une image !

[Code complet (Jupyter Notebook) et ressources sur Github]

Contexte réel

L’image à traiter est une image type de mes expériences en doctorat sur les dynamiques de mélange de particules solides. Deux populations de particules cylindriques creuses (grandes et petites) sont mélangées dans un domaine circulaire. Des images du milieu sont prises régulièrement pour étudier le changement de positions des particules. Une image type à un instant t sortie par la caméra ressemble à ceci, c’est tout ce qu’on a comme information d’entrée.

Image brute

Il faut donc qu’on puisse retrouver toutes les particules dans un premier temps ! Rentrons tout de suite dans le vif de sujet.

Prétraitement (la base de tout bon traitement)
Normalisation

Cette image est codée sur du 64 bits i.e. les pixels ont des valeurs d’intensité entre 0 et 65536 (2^{16}). Tout bon data scientist te dira qu’il est préférable de travailler avec des datas normalisées. Dans notre cas, nous allons employer une normalisation min-max de sorte à avoir de nouvelles valeurs qui vont de 0 à 1. Le principe ? Si img représente l’ensemble des intensités de pixel de l’image, la valeur normalisée norm pour chaque pixel sera :

norm = \frac{img - min(img)}{max(img) - min(img)}

Comme on aime être autonome, créons nous-même notre petite fonction pour faire ça :

# Fonction normalisation de l'intensité d'une image entre 0 et 1
def normalise(img):
    norm = (img - np.min(img)) / (np.max(img) - np.min(img))
    return norm

# Normalisation
image_normalise = normalise(image)

Je rappelle que c’est juste une mise à l’échelle pour fixer les limites min et max de nos intensités d’image, cela ne change rien à l’aspect. Tu peux vérifier en affichant l’image obtenue après la normalisation.

Masque et rognage

Deuxième étape, tu as dû peut-être remarquer que pour la détection, seule la partie centrale de l’image avec les particules nous intéresse ici : nous allons donc nous débarrasser des « parties inutiles ». Comment ? je te dis tout. On va d’abord commencer à récupérer le centre et le rayon du domaine central en cliquant sur 3 points sur le pourtour, puis on va créer un masque pour cacher la partie externe au domaine central (j’ai écrit deux mini-articles sur chacun de ces sujets si tu veux mieux comprendre).

# Fonction génération de masque circulaire à partir d'une image, d'un centre et d'un rayon
def genere_mask (image, centre, rayon):
    shape = np.shape(image)
    basex = np.arange(0,shape[0]); basey = np.arange(0, shape[1])
    [i0,j0] = centre
    mask = (basex[:,np.newaxis] - i0)**2 + (basey[np.newaxis,:]-j0)**2 < rayon**2
    return mask

# Création d'un masque sur les régions non-utiles (extérieur du domaine circulaire central)
mask = genere_mask(image_normalise,centre,rmax)
image_normalise[~mask] = 1

Petite touche finale : on découpe l’image pour garder uniquement le carré circonscrit au domaine central.

# Fonction détourage carré d'une image 
def crop_image (image,centre, rayon):
    [i0,j0] = centre
    image_cropped = image[int(i0 - rayon): int(i0 + rayon), int(j0 - rayon): int(j0 + rayon)]
    return image_cropped

# Détourage de l'image 
image_detoure = crop_image(image_normalise, centre,rmax)
Image rognée

Les plus attentifs auront remarqué que vu l’arrière-plan « clean » et uniforme de l’image brute, on n’avait pas vraiment besoin de créer un masque mais je voulais te le montrer quand même, il peut arriver que tu aies des images brutes moins « clémentes » et dans ces cas, tu n’échapperas pas au masque ! Autant savoir comment l’utiliser.

Amélioration de contraste

Bon, là, l’image rognée est plutôt relativement bien contrastée (dis-moi si tu veux savoir comment je l’ai obtenue expérimentalement), mais si on est un peu tatillon et qu’on regarde de très près, on peut constater qu’il y a un peu plus de lumière dans la zone centrale de l’image que dans ses parties périphériques, c’est minime mais autant essayer d’uniformiser ça tout de suite. Python propose une fonction adjust_gamma  dans la librairie skimage.exposure. C’est une correction gamma de l’intensité de l’image. Sans rentrer dans les détails, contentons-nous juste de l’appliquer ici pour uniformiser un peu les intensités de l’image.

# Amélioration du contraste, ajustement gamma
image_gamma = exposure.adjust_gamma(image_detoure,gamma)

Eh, avoue, il y a une petite différence quand-même, les tons sont plus ressortis et c’est un peu plus uniforme. Allez, passons à l’étape suivante 🙂!

Binarisation

Jusque-là, on a apprêté l’image pour la rendre plus « jolie ». Mais ne perdons pas de vue l’objectif initial : détecter les particules dessus et leur position. L’image obtenue après l’étape de la correction gamma est une image en niveaux de gris. Maintenant, on essayera de dire à Python : « les zones qui représentent une particule, affecte-leur une valeur de 0 et celles qui ne le sont pas, va pour la valeur 1  » : on dit qu’on binarise ou qu’on seuille l’image. Tu l’as certainement deviné, les particules sont noires (sauf la tige) et le reste de l’image est blanc. On peut donc :

  • Soit se fixer un seuil global S et on considérera que si dans un pixel, l’intensité est inférieure à S (tendance vers le noir), c’est une particule, on affecte donc 0 et sinon on affecte 1 (intensité supériure à S) : c’est un seuillage global
  • Soit parcourir l’image localement et dans chaque petite région, décider de la valeur de seuil qui sépare le mieux les parties blanches des parties noires : c’est un seuillage local et ça a l’air bien mieux et plus précis! La fonction python qui fait ça, c’est skimage.filters.threshold-local : 
# Matrice de seuils locaux
seuil_local = threshold_local(image_gamma, block_size_threshold)

# Binarisation de l'image
image_binaire = image_gamma > seuil_local

Une façon de vérifier si ta binarisation ne te fait pas perdre en information, c’est de vérifier que l’image obtenue après n’est pas trop différente de l’image initiale, dans notre cas, on n’est plutôt pas mal. Ne regarde pas la tige ou l’arrière-plan, ce qui nous intéresse c’est les particules, tu peux zoomer chez toi et tu verras que la binarisation est fidèle à l’image originale.

Technique du watershed ou ligne de partage des eaux

C’est une technique de segmentation assez utilisée en traitement d’images traditionnel. Comme le nom l’indique, imagine qu’on verse de l’eau dans le domaine et que les parties en noir agissent comme des murs. Dans un tel scénario, l’eau va stagner entre les murs et former des bassins versants (watershed en anglais) ! Génial… mais attends, pourquoi est-ce que c’est génial ? Nos particules, elles, sont des surfaces fermées n’est-ce-pas ? Cela veut dire qu’on aura de l’eau stagnante (bassin) dans l’espace intérieur des particules et vu qu’en plus, les particules sont à géométrie circulaire, un tel bassin sera lui aussi circulaire ! On pourra donc désormais retrouver nos particules en repérant tous les bassins ayant la forme d’un disque ! tu comprends maintenant pourquoi je disais que c’était génial ? Il faut quelques lignes de code pour retrouver les bassins (je détaille ça dans un autre article si tu veux connaître tous les ingrédients) 🙂 !

# Technique du watershed pour créer des bassins, voir mes autres articles pour l'explication pas-à-pas
distance = ndi.distance_transform_edt(image_binaire)
l_min = morphology.local_minima(-distance)
markers, _ = ndi.label(l_min)
labels = watershed(-distance, markers = markers, mask=image_binaire)
labels, _ = ndi.label(labels)
Exemple de bassins obtenus par technique de watershed

En gros, si on regarde une partie quelconque de l’image après cette étape, on a tous les bassins qui sont colorés : ceux correspondant à l’intérieur des particules ont une forme circulaire régulière et les autres bassins correspondant aux espaces interstitiels par exemple ont des formes très irrégulières. L’avantage, c’est que Python nous donne accès aux propriétés de ces bassins, leur aire, périmètre, grand axe, petit axe…etc. Nous allons utiliser certaines de ces propriétés pour retrouver nos particules et ensuite les classifier en petites ou grandes particules.

# Récupération des propriétés des régions labellisées créées par le watershed
regions = regionprops(labels)
Détection et classification

Ici, globalement, ce qu’on fait, c’est qu’on essaie de formuler mathématiquement ce qu’est un cercle/disque. Moi j’ai choisi les critères ci-dessous, mais tu es libre de définir tes propres critères si tu en vois d’autres 😉 :

  • Convexité de la forme : les formes irrégulières ici sont plutôt concaves i.e. « des formes avec des pics ». Un cercle n’a pas de pic : il est convexe donc si on calcule le ratio entre l’aire du bassin circulaire et l’aire de son équivalent convexe, on doit être proche de 1 pour les particules et de 0 pour les formes irrégulières. Python propose déjà un paramètre appelé solidity qui fait ça. Bon bien sûr, le monde n’étant toujours pas parfait (oui, c’est vrai), on n’obtient pas exactement 1 pour nos particules, on peut néanmoins fixer un seuil à 0.95
  • Différence de longueur d’axes : on peut aussi remarquer que les formes irrégulières sont souvent très allongées i.e. leur grand axe est grand devant le petit axe. On peut donc mettre un critère sur la différence entre les deux axes. Python donne déjà accès aux valeurs de ces longueurs. Ici on fixe un seuil à 4.
  • Classification en grandes ou petites particules : après avoir éliminé les bassins qui ne sont pas des particules, il faut désormais classifier les bassins restants : facile ! tu as peut-être deviné, on peut utiliser l’aire du bassin. Mais pour minimiser l’effet des contours irréguliers, on peut prendre plutôt l’aire du rectangle circonscrit au bassin. On fixera des seuils min et max pour les petites comme pour les grandes particules. Encore, ici je précise, c’est du freestyle, tu décides de ce qui marche le mieux pour tes images.
# Détection des particules
labels_number = np.arange(1,np.max(labels)+1) # Récupération des numéros de toutes les régions labellisées
# la liste "taille" va stocker les tailles pour chaque particule, petite : 1, grande : 2
taille = []; 
# la liste "x" stocke la coordonnée x des particules
x = []; 
# la liste "y" stocke la coordonnée y des particules
y = []; 
# la liste 's' stocke la solidité des particules
s = [];


#Première sélection basée sur les attributs géométriques
for label in labels_number :
    region = [region for region in regions if region.label == label][0]

    solidity = region.solidity # solidité de la région 
    centroid = region.centroid # coordonnées du centre de la région
    b = region.bbox # coordonnées des coins du rectangle circonscrit à la région
    b_area = region.bbox_area # aire du rectangle circonscrit à la région
    dxy = abs(abs(b[2] - b[0]) - abs(b[3] - b[1])) # différence entre le grand et le petit axe de la région


    if(solidity >= solidity_min and dxy <dxymax ): # Elimination des formes non circulaires
        if (est_compris_entre(b_area, p_min_area, p_max_area)): # Récupération petites particules
            x.append(centroid[1]); y.append(centroid[0])
            s.append(solidity)
            taille.append(1) # 1 = petite particule
        else :
            if (est_compris_entre(b_area, g_min_area, g_max_area)): # Récupération grosses particules
                x.append(centroid[1]); y.append(centroid[0])
                s.append(solidity)
                taille.append(2) # 2 = grande particule
Touche finale : visualisation graphique

Nous sommes presque à la fin. Détection ? fait ! Classification ? fait ! A la fin de cette étape, on a récupéré les coordonnés des centres des particules et stocker aussi leur taille (on a affecté 1 pour une petite particule par exemple et 2 pour une grande particule). Maintenant il ne reste qu’à représenter tout ça sur une jolie figure et à mettre les bonnes couleurs 😊.

Image finale

PS : Comme dans tous les problèmes réels, il peut arriver qu’on classifie ou élimine à tort certaines particules, on peut utiliser un dernier critère sur les distances entre particules pour essayer de corriger ça. Si ça t’intéresse, je peux écrire un article dessus mais je crois que celui-ci est déjà suffisamment long comme ça.

Bon, maintenant qu’on a les coordonnées des particules, on peut essayer à l’instant t, de les mettre en mouvement puis essayer de retrouver quelle particule s’est déplacée où à l’instant t+1, un peu de tracking, ça te dit ?  

Partager

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *