[SIFT] Développement de Taylor

**Noeim** · 13/02/2017, 19h48

Bonjour à tous !

Je cherche actuellement à programmer (en python) la méthode SIFT décrite ici :http://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf (page 10)

Mais l'utilisation du développement de Taylor n'est pas détaillé, j'ai donc cherché un peu partout des précisions mais il y a toujours des points obscures.

Si je comprend bien on considère que D, un DoG (différence de Gaussienne), dépend de x,y,sigma : D(x,y,sigma)
D'où mon premier problème : Puisque D est une différence de Gaussienne, donc de deux images de valeurs de sigma différentes, comment définir sigma dans D(x,y,sigma) ?

Mon second problème survient lorsque je cherche à faire la dérivée seconde de D par rapport à sigma : j'obtiens du D(x,y,sigma+2) et du D(x,y,sigma-2) sauf que ces valeurs de sigma ne sont pas définies(le problème est le même avec une dérivée simple). Peut-être que je devrais faire en sorte d'obtenir des valeurs connues à la place de sigma+2 et sigma-2 mais encore une fois je trouve que ça n'a pas de sens de définir un sigma sur un DoG... Ce qui revient, un peu, à ma première question.

Merci d'avance pour votre aide.

**pseudocode** · 13/02/2017, 20h47

Envoyé par Noeim

Mon second problème survient lorsque je cherche à faire la dérivée seconde de D par rapport à sigma : j'obtiens du D(x,y,sigma+2) et du D(x,y,sigma-2) sauf que ces valeurs de sigma ne sont pas définies(le problème est le même avec une dérivée simple).

Marrant. Il y a 5 ans, on m'avait posé la même question.

Et la réponse reste la même:

Il n'y a pas de 's+2' dans la formule des différences finies. c'est 's+1', 's,' et 's-1'.

D'ailleurs, toute la discussion de l'époque devrait répondre à tes questions.

Bon courage.

**Noeim** · 13/02/2017, 23h43

Merci d'avoir répondu si vite.

Effectivement ça m'a beaucoup aidé !
Et donc si j'ai bien compris, D ne dépend pas vraiment de sigma, la troisième variable c'est juste comme un indice "i" du DoG en question : D(x,y,i) avec i appartenant à {1,2,3,4} pour une octave (si s=2)
Pour la première octave ça donne :
D(x,y,1)= L(x,y,k*sigma)-L(x,y,sigma)
D(x,y,2)= L(x,y,k*k*sigma)-L(x,y,k*sigma)
...

Et comme les dérivées nous imposent d'avoir au moins un DoG d'indice inférieur et un autre d'indice supérieur on ne peut appliquer ce développement de Taylor qu'aux 2 DoG situé au milieu de notre octave (si s=2)

Oh et puisque j'y suis, j'ai encore 2 petites questions:

1) Lorsqu'on construit la pyramide, pour la deuxième octave par exemple, le floue gaussien s'effectue : à partir de la première image de l'octave ou, à partir de l'image initiale mais que pour un pixel sur deux ?
Je pencherais plutôt pour la deuxième possibilité, parce que sinon ça risque de fausser l'échelle... mais je préfère être sûr.

Selon les sources, les DoG sont construit en conservant les valeurs négatives de la différence d'image et d'autres non (les DoG si sombre hors des contours me font penser ça)
2) Faut-il conserver les valeurs négatives ou les mettre à 0 ou prendre les valeurs absolue ?

**pseudocode** · 14/02/2017, 01h25

Envoyé par Noeim

1) Lorsqu'on construit la pyramide, pour la deuxième octave par exemple, le floue gaussien s'effectue : à partir de la première image de l'octave ou, à partir de l'image initiale mais que pour un pixel sur deux ?
Je pencherais plutôt pour la deuxième possibilité, parce que sinon ça risque de fausser l'échelle... mais je préfère être sûr.

Le première image de l'octave c'est un downsample (un pixel sur deux) de la dernière image de l'octave précédent.

Faire un downsample ca crée une image 4 fois plus petite... c'est équivalent, en terme d'échelle, à l'image d'origine avec un sigma² multiplié par 4.
(Il suffit de remplacer x par x/2 et y par y/2 dans la formule de la gaussienne pour s'en rendre compte).

Et comme sur un octave on fait justement varier l'échelle de 1 à 4, alors le downsample de la dernière image (avec l'échelle = 4) devient la première image de l'octave suivant (avec une échelle ramenée à 1). On peut ainsi recommencer le processus sur le nouvel octave en faisant varier l'échelle de 1 à 4, etc.

