TP3: Assemblage de photos

Sélection manuelle de points

La première étape du travail consiste à sélectionner des points de correspondance entre les images. Il suffira pour cela de clique sur des points de correspondance entre deux images. Ce procédé doit être répété pour chaque image. Si le panorama est constitué de 10 images il faudra sélectionner 9 ensembles de points de correspondance. Pour que l'estimation de l’homographie soit bonne il vaut mieux sélectionner plus de 4 points de correspondance.

Calculer l'homographie

A cette étape, on cherche une homographie H qui permette de déformer une image pour qu'elle se superpose avec une autre. Pour cela, on travaille avec deux jeux de coordonnées de n points correspondants P1 et P2. Il est important de souligner qu'il faut que P1 et P2 soient exprimées en coordonnées homogènes, ce sont donc des matrices nx3 dont la première colonne représente les coordonnées x, la deuxième ligne les coordonnées y et la troisième ligne n'est composé que de 1. On cherche donc H telle que P1H = P2. La matrice H a 8 inconnues (9 éléments dont un facteur d'échelle). On pose le problème sous la forme Ah = b avec h un vecteur contenant les 8 inconnues. Pour trouver A et b il suffit de poser les équations qui sont donné par les points de correspondances. On trouve que les matrices seront de la forme:

 
0 0 0 -P1(1,1) -P1(1,2) -1 P1(1,1)P2(1,2) P1(1,2)P2(1,2)
P1(1,1) P1(1,2) 1 0 0 0 -P1(1,1)P2(1,1) -P1(1,2)P2(1,1)
... ... ... ... ... ... ... ...
 
 h =
 
-y2
x2
 

Une façon d'aborder le problème est de chercher h qui minimise Q = ||b - A*h||^2. On a :

Q = (b-A*h)'*(b-A*h) = b'*b -b'*A*h-h'*A'*b+h'*A'*A*h.

Pour trouver H qui minimise Q, on cherche à résoudre dQ/dh=0:

