Invariance à la Rotation avec SIFT (Matlab)

Version imprimable

Bonjour,

J'ai créé un programme composé en trois étapes :
-Détecter les points d’intérêts de deux images avec Harris
-Trouver les descripteurs de Harris
- Matcher les points identiques.

Les points de harris sont invariants à la rotation mais pas ceux du descripteur SIFT, et c'est là ou réside mon problème. J'ai créée une fonction matlab pour générer mon descripteur et mon invariance à la rotation par recalage d'histogramme mais cela ne fonctionne pas... Pouvez vous m'aider ?
Ci - après mon code source :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
function [H] =calculSIFTBloc(ptsHarris,IgPad,IorPad,k,Mg,kIndice)
%Créer l'histogramme d'orientation pour un point et un bloc donné
% @param ptsHarris : la liste des points
% @param IgPad : le module du gradiant paddé
% @param IorPad : l'image discrétisé en 8 Bins paddé
% @param k : le numéro de bloc
% @param Mg : the Gaussian SIFT mask
% @param kIndice : matrice d'indexation de bloc
 
% Coordonées du points de Harris
x=ptsHarris(1,1);
y=ptsHarris(1,2);
 
% Size of padding
t=8;
 
%% Trouve le voisinage de 16x16 au point indiqué
 
for i=8:(size(IorPad,1)-8)
    for j=8:(size(IorPad,2)-8)
 
        if (i-8==x && j-8==y)
 
            blocIorPad(:,:)=IorPad((i-t):(i+t -1),(j-t):(j+t -1));
            blocIgPad(:,:) =IgPad ((i-t):(i+t -1),(j-t):(j+t -1));
 
            %%% Mode DEBUG %%%
            save('blocIorPad.txt','blocIorPad', '-ascii');
            save('blocIgPad.txt','blocIgPad', '-ascii');
 
        end
    end
end
 
%%% Mode DEBUG %%%
load('blocIorPad.txt');
load('blocIgPad.txt');
 
%% Invariance a la rotation
 
% Calcul l'histogramme discrétisé en 8 bin 
HRot=zeros(1,8);
 
 
 MatH=blocIgPad.*Mg;
 
for i=1:16
    for j=1:16
 
        valH1=blocIorPad(i,j);
 
 
        if (valH1 ~= 0)
            HRot(valH1)=HRot(valH1)+MatH(i,j);
        end
 
    end
end
 
% Trouve le maximum du descripteur de Ig
maxHRot=max(max(HRot));
 
 
for i=1:8
    if (HRot(i)==maxHRot)
 
        thetaiIor=i;
 
    end
end
%[val,pos] = max(C(:));
 
%% Calcul le descripteur SIFT en fonction du sous bloc défini
 
% Trouve les indices du sous-bloc
valInd=kIndice(:,:,k);
valx=valInd(1,1);
valy=valInd(1,2);
 
% Trouve les sous-blocs
 ssBlocIorPad(:,:)=blocIorPad(valx:valx+3,valy:valy+3);
 ssBlocMatCste(:,:)=MatH(valx:valx+3,valy:valy+3);
 
 
 
% Calcul le descripteur SIFT NON Invariant à la rotation
HNonrect=zeros(1,8);
 
for i=1:4
    for j=1:4
        val1=ssBlocIorPad(i,j);
 
          if (val1 > 0)
             HNonrect(val1)=HNonrect(val1)+ssBlocMatCste(i,j);
          end
    end
end
 
%Recalage Géométrique
Hrect=[HNonrect HNonrect];
 
H=Hrect(thetaiIor:thetaiIor+7);
end

09/05/2011, 22h56
TheDwarf01

Ne me dite pas qu'aucun modérateur ou rédacteur du forum ne sait résoudre mon problème .... :cry:

Ps: Si vous pouvez éviter la référence à l'article de Lowe :ccool:
10/05/2011, 11h20
pseudocode

Citation:

Envoyé par TheDwarf01

Ps: Si vous pouvez éviter la référence à l'article de Lowe :ccool:

C'est dommage, car si on pouvait faire référence à l'article de Lowe (chap 5, page 13), on pourrait expliquer que le descripteur SIFT est invariant aux rotations. En effet, on fait toujours une rotation du voisinage pour qu'il coïncide avec l'orientation du gradient.
10/05/2011, 13h57
TheDwarf01