J'avais posté un exemple de code dans la rubrique "contribuez" sur la construction des pyramides gaussiennes:
https://www.developpez.net/forums/d7...ne-scale-space

Note: cet exemple ne calcule pas les 2 images supplémentaires nécessaires pour calculer les DoG aux extrémités de l'octave.

Selon les sources, les DoG sont construit en conservant les valeurs négatives de la différence d'image et d'autres non (les DoG si sombre hors des contours me font penser ça)
2) Faut-il conserver les valeurs négatives ou les mettre à 0 ou prendre les valeurs absolue ?

D'une manière générale il faut conserver les valeurs négatives.

Mais je ne me souviens plus trop du détecteur de Lowe donc je ne sais pas si on peut "optimiser" et garder seulement les amplitudes.
Dans le doute, conserve les valeurs négatives.

**Noeim** · 18/02/2017, 15h41

Bonjour et encore merci pour votre aide.

Mais il se trouve que les échelles me posent toujours problème.
En effet une fois qu'on calcul le ^x (on note X) à partir du maximum (x,y,s) trouvé précédemment, on doit vérifier que ses coordonnées X=(^x,^y,^s) soit comprises entre -0.5 et 0.5 et si ce n'est pas le cas on modifie la valeur de x,y ou s de + ou - 1 et on recalcule X à partir de ces nouvelles valeurs.
Mais pour s qu'est-ce que c'est que s-1 (ou s+1) ? Comment refaire le calcul de X après une modification de s ?
J'ai l'impression de tournée autour du même problème à chaque fois mais même si je comprend qu'on cherche une position réelle, on n'a pas de quoi calculer les dérivées avec ces valeurs de s...

Et je pense que j'aurais le même problème quand il faudra, pour traiter l'invariance par rapport à l'illumination, calculer |D(X)| ? Par une approximation aux valeurs connues les plus proches ? Et, encore, comment arrondir s (on ne connait que très peu de valeur de D(s)) ?

En espérant être compréhensible.
Merci d'avance !

**pseudocode** · 19/02/2017, 02h22

Envoyé par Noeim

En effet une fois qu'on calcul le ^x (on note X) à partir du maximum (x,y,s) trouvé précédemment, on doit vérifier que ses coordonnées X=(^x,^y,^s) soit comprises entre -0.5 et 0.5 et si ce n'est pas le cas on modifie la valeur de x,y ou s de + ou - 1 et on recalcule X à partir de ces nouvelles valeurs.
Mais pour s qu'est-ce que c'est que s-1 (ou s+1) ? Comment refaire le calcul de X après une modification de s ?
J'ai l'impression de tournée autour du même problème à chaque fois mais même si je comprend qu'on cherche une position réelle, on n'a pas de quoi calculer les dérivées avec ces valeurs de s...

C'est peut être plus compréhensible si on change les notations:

1. soit une fonction tridimensionnelle D(x,y,s).
2. soit un point P(x,y,s) et un petit déplacement autour de ce point: dp(dx,dy,ds)

3. Le développement de taylor de D autour du point P nous donne: D(P+dp) = V + G.dp + 0.5*H.dp²

(avec V=valeur de D au point P, G=Gradient de D au point P et H=Hessian de D au point P)

4. Si "P+dp" est le point qui maximise D, alors la dérivée de D en ce point est nulle:
d( D(P+dp) )/dp = 0
d( D + G.dp + 0.5*H.dp² ) / dp = 0
0 + G + H.dp = 0
dp = - H^-1 . G

=> on estime G et H au point P par les différences finies (= les 27 valeurs de D(x-1|x|x+1 , y-1|y|y+1 , s-1|s|s+1) )
=> on calcule dp en resolvant le système "G + H.dp = 0" (système 3 équations à 3 inconnues: dx,dy,ds)
=> on estime le point qui maximise D: P+dp

**Noeim** · 19/02/2017, 13h36

D'accord jusque là je suis. Sauf peut-être pour G qui est une gaussienne de D ? Dans l'expression que j'avais c'était, à la place de G, le vecteur colonne :
| dD/dx |
| dD/dy |
| dD/ds |

Mais si une des valeurs de dp n'est pas comprise entre -0.5 et 0.5
je pensais qu'il fallait réitérer le calcul que tu as développer avec un nouveau point P.
Peut-être que j'ai mal compris ce passage du document de Lowe :

