Objectif

Nous allons pratiquer la programmation dynamique à travers différent un algorithme d’analyse des formes en traitement d’images.

On cherche à identifier les structures foncées de l’image dont dont la taille est inférieure à N (i.e. les structures qui ne tiennent pas dans un carré de N × N.

Pour cela, on effectue une fermeture par élément structurant carré de taille N × N et on calcule le résidu par rapport à l’image d’origine. Par exemple, le résidu après fermeture par un carré de taille N × N donne:

Si on zoom, on se rend compte alors que uniquement les petits éléments noirs (ceux qui ne contiennent pas un carré de 10 × 10 sont présents. Les éléments sombres plus grands n’apparaissent pas.

Une fermeture est une opération morphologique coomposée d’une dilatation suivie d’une érosion. Une dilation affecte pour chaque pixel d’entrée (x,y) la valeur du maximum dans une fenêtre autour de (x,y) (par dualité, l’érosion affecte le minimum).

Exemple de dilatation d’une image par un carré 3 × 3

Travail

Nous allons implémenter différentes versions de ces algorithmes dont la plus efficace utilisera la programmation dynamique.

Les fonctions utisent l’API image suivante:

struct image2d
{
    uint8_t* buffer; // Pointer to the data buffer
    int width;       // Width of the image
    int height;      // Height of the image
    int stride;      // Number of bytes between two lines
};

// Charge depuis un fichier png
image2d* load_image(const char* filename);

// Sauvegarde dans un fichier png
void save_image(const char* filename, image2d* ima);

// Crée une image niveau de gris de taille width * height 
image2d* create_image(int width, int height);

// Détruit une image allouée avec `create_image`
void     destroy_image(image2d* image);

Récupérer l’archive algo-tp01.tar.gz puis builder là:

$ tar xf algo-tp01.tar.gz && cd algo-tp01
$ mkdir build && cd build
$ cmake ..

Cela crée deux exécutables: * process qui exécute l’algorithme de fermeture par un carré de taille K sur une image et sauvegarde le résultat (pour l’instant le résultat est tout noir). On l’utilise ainsi:

  ./process input.png K output.png

• bench qui permettra de mesurer la performance des algorithmes. Cela affiche:

    ----------------------------------------------------
    Benchmark             Time           CPU Iterations
    ----------------------------------------------------
    BMBase/V1/3           3 ms          3 ms        242
    BMBase/V1/11          3 ms          3 ms        241
    BMBase/V1/51          3 ms          3 ms        243
    BMBase/V2/3           3 ms          3 ms        241
    BMBase/V2/11          3 ms          3 ms        240
    BMBase/V2/51          3 ms          3 ms        239
    BMBase/V3/3           3 ms          3 ms        241
    BMBase/V3/11          3 ms          3 ms        241
    BMBase/V3/51          3 ms          3 ms        240

Implémentation naïve

Implémenter naïvement l’algorithme de fermeture par un carré de taille K × K. La fonction est à coder dans closing.c de prototype: image2d* closing_1(const image2d* input, int K).

L’algorithme naïf de dilatation fait simplement une double boucle:

Pour chaque pixel (x,y):
    m = -∞
    Pour chaque voisin (p,q) autour de (x,y) dans une fenêtre de taille *K*
        Si (p,q) n'est pas en dehors de l'image et f(p,q) > m:
            m = input(p,q)
    output(x,y) = m

L’érosion s’écrit de manière similaire et la fermeture est simplement la succession des deux opérations.

• Question quelle est la compléxité de cet algorithme en fonction de la taille du carré K et du nombre de pixels N ?

Implémentation moins naïve

Le carré est en fait un élément structurant séparable, ce qui veut dire qu’éroder avec un carré de taille K × K revient à éroder par un segment horizontal de taille K, puis éroder par un segment vertical de même taille. Une érosion 2D devient donc une succession de deux érosions 1D. Cette propriété est valable aussi pour la dilatation. Implémenter cette version dans la fonction closing_2.

• Question quelle est la compléxité de cet algorithme en fonction de la taille du carré K et du nombre de pixels N ? Vérifier le empiriquement avec les benchmarks.

Implémentation dynamique 1

L’implémentation précédente se résume donc à effectuer des érosions/dilatations 1D. Autrement dit, cela revient à trouver le maximum ou le minimum sur une fenêtre glissante. On remarquera que l’on se décale à chaque fois de 1 pixel. Lorsqu’on calcule le maximum sur l’ensemble [-K+x, -(K-1)+x, -(K-2)+x, ..., -1+x, x, x+1, ..., x+(K-1), x+K], on a juste ajouté l’élément x+K et enlevé l’élément -(K+1)+x. On peut donc mettre à jour dynamiquement cet ensemble à chaque itération en maintenant toujours le max.

1. Sachant que les valeurs sont comprises entre 0 et 255, proposer une structure de données permettant de représenter l’ensemble des valeurs contenues dans la fenêtre avec un ajout/suppression efficace et un accès efficace au maximum de l’ensemble.
2. Implémenter l’algorithme
3. Quelle est la compléxité de cet algorithme. La vérifier empiriquement avec les benchmarks.

Implémentation dynamique 2

Il est possible d’implémenter la recherche du maximum sur fenêtre glissante plus efficacement encore. Simplifions dans un premier temps le problème. Soit A[i] un vecteur d’entier, on cherche à calculer M[i] = Max A[i:i+K] (ie le max sur une fenêtre de K éléments). On découpe le vecteur A en partie de K élements, on pose \[ a = \lfloor[ a \rfloor] \] et on définit les vecteurs L[i] = Max A[a.K:i] (left) et R[i] = Max A[i : (a+1).K] (right) qui sont les maxima locaux à partir des bordures gauches et droites.

1. Exprimer M[i] en fonction de L[i] et R[i]. En déduire la complexité en terme de comparaisons du calcul de M[i].
2. Faire le lien entre M[i] et notre problème originel qui consiste à calculer le maximum (et le minumum) sur A[i-K : i+K+1] (i.e. comment ‘recentrer’ la fenêtre).
3. Implementer l’algorithme et mesurer ses performances.