L'invariance à la rotation dans l'article de Lowe prend en considération le multi-échelle et un histogramme discrétisé en 36bins.

Dans mon cas mon orientation est discrétisé en 8bins ce qui n'influence pas en grande partie les résultats et fait juste appel à un détecteur de Harris Mono Echelle.

Ce que je voudrais savoir c'est si le recalage géométrique se fait bien de cette façon :

-Prendre une fenêtre de voisinage 16x16 dans ma matrice discrétisé en 8 bins au point de harris détecter => blocIorPad
-Prendre une fenêtre de voisinage 16x16 dans ma matrice de module du gradient au point de harris détecté => blocIgPad
-Créer un histogramme pour Ig pour un voisnage 16*16 => trouver le thetai position de la valeur maximal de mon histo 1<=thetai<=8

-Diviser mon bloc "blocIorPad" et "blocIgPad" en 16 sous blocs de 4x4 et calculer un histogramme pour un numéro de bloc donnée donc un 1<=kbloc<=16.
-Recaler l'histogramme de mon kbloc en fonction de mon thetai de cette maniere :

Citation:

%Recalage Géométrique
Hrect=[Histokbloc Histokbloc ];

H=Hrect(kbloc:kbloc+7);

A NOTER : Mon descripteur de Harris fonctionne bien avec deux images translaté (recalage géométrique inclus dans le code de compilation), mais me donne un matching à 100% faux dans le cas d'une image ayant une rotation avec même un angle de 5°... Or je devrais au moins trouver un matching correct à 80% min.

Donc auriez vous une idée pour un recalage géométrique correct ?
10/05/2011, 14h22
pseudocode

Citation:

Envoyé par TheDwarf01

L'invariance à la rotation dans l'article de Lowe prend en considération le multi-échelle et un histogramme discrétisé en 36bins.

Le multi-échelle n'intervient que dans l'étape de localisation des keypoints, pas dans la phase d'orientation.

L'histogramme est effectivement discrétisé en 36 bins de 10°, mais la valeur de l'orientation dominante est trouvée par un fit sur une parabole : ca donne donc une valeur précise de l'orientation et pas l'index dans l'histogramme.

Citation:

Dans mon cas mon orientation est discrétisé en 8bins ce qui n'influence pas en grande partie les résultats et fait juste appel à un détecteur de Harris Mono Echelle.

Ce que je voudrais savoir c'est si le recalage géométrique se fait bien de cette façon :

-Prendre une fenêtre de voisinage 16x16 dans ma matrice discrétisé en 8 bins au point de harris détecter => blocIorPad
-Prendre une fenêtre de voisinage 16x16 dans ma matrice de module du gradient au point de harris détecté => blocIgPad
-Créer un histogramme pour Ig pour un voisnage 16*16 => trouver le thetai position de la valeur maximal de mon histo 1<=thetai<=8

-Diviser mon bloc "blocIorPad" et "blocIgPad" en 16 sous blocs de 4x4 et calculer un histogramme pour un numéro de bloc donnée donc un 1<=kbloc<=16.
-Recaler l'histogramme de mon kbloc en fonction de mon thetai de cette maniere :

Il faut faire une rotation du voisinage avant de diviser le bloc et sous-blocs.

http://www.cg.tu-berlin.de/fileadmin...rientation.png

Pour le reste, le principe m'a l'air bon. Pour le code Matlab, je ne sais pas : je ne lis pas ce langage. :aie:

Citation:

A NOTER : Mon descripteur de Harris fonctionne bien avec deux images translaté (recalage géométrique inclus dans le code de compilation), mais me donne un matching à 100% faux dans le cas d'une image ayant une rotation avec même un angle de 5°... Or je devrais au moins trouver un matching correct à 80% min.

Donc auriez vous une idée pour un recalage géométrique correct ?

Le détecteur de Harris "de base" n'est pas très robuste. Comme le coté multi-echelle de SIFT ne semble pas t'intéresser, jette un coup d'oeil sur le détecteur FAST : il est très rapide, très simple à programmer et beaucoup plus robuste.
16/05/2011, 16h56
TheDwarf01

Citation:

Il faut faire une rotation du voisinage avant de diviser le bloc et sous-blocs.

Je tourne le voisinage autour de mon point d'intérêt de l'angle d'orientation correspondant à l'indice de mon theta i dans mon histogramme d'orientation sur le voisinage de 16x16 dans mon cas c'est bien ça ?

