TP4: Assemblage de photos

- Création de mosaïques d'images manuelle et automatique -

par Olivier Duguay

GIF-4105/7105 Photographie Algorithmique, Hiver 2015

1-Objectifs du projet 2-Description de l'algorithme d'appariement manuel 3-Description de l'algorithme d'appariement automatique 4-Résultats

1-Objectifs du projet

Dans le cadre de ce quatrième projet du cours de photographie algorithmique, il nous faut implémenter une application qui permet de créer un panorama d'images à l'aide de prises de vue se chevauchant. Pour effectuer l'appariement d'image, il faudra d'abord déterminer des points d'intérêts commun entre les images. Pour ce faire, il faudra utiliser une technique de sélection manuelle de points et une autre qui sélectionnera des points d'intérêts automatiquement.

Pour plus de détails, voyez l'énoncé complet ICI.

2-Description de l'algorithme d'appariement manuel

Étape 1 - Appariement de caractérisque

Tout d'abord, il faut sélectionner manuellement des points d'intérêts entre 2 images. Habituellement, on utilisera l'image du centre comme image de départ et on s'assurera que les points d'intérêts sélectionnés dans les autres images (de part et d'autre du panorama) sont présents dans cette photo. Cependant, il est également possible de sélectionner des points d'intérêts que pour les paires d'images adjacentes et de faire le calcul de l'homographie «en cascade». Il est important de bien déterminer les pixels (point x,y) qui correspondent dans chaques images puisqu'une simple erreur peu engendrer une grande erreur qui en résulterait d'une déformation peu réaliste de l'image vers le panorama.

Il est important de sélectionner certains nombres de points afin d'obtenir de meilleurs résultats comme nous le verrons dans une prochaine étape. Pour ce faire, une fonction à été créer afin de faire la sélection de points un-à-un sur deux images et de sauvegarder la sélection sur un fichier texte au nom des deux images combinées. Il faut d'abord savoir combien de points on compte utiliser.

saveAndDisplayFeatures(imgfilename1, imgfilename2, n)

Étape 2 - Calculer l'homographie

Puis, une fois les points d'intérêts déterminées, il faut déterminer la matrice d'homographie d'images «H» qui permettra la déformation d'une image vers l'autre. Pour ce faire, il faut créer un système de «n» équations linéaires de la forme Ah = b où h correspondera à un vecteur des 8 inconnus de la matrice H. Comme la matrice de déformation est de type projective, c'est la raison pour laquelle on se retrouve avec 8 inconnus (car le 9e élément de la matrice pourra être simplifier et être mis à 1 puisqu'il correspond à un facteur d'échelle). Pour résoudre de le système d'équations, il faudra donc au minimum 4 points d'intérêts.

Cependant, comme il a déjà été mentionné, il faut sélectionner le plus de points possible, le plus précisément possible entre les images et de s'assurer que ceux-ci ne sont pas tous trop près spatialement. Des points tous sur la même ligne de l'image ou un amassi de points au même endroit n'est pas tolérable et de prendre un minimum de points également. Dans la figure suivante, on voit le résultat que l'on obtient pour une sélection de quelques points le long d'une clotûre en bordure d'image.

Mauvais choix de points Déformation non-désirée Bon choix de points Déformation désirée
Mauvais choix de points Déformation non-désirée Bon choix de points Déformation désirée


Pas très satisfaisant si on s'en tient au minimum !

Il faut donc chercher un homographie en surestimant celle-ci à l'aide de plusieurs points qui seront utilisées pour résoudre plusieurs systèmes d'équations et de résoudre ceux-ci en utilisant la méthode des moindres carrée. Dans Matlab, il suffit donc de créer la matrice d'équations linéaires en utilisant les points d'intérêts de l'image source que l'on utilise sous forme homogène. La matrice de sortie correspond aux points de l'image de destination et il faut donc trouver les inconnus qui permettent cette transformation de façon fidèle et juste. On utilisera l'opérateur mldivide «\» et il faut également savoir qu'il est plus simple de travailler sous forme vectorielle pour définir les équations linéaires et d'utiliser celle-ci sous forme d'une matrice par la suite.

Il nous à donc fallu créer la fonction ;

H = computeH(im1_pts,im2_pts)