dQ/dh=0 <=> 2A'*A*h-2A'b = 0 <=> A'*A*h = A'*b <=> h = inv(A'*A)*A'*b.

Il suffit donc ici d'inverser une matrice 3x3 pour résoudre le problème et obtenir les coefficients de H. La matrices H seront calculées ainsi pour toutes les paires de points correspondants. Les points obtenus en multipliant un ensemble de points par H auront besoin d'être normalisés, c'est à dire qu'il faut que la troisième colonne des coordonnées homogènes de l'ensemble de points ne soit constituée que de 1. Pour cela on divise chaque colonne par la troisième colonne (il ne faut pas considérer les éventuels points dont la troisième coordonnée est nulle)

Déformation projective

Maintenant que les matrices H sont calculées, nous pouvons déformer une image pour la superposer à une autre image. Pour cela, nous avons créé une fonction warpImage qui, à partir d'une image et une homographie, renvoie l'image transformée par l'homographie. La première chose à connaître pour déformer une image est la taille de l'image déformée. Pour cela il suffit de calculer la projection des coins de l'image et de calculer la différence entre le minimum et le maximum des coordonnées. Nous avons choisi d'effectuer la transformation inverse, c'est à dire que pour chaque pixel de l'image de sortie, on va regarder la valeur correspondante dans l'image d'entrée. Pour cella on crée un vecteurs contenant toutes les coordonnées de l'image de sortie et on calcule sa transformation par H-1. On ne considèrera que les points dont la transformation est à l'intérieur de l'image de départ. Les coordonnées ne seront presque jamais des nombres entiers, il faudra alors calculer une interpolation linéaire sur l'image pour estimer les valeurs des canaux R, G et B pour des valeurs subpixelique. Ci-dessous nous présentons deux déformations de Lenna. La première est obtenue avec une homographie qui correspond à une rotation de 45° et la seconde est une homographie faisant intervenir plusieurs paramètres.


Les déformations d'images ont beaucoup d'applications différentes. Par exemple on pourra transformer les images pour améliorer la lisibilité de zones dont le plan n'est pas parallèle au plan de l'image. Ci-dessous nous présentons un exemple d'homographie pour améliorer la lisibilité d'une plaque d'immatriculation. Ce traitement peut être intéressant pour la lecture automatique des plaques d'immatriculations.



Mosaïque d'images

Maintenant que l'on connait les homographies pour passer d'une image à la suivante sur un ensemble d'image, on est capable de calculer les homographies pour passer d'une une image à n'importe qu'elle autre de l'ensemble d'image. En effet pour un ensemble de n images nous avons calculé les homographies Hi->i+1 ; donc l'homographie Hi->i+2 est donnée par Hi->i+2 = Hi->i+1*Hi+1->i+2 et l'homographie de l'image i+2 à l'image i est donné par l'inverse de cette matrice. Nous avons choisi de projeter toutes les images vers l'image centrale. Cependant pour projeter toutes les images sur un même plan il est nécessaire que le champ de vue du panorama soit inférieur à 180°. Les panorama du TP ont tous un champ de vue suppérieur à 180°. Nous les avons donc séparées en deux pour les deux premiers panoramas et nous avons pris 1/3 des images du panorama 3. De plus pour que le traitement soit plus rapide et pour éviter les complications dues à des mauvaises estimations d'homographies, nous projetons, l'image finale vers une image dont la dimension maximum est 2048 pixels. Nous n'avons effectué ce traitement que sur le panorama 2 et les panoramas des crédits supplémentaires. Voici les résultats obtenus pour les panoramas du TP.

Panorama 1



Panorama 2



Panorama 3

Ici nous avons coupé l'image pour une meilleure visibilité



Appariement de caractéristiques automatique

Pour éviter l'appariement manuel de points qui est une étape assez longue, nous proposons d'effectuer un appariement automatique des images. La première étape consiste à trouver des points d'intérêts dans les images. Ensuite il faut créer des descripteurs pour chaque point d'intérêt et enfin chercher les meilleurs appariements possibles. Nous décrivons ci-dessous chacune de ces étapes.

Détection de caractéristiques de coins dans une image


Pour sélectionner des points d'intérêts nous avons utilisé la fonction de Harris qui détecte les coins dans l'image. Plus exactement, la fonction de Harris calcule pour chaque pixel X de l'image: fHM(X) = det M(X) / tr M(X) avec M qui est le produit dyadique des gradients de X lissé par une Gaussienne. Enfin, les maximums locaux de fHM sont gardés en tant que points d'intérêts. Ce traitement permet d'obtenir beaucoup de points d'intérêts. Cependant il sera préférable de diminuer ce nombre de points d'intérêt pour réduire les temps de calculs d'appariement. La première idée est de garder les points d'intérêts ayant obtenus les meilleurs score pour fHM. L'inconvénient de cette technique est que l'on risque de garder des points très proches et donc difficilement différenciables lors de la recherche d'appariements. Nous avons donc programmé la répression maximale non adaptative (ANMS) proposé dans l'article "Multi-Image Matching using Multi-Scale Oriented Patches" de Brown et al. Cette méthode permet de garder des points d'intérêts bien espacés dans l'image. La première étape de de cette méthode est de garder le point maximisant le poids fHM, ensuite pour chaque points d'intérêts, on calcule r, la distance minimum aux autres points en ne considérant que ceux qui ont un poids 1.1 fois supérieur a celui du pixel considéré. On ne garde ensuite les pixels qui ont obtenus les plus grandes distances r.

500 points maximisant fHM 300 points par ANMS
r = 16.3

Extraction de descripteurs

Pour chaque point d'intérêt d'une image on va extraire un descripteur. Pour cela on va considérer une fenêtre 40x40 autour de chaque point caractéristique que l'on va échantillonner pour obtenir un groupement de pixel 8x8. On normalise ces 64 pixels en sous raillant par leur moyenne puis en divisant par l'écart type. Cette normalisation va permettre un appariement robuste aux changements d'intensité lumineuse. L'image ci-dessous montre un exemple de descripteur pour un point d'intérêt particulier.



Appariement de descripteurs de deux images

Lorsque les descripteurs de points caractéristiques de deux images ont été calculés, on recherche les correspondances en comparants les descripteurs de chaque image. Pour cela on calcule les erreurs pour tous les appariements possibles. Cette opération nous donne pour chaque point d'intérêt de l'image 1 les points les plus ressemblants dans l'image 2. On peut donc associer chaque point de l'image 1 à son meilleur correspondant. Cependant il est probable que certains points n'aient pas de correspondants, dans ce cas-là, l'erreur la plus faible pour ce point sera importante. On peut donc décider de ne garder que les appariements dont l'erreur est inférieure à un seuil. Il est cependant plus efficace de ne garder les appariements dont le ratio entre l'erreur du meilleur appariement et l'erreur du deuxième meilleur appariement est inférieur à un seuil. On choisira ce seuil à 0.1. Un exemple des appariements sélectionnés par cette méthode est montré ci-dessous:



Calcul robuste de l'homographie

Jusqu'à présent, nous avons calculé des points de correspondance entre les deux images, nous pouvons donc calculer l'homographie à partir de ces correspondances comme vu précédemment. Cependant, il se peut qu'il y ait encore de mauvais appariement entre deux images, ce qui donnera une mauvaise estimation de l'homographie. On va donc chercher l'homographie qui maximise le nombre d'appariements corrects. Pour cela on utilise une méthode appelée RanSaC dont la procédure est la suivante. On choisit en premier lieu un nombre d'itération. Plus ce nombre sera grand, plus l'algorithme sera robuste. A chaque itération on choisit au hasard 4 correspondances de points distinctes. On calcule l'homographie à partir de ces correspondances et calcule le nombre de points qui vérifient SSD(p',Hp) < eps avec SSD la somme des différences au carré et eps un seuil que nous avons choisi à 1 pour ne garder que les appariements qui sont justes à un pixel près. Lorsqu’on applique ce traitement aux points appareillés obtenus dans l'exemple précédent, on supprime 5 appariements (voir image ci-dessous).