Mais dois je aussi le tournée pour le module du gradient de l'image ?

D'après wikipédia la rotation doit être de sens inverse pourquoi cela ?
16/05/2011, 17h18
pseudocode

Citation:

Envoyé par TheDwarf01

Je tourne le voisinage autour de mon point d'intérêt de l'angle d'orientation correspondant à l'indice de mon theta i dans mon histogramme d'orientation sur le voisinage de 16x16 dans mon cas c'est bien ça ?

Ta phrase est un poil compliquée, mais je pense que c'est cela.

Citation:

Mais dois je aussi le tournée pour le module du gradient de l'image ?

On ne tourne pas le "module" (= la grandeur) mais l'angle. Ca revient à faire une soustraction de l'angle du gradient avec l'angle de la région. Si ton gradient à un angle de 25°, et que la région est également orientée à 25°, alors ton gradient forme un angle de 0° par rapport à la région.

Citation:

D'après wikipédia la rotation doit être de sens inverse pourquoi cela ?

Généralement, on "tourne" la région pour ramener que son orientation à 0°. Ce qui revient à faire tourner toute l'image en sens inverse. Ca permet de raisonner dans la région en terme de coordonnées verticales/horizontales habituelles (=> le point haut/gauche est en (-8,+8)).
21/05/2011, 19h32
TheDwarf01

J'ai un petite impression que mon erreur viendrai de là :

Si je discrétise l'orientation en 8 bins de cette façon :
orientation = [
1,0 ;
cos(pi/8), -sin(pi/8) ;
cos(pi/4) , -sin(pi/4) ;
cos(3*pi/8) , -sin(3*pi/8) ;
0, -1 ;
cos(5*pi/8) , -sin(5*pi/8) ;
cos(3*pi/4) , -sin(3*pi/4) ;
cos(7*pi/8) , -sin(7*pi/8) ;
-1,0 ];

Que je calcul mon histogramme sur mon patch de 16x16 sur l'image discrétisé pour un point d'intérêt et que par exemple mon indice imax = 4. Je dois donc tourner ce patch de :
Histo Ior =[ 1 2 20 44 2 7 8 9 ]
Angle degrès = [0 22.5 45 67.5 90 112.5 135 157.5] donc de -67.5°

J'ai aussi lu que lorsque je calcule mes 16 histogrammes, je devais retirer pour chaque pixel la distance euclidienne (du point d'intérêt par rapport au pixel en question). Est ce exact ?
22/05/2011, 11h22
pseudocode

Citation:

Envoyé par TheDwarf01

sur mon patch de 16x16 sur l'image discrétisé pour un point d'intérêt et que par exemple mon indice imax = 4. Je dois donc tourner ce patch de :
Histo Ior =[ 1 2 20 44 2 7 8 9 ]
Angle degrès = [0 22.5 45 67.5 90 112.5 135 157.5] donc de -67.5°

Lowe conseille de trouver l'index (imax) de l'orientation par fitting d'une parabole --> Dans ton exemple, il faut trouver la parabole qui passe par les points (3,2) (4,44) (5,2), et de là chercher son maximum. Ca nous donne un "imax" (et donc un angle) beaucoup plus précis.

Citation:

J'ai aussi lu que lorsque je calcule mes 16 histogrammes, je devais retirer pour chaque pixel la distance euclidienne (du point d'intérêt par rapport au pixel en question). Est ce exact ?

Oui, c'est ce qui est expliqué au chapitre 6.1 du papier de Lowe. Toutes les valeurs de magnitude dans le patch 16x16 sont pondérées par une gaussienne centrée sur le keypoint
23/05/2011, 00h11
TheDwarf01

Citation:

Envoyé par pseudocode

Dans ton exemple, il faut trouver la parabole qui passe par les points (3,2) (4,44) (5,2), et de là chercher son maximum. Ca nous donne un "imax" (et donc un angle) beaucoup plus précis.

Les valeurs de Histo Ior =[ 1 2 20 44 2 7 8 9 ] sont les nombres de point d'orientation de 0° à 157.5°.

Si j'ai bien compris je dois un faire un fitting des paraboles passant par les orientations +-1 par rapport a mon orientation maximal et trouver la valeur la plus grande entre ces 3 ajustements (fitting)... :? Mais pour avoir un angle plus précis il aurait fallu que je discrétise mon images Ior en 36 bins et non 8.