qui permet de déterminer l'homographie et donc, pour chaque points, on il faut reconstituer le système d'équations ce qui mènera à la matrice suivante:

 
0 0 0 -P1x -P1y -1 P1x * P2x P1y * P2y
P1x P1y 1 0 0 0 -P1x * P2x -P1y * P2y
... ... ... ... ... ... ... ...
 
h =
 
-P2y
P2x
...
 


Comme il à été mentionné précédemment, il est plus simple de toujours calculer l'homographie pour un même plan (image de référence) et de déformer les autres images par rapport à celui-ci. Cependant, une meilleure technique serait de calculer l'homographie pour chaque paires d'images adjacentes dans le panorama et d'établir l'homographie en cascade. Pour ce faire, par exemple si on 5 images et que la 3e est celle du centre, normalement, on choisirait des points d'intérêts pour les paires d'images suivantes (1-3, 2-3, 4-3, 5-3) pour ainsi toujours trouver des homographies en fonction de l'image centrale et donc de déformer vers l'image 3. Cependant il est possible d'établir les correspondances suivante (1-2, 2-3 ,4-3, 5-4) et de déterminer les homographies des images non-adjacentes à l'image du centre par le biais des autres homographies. Ainsi, par exemple, pour l'image 1, on calculerait l'homographie en cascade de la manière suivante H12 * H23 = H13.

Étape 3 - Déformation projective

Résultat de déformation d'image
Résultat de déformation d'image
Une fois toutes les matrices d'homographies «H» de calculées pour les paires d'images, on peut déterminer la déformation de l'image de destination vers sa source. Dans le cadre du projet, nous avons choisi d'utiliser la transformation inverse pour calculer la déformation car cela semblait donner de meilleurs résultats (ce qui pourrait être sujet à l'argumentation). Il ne faut pas oublier de calculer H dans la bonne direction (source->destination VS. destination->source) si on utilise la transformation inverse.

La déformation peut aisémemt se calculer à l'aide de la fonction Matlab interp2 qui sera utilisé sur tout les canaux (RGB) de l'image à l'aide de la fonction ;

imwarped = warpImage(allImages, allMatrixH ,X,Y)

De plus les coordonnées d'interpolation seront donc établie à l'aide de l'homographie en utilisant notre fonction ;

[pts_x2, pts_y2] = useH(matrixH, pts_x, pts_y)

Il est à noter, que pour établir la déformation et surtout la nouvelle image, il faut savoir le format (taille) de l'image résultante qui correspond à la mosaïque d'images. Pour déterminer la format du panorama final, la fonction suivante à été codé ;

[X, Y] = findPanoSize(allImages, allMatrixH)

Cette fonction permet, à l'aide du lot d'images et des homographies, de déterminer l'emplacement des quatres coins de l'image résultante et par le fait même, de connaître sa dimension. En effet, en effectuant un «forward-warping» que sur les coins des images et il est donc possible de déterminer les coins minimaux et maximaux en X et Y quicorrespondent aux quatres coins du panorama.

Image destination (référence) Image source (à déformer)
Image destination (référence) Image source (à déformer)

Étape 4 - Fusionner les images en une mosaïque

Pour fusionner les images entre elles, il ne suffit pas toujours de les superposer car on remarque trop l'effet de «collage» de photos. Une manière de contrer cet effet, est d'utiliser une moyenne pondérée sur chacuns des pixels qui se chevauchent. Une fonction à donc été programmé dans le but d'accomplir cette tâche :

[panoramaFinal] = createPano(imgWarped)

En y insérant les images déformées, cette fonction additionne les pixels qui se chevaucheront et détermine le nombre de chevauchement par pixels afin d'obtenir la moyenne pondérée des images. Il suffit de superposer les images par la suite.

3-Description de l'algorithme d'appariement automatique

Étape 1 - Détection de caractérisques de coins dans une image

Points d'Harris
Points d'Harris

Dans un premier temps, pour trouver les coins caractérisques d'une image, on utilisera l'algorithme d'Harris qui est un classique en matière de détection de coins. Il es possible de faire la détection à plusieurs échelle bien que dans le but de ce projet ce n'est pas nécessaire que nous utiliseront la fonction existante suivante:

function [newy,newx,newv] = harris(imrgb, nbPtsHarris, robust)

Cette fonction a été modifiée légèrement de son originale puisqu'elle comprend maintenant une section qui permet de choisir le nombre de points à retenir pour faire l'homographie. De plus, elle implémente maintenant la répression adaptive non-maximale qui permet de faire la sélection le plus judicieusement possible, c'est à dire de manière uniforme et pour des points très distincts. Cette technique est un peu complexe mais fonctionne très bien. On tri d'abord les coins/points selon leurs «forces» et le maximum global sera notre points d'entrée. En formant un certain rayon autour de ce point qu'on élargit peu à peu, on obtient d'autres points maximaux qui doivent respecter le critère de robustesse (variable d'entrée de la fonction). On trouve donc le rayon minimal en utilisant la fonction pdist2 et ainsi obtenir le nombre de points désirés avec des «forces» suffisament élevés.