If the offset ˆx is larger than 0.5 in any dimension, then it means that the extremum lies closer to a different sample point. In this case,
the sample point is changed and the interpolation performed instead about that point.

Je donne un exemple du cas qui me pose problème :
On prend le point P(x,y,s) de D on cherche P+dp qui va maximiser D
On fait le calcul de dp= -H^(-1) . G
qui nous donne dx,dy,ds
et on suppose que seul |ds|>0.5 . On pose donc P'=(x,y,s+1)
On recalcule dp'.
Sauf que je ne sais pas comment calculer H et G pour ce nouveau point puisque qu'il a changé d'échelle et qu'on n'a pas calculer D pour cette nouvelle échelle.
Peut-être que c'est justement un lien entre G et D que je n'ai pas fais ... je sais pas. Mais je ne le vois vraiment pas dans le document de Lowe.

**pseudocode** · 19/02/2017, 18h59

Envoyé par Noeim

Sauf peut-être pour G qui est une gaussienne de D ? Dans l'expression que j'avais c'était, à la place de G, le vecteur colonne :

heu... oui. C'est moi qui me suis trompé dans mon message.

G c'est bien sur le Gradient de D au point P.

Mais si une des valeurs de dp n'est pas comprise entre -0.5 et 0.5
je pensais qu'il fallait réitérer le calcul que tu as développer avec un nouveau point P.

Si |dp| est plus grand que 0.5, ca veut dire que "P+dp" est plus proche d'un autre point Q de la pyramide.
On peut donc refaire les calculs en cherchant "Q+dq", avec de bonne chance que |dq| soit moins grand que 0.5

Note: P et Q sont des points de la pyramide => ils ont des coordonnées entières.

Sauf que je ne sais pas comment calculer H et G pour ce nouveau point puisque qu'il a changé d'échelle et qu'on n'a pas calculer D pour cette nouvelle échelle.

On a les valeurs de D pour toutes les coordonnées x,y,s entières. Ce sont toutes les valeurs de la pyramide.

Peut-être que c'est justement un lien entre G et D que je n'ai pas fais ... je sais pas.

G (Gradient) = Vecteur des dérivées partielles d'ordre 1.
H (Hessian) = matrice des dérivées partielles d'ordre 2.

En utilisant les différences finies:

Gradient:
D(x,y,s)/dx = ( D(x+1,y,s) - D(x-1,y,s) ) /2
D(x,y,s)/dy = ( D(x,y+1,s) - D(x,y-1,s) ) /2
...

Hessian:
D²(x,y,s)/dx² = ( D(x+1,y,s) - 2*D(x,y,s) + D(x-1,y,s) )
D²(x,y,s)/dxds = ( D(x+1,y,s+1) - D(x+1,y,s-1) - D(x-1,y,s+1) + D(x-1,y,s-1) ) / 4
..

**Noeim** · 19/02/2017, 23h27

Ah merci ! Je pense avoir compris la partie calcul

mais c'est dans son application que j'ai encore un soucis.

heu... oui. C'est moi qui me suis trompé dans mon message.

Ah ok c'est rassurant

Sauf que je ne sais pas comment calculer H et G pour ce nouveau point puisque qu'il a changé d'échelle et qu'on n'a pas calculer D pour cette nouvelle échelle.

On a les valeurs de D pour toutes les coordonnées x,y,s entières. Ce sont toutes les valeurs de la pyramide.

Je ne suis pas sûr de comprendre ce que c'est de passer de D(x,y,s) à D(x,y,s-1) ?
Je m'explique : avec s=2 cad k=2^(1/s)
On a pour chaque octave 3+s-1=4 différences de gaussienne.
Octave 1 : D(1,1) D(1,2) D(1,3) D(1,4)
Octave 2 : D(2,1) D(2,2) D(2,3) D(2,4)
Octave 3 : D(3,1) D(3,2) D(3,3) D(3,4)

Donc si je cherche P+dp dans D(1,2) et que j'obtiens ds<-0.5 je devrais alors, si j'ai bien compris, calculer Q+dq avec les mêmes x et y que P mais dans D(1,1) ?
Parce que là sa me pose un problème dans les calculs de G et H : il n'existe pas de D(1,0). Je devrais alors supprimé ce point ?

**pseudocode** · 20/02/2017, 00h41

Envoyé par Noeim