Je cherche juste à implémenter un SIFT très simplifié m'amenant à un matching correct dans un peu plus de 80% des points matchés. C'est pour celà que j'ai discrétisé seulement en 8bins Ior.
23/05/2011, 13h24
pseudocode

Citation:

Envoyé par TheDwarf01

Les valeurs de Histo Ior =[ 1 2 20 44 2 7 8 9 ] sont les nombres de point d'orientation de 0° à 157.5°.

Si j'ai bien compris je dois un faire un fitting des paraboles passant par les orientations +-1 par rapport a mon orientation maximal et trouver la valeur la plus grande entre ces 3 ajustements (fitting)... :? Mais pour avoir un angle plus précis il aurait fallu que je discrétise mon images Ior en 36 bins et non 8.

Oui, exact. Normalement dans SIFT c'est un histo 36 bins pour trouver l'orientation du patch.

Citation:

Je cherche juste à implémenter un SIFT très simplifié m'amenant à un matching correct dans un peu plus de 80% des points matchés. C'est pour celà que j'ai discrétisé seulement en 8bins Ior.

Il ne vaut mieux pas sacrifier l'estimation de l'orientation, car le descripteur en lui même n'est pas invariant à la rotation (les 16 sous-régions 4x4 ont une position précise). C'est pour cela qu'il vaut mieux orienter le patch avec précision, surtout que 8 ou 36 bins, ca ne fait pas une si grande différence.
26/05/2011, 21h55
TheDwarf01

Pour faire une rotation de mon patch manuellement sans passer par la fonction de matlab :

- Comment a partir d'un pixel je peu trouver le pixel à X° en passant par une matrice de rotation sachant que ça peut ne pas tomber sur un pixel parfaitement.

- Et comment retracer mon patch sachant que je traverse des pixels en diagonal ?
26/05/2011, 22h22
pseudocode

Citation:

Envoyé par TheDwarf01

- Comment a partir d'un pixel je peu trouver le pixel à X° en passant par une matrice de rotation sachant que ça peut ne pas tomber sur un pixel parfaitement.

On ne trouve pas "le" pixel correspondant car, comme tu le dis, le résultat du calcul nous donne des coordonnées non-entières. On tombe donc entre 4 pixels -> il faut interpoler pour trouver la valeur.

Citation:

- Et comment retracer mon patch sachant que je traverse des pixels en diagonal ?

:koi: ?
26/05/2011, 23h00
TheDwarf01

1 pièce(s) jointe(s)

-D'accord mais sachant qu'une matrice de rotation est défini comme suit :

M=(cos(A) -sin(A) )
sin(A) cos(A)

Sachant que moi j'ai un point (x,y) je l'adapte comment pour trouver les nouveaux coordonnées et en utilisant quel matrice la première en M où une matrice de rotation sous y , y ou z :

cf ici

-Sinon quand je souhaite reformer mon patch de 16x16 je passe de l'image de (cf ma piece jointe ) l'image de gauche a celle de droite et donc je doit bien trouver les pixels en 'diagonale'
26/05/2011, 23h19
pseudocode

Citation:

Envoyé par TheDwarf01

-D'accord mais sachant qu'une matrice de rotation est défini comme suit :

M=(cos(A) -sin(A) )
sin(A) cos(A)

Sachant que moi j'ai un point (x,y) je l'adapte comment pour trouver les nouveaux coordonnées et en utilisant quel matrice la première en M où une matrice de rotation sous y , y ou z :

C'est une multiplication vecteur 2d * matrice de rotation 2d

x' = x.cos(a) - y.sin(A)
y' = x.sin(a) + y.cos(A)

Citation:

-Sinon quand je souhaite reformer mon patch de 16x16 je passe de l'image de (cf ma piece jointe ) l'image de gauche a celle de droite et donc je doit bien trouver les pixels en 'diagonale'

:D Non, les pixels ne sont pas en diagonale. Le dessin est juste là pour expliquer le principe, mais un pixel reste un "point" et il n'a donc pas d'orientation.
26/05/2011, 23h42
TheDwarf01

Citation:

Envoyé par pseudocode

C'est une multiplication
:D Non, les pixels ne sont pas en diagonale. Le dessin est juste là pour expliquer le principe, mais un pixel reste un "point" et il n'a donc pas d'orientation.

Merci pour le lien wiki.

Ho lalala c'est bon j'ai compris désolé grosse fatigue en ce moment...