Dans l'image à droite, on peut voir les milliers de points en ROUGE qui correspondent aux points de Harris créer par la fonction, puis en BLEUE, les points sélectionnés en utilisant la technique de répression adaptive non-minimale.

Étape 2 - Extraire un descripteur pour chaque point caractéristique

Une fois les séries de points déterminés dans chacunes des images du panorama, il est important dans déterminer un descripteur. Les descripteurs sont des groupements de pixels de 8x8 pris autour des points caractéristique d'Harris qui décrivent en quelque sorte, l'identité du pixel dans le but de faire l'appariement. Bien qu'il est favorable d'utiliser l'invariance en rotation de ces descripteurs puisque la prise de photos n'est pas toujours parfaitement «droit» et sans rotation, cette technique n'a pas été utilisé. Cependant dans le but de rendre plus rapide les calculs, on utilise ces descripteurs sur une plus grande fenêtre de dimension 40x40 soit 5x plus grand et légèrement flou. La fonction suivante permet de faire cette étape :

[descripteurs] = getDescript(image, ptsHarris)

Comme indiquer cette fonction applique d'abord un flou sur l'image à l'aide de imfilter, puis on transforme un domaine spatiale de afin d'obtenir 5 pixels pour un (8x8 - 40x40) à l'aide de meshgrid. Finalement, il suffit d'interpoler les points en utilisant interp2 et surtout, de ne pas oublier de normaliser les points avec le gain et le biais (moyenne (mean) et variance sqrt(var(...)) des descripteurs).

Étape 3 - Apparier ces descripteurs de caractéristiques entre les deux images

À l'aide de la fonction suivante il est possible d'apparier les différents descripteurs provenant des différents lots de points de Harris déterminés auparavant.

[ptsDest, ptsSource] = matchDescript(descript1,descript2,ptsHarris1,ptsHarris2)

Cette fonction utilise la fonction de distance (pdist2) pour déterminer la distance entre les paires de descripteurs. Pour l'appariement il suffit d'utiliser le plus proche voisin ou d'utiliser la technique de Lowe qui calcule le rapport entre les 2 plus proches voisins du descripteur. Voir le fichier suivant pour plus de détails.

Étape 4 - Utiliser une méthode robuste (RANSAC) pour calculer l'homographie

Finalement, après avoir apparier les descripteurs caractérisques les uns avec les autres, il faut estimer l'homographie. Cette tâche peut être très longue si on utilise tous les points de Harris, ou même tous les points d'intérêts sélectionnées. Pour ce faire, il existe l'algorithme de RANSAC (RANdom SAmple Consensus) qui n'est pas très difficile à implémenter. Par itération (plus le nombre sera grand, plus l'algorithme sera robuste), on prend aléatoirement 4 points pour en déterminer l'homographie en plus de déterminer le nombre de points qui satisfait l'appariement. La tolérance d'erreur des pixels peut être modifié dans la fonction suivante qui effectue l'algorithme de RANSAC et qui retourne la meilleure homographie d'image :

[matrixH] = findRANSAC(ptsDest, ptsSource)

4-Résultats

Dans le but d'alléger le site web, les panoramas ont été simplifiés

Stade du Rouge et Or
Stade du Rouge et Or


Grand Axe de l'université
Grand Axe de l'université


Linges
Linges


Comme on le remarque sur cette dernière photo, lorsque la lumière varie beaucoup, il est difficile pour l'algorithme de bien détecter les points qui sont semblable et aussi le fait qu'il manque un peu de coins dans l'image n'aide pas beaucoup à avoir un beau panorama final.

Dominos
Dominos


Comme on le remarque sur la photo précédente, en présence de forts points d'intérêts/coins et le fait que ceux-ci soient en bons nombres, on trouve une correspondance quasi-parfaite.