Crédits supplémentaires

Photographies personnelles

Voici des exemples de panoramas à partir de rafales de 10 images aux chutes de Montmorency.

Mélange de saisons

Hiver


Crédit Andreas Wiesner



Été


Crédit Angel Marino



Panorama

Panorama avec deux fois la même personne

Descripteurs invariants par rotation

Nous avons programmé des descripteurs invariants par rotation. Pour cela, on regarde l'orientation du gradient au pixel considéré et on oriente la fenêtre du descripteur (40x40) selon cette orientation. Ce traitement donne moins d'appariements mais ceux-ci sont robustes aux changements d'orientations. La figure ci-dessous montre les appariements trouvés pour une image et sa rotation de 45° avant la recherche d'homographie par RanSaC. Lorsqu'on recherche les appariements sur les mêmes images sans l'amélioration des invariances par rotation on ne trouve aucun appariement.

Incrustation d'images et de vidéos

Nous avons créé une fonction permettant d'incruster une image dans une autre en sélectionnant où se trouverons les 4 coins de l'image à incruster dans l'autre image. Dans l'image ci-dessous nous avons incrusté un graffiti dans une des images des panoramas.



Nous avons aussi créé une fonction qui fait la même chose avec une vidéo. Dans l'image ci-dessous nous montrons un exemple d'incrustation de vidéo (premier plan) et d'incrustation d'image (second plan) sur des panneaux publicitaires. Pour que les bords de l'image soient adoucis, nous avons convolué le Masque correspondant à l'emplacement de l'image à incruster avec une gaussienne.