Donc si je cherche P+dp dans D(1,2) et que j'obtiens ds<-0.5 je devrais alors, si j'ai bien compris, calculer Q+dq avec les mêmes x et y que P mais dans D(1,1) ?
Parce que là sa me pose un problème dans les calculs de G et H : il n'existe pas de D(1,0). Je devrais alors supprimé ce point ?

Oui, si dans les calculs tu as besoin de la valeur de D en dehors de la pyramide, tu peux oublier le point.

Cela dit, dans ton exemple, je pense comprendre que D(1,0) c'est l'image originale (sans convolution par une gaussienne).

**Noeim** · 20/02/2017, 09h27

Ok ça marche !

Envoyé par pseudocode

Cela dit, dans ton exemple, je pense comprendre que D(1,0) c'est l'image originale (sans convolution par une gaussienne).

Euh... les D(1,1) D(1,2)... D(3,4) sont des différences de gaussiennes donc D(1,0) n’existe pas.
Dans mon exemple, si on pose G(1,1) et G(1,2) les deux premières images de la première octave alors D(1,1) = G(1,2) - G(1,1)

**pseudocode** · 20/02/2017, 10h32

Envoyé par Noeim

Ok ça marche !

Euh... les D(1,1) D(1,2)... D(3,4) sont des différences de gaussiennes donc D(1,0) n’existe pas.
Dans mon exemple, si on pose G(1,1) et G(1,2) les deux premières images de la première octave alors D(1,1) = G(1,2) - G(1,1)

Ah, ok pour ta notation. Du coup:

D(1,1) = G(1,2) - G(1,1)
D(1,0) = G(1,1) - G(1,0) = G(1,1) - Original

Envoyé par wikipedia

The scale parameter t=sigma² is the variance of the Gaussian filter and as a limit for t=0 the filter G becomes an impulse function such that L(x,y;0)=f(x,y), that is, the scale-space representation at scale level t=0 is the image f itself.

https://en.wikipedia.org/wiki/Scale_space

**Noeim** · 20/02/2017, 15h49

Envoyé par pseudocode

D(1,1) = G(1,2) - G(1,1)
D(1,0) = G(1,1) - G(1,0) = G(1,1) - Original

Ah d'accord ! Donc je me suis trompé dans la construction de la pyramide, décidément ...
Donc je reprend tout depuis le début avec, k= 2^1/2 et sigma0=1.6=s0
On construit, en notant L le résultat de la convolution par la gaussienne, D les DoG et en dessous des L les valeurs de "flou" (utilisé pour la gaussienne) appliqué à la première image de chaque octave :

Octave 1 : L(1,0)=(image originale) L(1,1) L(1,2) L(1,3) L(1,4)
flou---------1----------------------s0.k---s0.k²--s0.k³--s0.k⁴
D(1,0) = L(1,1) - L(1,0)
D(1,1) = L(1,2) - L(1,1)
D(1,2) = L(1,3) - L(1,2)
D(1,3) = L(1,4) - L(1,3)

Octave 2 : L(2,0)=(L(1,2) en prenant 1 pixel sur 2) L(2,1) L(2,2) L(2,3) L(2,4)
flou---------1------------------------------------s0.k¹--s0.k²--s0.k³--s0.k⁴
D(2,0) = L(2,1) - L(2,0)
D(2,1) = L(2,2) - L(2,1)
D(2,2) = L(2,3) - L(2,2)
D(2,3) = L(2,4) - L(2,3)

Octave 3 : L(3,0)=(L(2,2) en prenant 1 pixel sur 2) L(3,1) L(3,2) L(3,3) L(3,4)
flou---------1------------------------------------s0.k¹--s0.k²--s0.k³--s0.k⁴
D(3,0) = L(3,1) - L(3,0)
D(3,1) = L(3,2) - L(3,1)
D(3,2) = L(3,3) - L(3,4)
D(3,3) = L(3,4) - L(3,3)

(je sais que tu as donner un exemple de construction de pyramide mais comme je ne connais pas java...)

On recherche les maximas entre les 26 voisins pour D(1,1) D(1,2) et D(2,1) D(2,2) et D(3,1) D(3,2).
Et on fait les calculs de recherche de maxima avec les positions réelles sur les mêmes D.
Si le résultat nous amène à calculer dans un autre D que ceux juste au dessus on abandonne ce point.

J’espère que cette fois c'est bon

Et encore merci pour ton aide !