[image] géométrie projective (homography)

[pseudocode]

Le code suivant est une implémentation Java des formules du document

Projective Mappings for Image Warping (1989) by Paul Heckbert

Cette implémentation simple permet de transformer une zone quadrilatère quelconque en une zone rectangulaire.



Le programme principal:

Code java :

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
 
public static void main(String[] args) {
	Image source = new Image("yourimage.jpg");
 
	// 4 points qui définissent un quadrilatère convexe
	Point P1 = new Point(10,71);
	Point P2 = new Point(134,71);
	Point P3 = new Point(71,24);
	Point P4 = new Point(62,24);
 
	// résolution du système linéaire
	double[] system = getSystem(P1,P2,P3,P4);
	System.out.println(Arrays.toString(system));
 
	// création de l'image corrigée
	int W=640,H=480;
	Image target = new Image(W,H);
 
	// pour chaque pixel (x,y) de l'image corrigée
	for(int y=0;y<H;y++) {
		for(int x=0;x<W;x++) {
 
			// conversion dans le repère orthonormé (u,v) [0,1]x[0,1]
			double u = (double)x/W;
			double v = (double)y/H;
 
			// passage dans le repère perspective
			double[] P = invert(u, v, system);
 
			// copie du pixel (px,py) correspondant de l'image source 
			// TODO: faire une interpolation pour avoir de meilleurs résultats
			int px=(int)Math.round(P[0]);
			int py=(int)Math.round(P[1]);
			target.setrgb(x, y, source.getrgb(px, py) );
		}
	}
}

Et les algos de calcul issus du document dont j'ai donné le lien ci-avant

Code java :

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
 
public static double[] getSystem(Point... P) {
	double[] system = new double[8];
	double sx = (P[0].x-P[1].x)+(P[2].x-P[3].x);
	double sy = (P[0].y-P[1].y)+(P[2].y-P[3].y);
	double dx1 = P[1].x-P[2].x;
	double dx2 = P[3].x-P[2].x;
	double dy1 = P[1].y-P[2].y;
	double dy2 = P[3].y-P[2].y;
 
	double z = (dx1*dy2)-(dy1*dx2);
	double g = ((sx*dy2)-(sy*dx2))/z;
	double h = ((sy*dx1)-(sx*dy1))/z;
 
	system[0]=P[1].x-P[0].x+g*P[1].x;
	system[1]=P[3].x-P[0].x+h*P[3].x;
	system[2]=P[0].x;
	system[3]=P[1].y-P[0].y+g*P[1].y;
	system[4]=P[3].y-P[0].y+h*P[3].y;
	system[5]=P[0].y;
	system[6]=g;
	system[7]=h;
 
	return system;
}
 
public static double[] invert(double u, double v, double[] system) {
	double x = (system[0]*u+system[1]*v+system[2])/(system[6]*u+system[7]*v+1); 
	double y = (system[3]*u+system[4]*v+system[5])/(system[6]*u+system[7]*v+1);
	return new double[]{x,y};
}

[pseudocode]

Une implémentation plus complète pour OpenCV fournie par luxigo sur le forum.

Code C :

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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
 
/*
 L'algorithme vient de:
 
 Projective Mappings for Image Warping (1989)
 by Paul Heckbert
 in Fundamentals of Texture Mapping and Image Warping (Paul Heckbert,
 Master’s Thesis), U.C.Berkeley
 
 http://www.cs.utah.edu/classes/cs6964-whitaker/heckbert_projective.pdf 
 
 Le cas des parallèlogrames n'est pas pris en compte dans cette implémentation.
 
 Merci à pseudocode pour le débroussaillage :)
 
*/
 
double *getPerspectiveTransform(CvPoint2D32f *P) {
 
        double *H=calloc(18,sizeof(double));
        double *adj=H+9;
 
        double sx = (P[0].x-P[1].x)+(P[2].x-P[3].x);
        double sy = (P[0].y-P[1].y)+(P[2].y-P[3].y);
        double dx1 = P[1].x-P[2].x;
        double dx2 = P[3].x-P[2].x;
        double dy1 = P[1].y-P[2].y;
        double dy2 = P[3].y-P[2].y;
 
        double z = (dx1*dy2)-(dy1*dx2);
        double g = ((sx*dy2)-(sy*dx2))/z;
        double h = ((sy*dx1)-(sx*dy1))/z;
 
        // matrice de transformation
        double a=H[0]=P[1].x-P[0].x+g*P[1].x;
        double b=H[1]=P[3].x-P[0].x+h*P[3].x;
        double c=H[2]=P[0].x;
        double d=H[3]=P[1].y-P[0].y+g*P[1].y;
        double e=H[4]=P[3].y-P[0].y+h*P[3].y;
        double f=H[5]=P[0].y;
        H[6]=g;
        H[7]=h;
        H[8]=1;
 
        // matrice de transformation inverse (matrice adjointe)
        adj[0]=e-f*h;
        adj[1]=c*h-b;
        adj[2]=b*f-c*e;
        adj[3]=f*g-d;
        adj[4]=a-c*g;
        adj[5]=c*d-a*f;
        adj[6]=d*h-e*g;
        adj[7]=b*g-a*h;
        adj[8]=a*e-b*d;
 
        return H;
}
 
void projective_mapping(double *u, double *v, double *H) {
        double x = (H[0]*u[0]+H[1]*v[0]+H[2])/(H[6]*u[0]+H[7]*v[0]+1);
        double y = (H[3]*u[0]+H[4]*v[0]+H[5])/(H[6]*u[0]+H[7]*v[0]+1);
        u[0]=x;
        v[0]=y;
}
 
void inverse_projective_mapping(double *u, double *v, double *H) {
        double x = (H[0]*u[0]+H[1]*v[0]+H[2])/(H[6]*u[0]+H[7]*v[0]+1);
        double y = (H[3]*u[0]+H[4]*v[0]+H[5])/(H[6]*u[0]+H[7]*v[0]+1);
        double z = (H[6]*u[0]+H[7]*v[0]+H[8])/(H[6]*u[0]+H[7]*v[0]+1);
        u[0]=x/z;
        v[0]=y/z;
}
 
void perspectiveTransform(IplImage *img,CvPoint *I, double *H) {
 
        // conversion dans le repère orthonormé (u,v) [0,1]x[0,1]
        double u = (double)I->x/(double)img->width;
        double v = (double)I->y/(double)img->height;
 
        // passage dans le repère perspective
        projective_mapping(&u, &v, H);
 
        // pixel correspondant dans le quadrilatère
        I->x=(int)round(u);
        I->y=(int)round(v);
 
}
 
void inversePerspectiveTransform(IplImage *img,CvPoint *I, double *H) {
 
        double u=I->x;
        double v=I->y;
 
        // passage dans le repère orthonormé
        inverse_projective_mapping(&u, &v, H+9);
 
        // pixel correspondant dans le carré
        I->x=(int)round(u*(double)img->width);
        I->y=(int)round(v*(double)img->height);
 
}
 
/*
 
 L'image source est transformée pour tenir dans dst.
 
 Quand map n'est pas nul c'est un pointeur sur un tableau de
 dimension sizeof(CvPoint)*dst->width*dst->height
 qui contiendra les coordonées d'origine des pixels.
 
*/
 
void perspectiveTransformImage(IplImage *src,CvPoint2D32f *P,IplImage *dst,CvPoint *map,double *H) {
 
        double x;
        double y;
        double width=dst->width;
        double height=dst->height;
 
        // pour chaque pixel de l'image cible   
        for(y=0;y<height;y++) {
                for(x=0;x<width;x++) {
 
                        // conversion dans le repère orthonormé (u,v) [0,1]x[0,1]
                        double u = x/width;
                        double v = y/height;
 
                       // passage dans le repère perspective
                        projective_mapping(&u, &v, H);
 
                        // pixel correspondant de l'image source 
                        int sx=(int)round(u);
                        int sy=(int)round(v);
 
                        // pixel source 
                        CvScalar rgb=cvGet2D(src,sy,sx);
 
                        // pixel destination        
                        int dx=x;
                        int dy=y;
                        cvSet2D(dst,dy,dx,rgb);
 
                        // conserver position pixel source
                        if (map) {
                                size_t index=dy*dst->width+dx;
                                map[index].x=sx;
                                map[index].y=sy;
                        }
                }
        }
}

[mina777]

salut,je veux savoir quels sont les coefficients de la matrice homgraphique qui permet le passage d'une image non rectifiée à une image rectifiée(càd qcq ils présentent).
merci

[pseudocode]

Envoyé par mina777

salut,je veux savoir quels sont les coefficients de la matrice homgraphique qui permet le passage d'une image non rectifiée à une image rectifiée(càd qcq ils présentent).
merci

Je n'ai pas compris ta question. Est-ce que tu as lu le document PDF mentionné dans le premier post ?

[Dut]

Envoyé par pseudocode

Est-ce que tu as lu le document PDF mentionné dans le premier post ?

=> http://www.eng.utah.edu/~cs6964-whit...projective.pdf

Objet non trouvé!

[pseudocode]

Envoyé par Dut

=> http://www.eng.utah.edu/~cs6964-whit...projective.pdf

Ah! Heureusement il reste une copie du document sur citeseerx.

J'ai corrigé le lien.

[sniperMdAb]

Joli!!

Est-ce que ça pourrai être fait sous flash8 si on a les coordonnées des 4 points du premier quadrilatère?
En fait ce qu'il me faut c'est de déterminer les coordonnées d'un point du rectangle a partir de ses coordonnées dans un quadrilatère quelconque.

Ce truc me rend fou depuis quelque temps, je ne dort plus la nuit

merci d'avance pour votre aide!

[pseudocode]

Envoyé par sniperMdAb

Joli!!

Est-ce que ça pourrai être fait sous flash8 si on a les coordonnées des 4 points du premier quadrilatère?
En fait ce qu'il me faut c'est de déterminer les coordonnées d'un point du rectangle a partir de ses coordonnées dans un quadrilatère quelconque.

Ce truc me rend fou depuis quelque temps, je ne dort plus la nuit

merci d'avance pour votre aide!

nicoptere à déja fait le boulot pou toi sur son blog : http://fr.nicoptere.net/?p=351

[sniperMdAb]

WOW! Merci beaucoup l'ami c'est E X A C T E M E N T ce que je voulais!

[ToTo13]

Bonjour,

histoire de compliquer un peu les choses... est ce qu'il y a une version pour quatre points quelconques (donc pas les coins d'un quadrilatère) ?

[pseudocode]

Envoyé par ToTo13

Bonjour,

histoire de compliquer un peu les choses... est ce qu'il y a une version pour quatre points quelconques (donc pas les coins d'un quadrilatère) ?

Du genre "les 4 milieux des cotés" ? Je ne pense pas. Dans ce cas particulier on aurait une infinité de parallélogrammes possibles comme quadrilatère de départ. Et donc une infinité de déformation "skew" possible.

L'idée de base c'est quand meme que 3 des 4 points forment une base otho, et que le 4eme point serve a mesurer la déformation.

[ToTo13]

Envoyé par pseudocode

Du genre "les 4 milieux des cotés" ? Je ne pense pas. Dans ce cas particulier on aurait une infinité de parallélogrammes possibles comme quadrilatère de départ. Et donc une infinité de déformation "skew" possible.
L'idée de base c'est quand meme que 3 des 4 points forment une base otho, et que le 4eme point serve a mesurer la déformation.

mouai... ça aurait été trop beau...

En fait j'ai réussi à faire un tracking robuste d'un objet référence dans un flux, mais les points caractéristiques que je trouve se situent un peu partout dans l'objet et non pas seulement sur les bords. Bien entendu, j'ai l'appariement entre les points sur l'objet traqué et les points sur l'objet référence.
Donc la disposition des points que je trouve ne forme pas un quadrilatère.

[pseudocode]

Envoyé par ToTo13

mouai... ça aurait été trop beau...

En fait j'ai réussi à faire un tracking robuste d'un objet référence dans un flux, mais les points caractéristiques que je trouve se situent un peu partout dans l'objet et non pas seulement sur les bords. Bien entendu, j'ai l'appariement entre les points sur l'objet traqué et les points sur l'objet référence.
Donc la disposition des points que je trouve ne forme pas un quadrilatère.

C'est clair que de la géométrie projective simple (comme ici) ne permet pas de faire de l'estimation de pose 3D. Ici c'est simplement des déformations 2D.

[ToTo13]

Envoyé par pseudocode

C'est clair que de la géométrie projective simple (comme ici) ne permet pas de faire de l'estimation de pose 3D. Ici c'est simplement des déformations 2D.

Oui, c'est vrai, mais comme l'affichage est 2D, je dois peut être pouvoir truander le système (ok, j'ai l'habitude).
Dis moi ce que tu en penses :
- j'ai mes points traqués, répartis un peu partout au milieu de l'objet.
- mais j'ai aussi les références/appariement avec l'objet d'origine (donc non déformé par les différentes orientations 3D).
Donc avec cet ensemble (formé d'au moins quatre points), je peux sans doute extrapoler les coordonnées des quatre coins de mon objet, avec une bonne approximation.
Faisable ou je suis complètement utopiste ?

[pseudocode]

Envoyé par ToTo13

Oui, c'est vrai, mais comme l'affichage est 2D, je dois peut être pouvoir truander le système (ok, j'ai l'habitude).
Dis moi ce que tu en penses :
- j'ai mes points traqués, répartis un peu partout au milieu de l'objet.
- mais j'ai aussi les références/appariement avec l'objet d'origine (donc non déformé par les différentes orientations 3D).
Donc avec cet ensemble (formé d'au moins quatre points), je peux sans doute extrapoler les coordonnées des quatre coins de mon objet, avec une bonne approximation.
Faisable ou je suis complètement utopiste ?

Bah, heu... oui c'est faisable. C'est le principe de l'

estimation de pose

estimation de pose

.

Le truc c'est que je ne vois pas le rapport avec cette contribution.
Il y a des méthodes beaucoup plus adaptée a ton problème (cf. Long Quan and Zhongdan Lan)

[ToTo13]

Envoyé par pseudocode

Le truc c'est que je ne vois pas le rapport avec cette contribution.
Il y a des méthodes beaucoup plus adaptée a ton problème (cf. Long Quan and Zhongdan Lan)

Merci infiniment.
Je pensais (à tord) que la matrice d'homographie pourrait m'aider dans mon problème. Il est vrai que je ne connais rien au domaine :-(

Dans mon cas précis, c'est exactement cet article qu'il me faut.
Reste plus qu'à ce que j'arrive à le développer ou à trouver une bonne implémentation.

En tout cas, merci

[tommy94]

Bonjour et merci de m'aider sur ce petit projet
désolé de réouvrir ce topic 2 ans après mais j'ai un gros problème :/
voilà alors j'ai tapé le code que tu as mis malheureusement je ne comprends pas trop pourquoi celà ne marche pas :'(

voilà les paramètres :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
P1.x = 449.549011
P1.y = 362.789001
P2.x = 591.632996
P2.y = 374.233002
P3.x = 588.617981
P3.y = 60.318100
P4.x = 513.554016
P4.y = 53.292801
system[0] : 260
system[1] : -3606.67
system[2] : 400
system[3] : 0
system[4] : -865.444
system[5] : 253
system[6] : 0
system[7] : -6.77778

le résultat en image : http://imageshack.us/photo/my-images...enshottzl.png/

dont voici le 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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#include "cv.h"
#include "highgui.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <time.h>
#include <ctype.h>


double *mygetPerspectiveTransform(CvPoint2D32f P1, CvPoint2D32f P2, CvPoint2D32f P3, CvPoint2D32f P4){
 
        double *H;
	H = new double [18];
        double *adj=H+9;
 
        double sx = (P1.x-P2.x)+(P3.x-P4.x);
        double sy = (P1.y-P2.y)+(P3.y-P4.y);
        double dx1 = P2.x-P3.x;
        double dx2 = P4.x-P3.x;
        double dy1 = P2.y-P3.y;
        double dy2 = P4.y-P3.y;
 
        double z = (dx1*dy2)-(dy1*dx2);
        double g = ((sx*dy2)-(sy*dx2))/z;
        double h = ((sy*dx1)-(sx*dy1))/z;
 
        // matrice de transformation
        double a=H[0]=P2.x-P1.x+g*P2.x;
        double b=H[1]=P4.x-P1.x+h*P4.x;
        double c=H[2]=P1.x;
        double d=H[3]=P2.y-P1.y+g*P2.y;
        double e=H[4]=P4.y-P1.y+h*P4.y;
        double f=H[5]=P1.y;
        H[6]=g;
        H[7]=h;
        H[8]=1;
 
        // matrice de transformation inverse (matrice adjointe)
        adj[0]=e-f*h;
        adj[1]=c*h-b;
        adj[2]=b*f-c*e;
        adj[3]=f*g-d;
        adj[4]=a-c*g;
        adj[5]=c*d-a*f;
        adj[6]=d*h-e*g;
        adj[7]=b*g-a*h;
        adj[8]=a*e-b*d;
 
        return H;
}

double * invert(double u, double v, double* system) {
        double *H;
	H = new double [2];
	H[0] = (system[0]*u+system[1]*v+system[2])/(system[6]*u+system[7]*v+1); 
	H[1] = (system[3]*u+system[4]*v+system[5])/(system[6]*u+system[7]*v+1);
	return H;
}


int main(int argc, char* argv[]){

IplImage* image_src = cvLoadImage(argv[1],1);
IplImage* image_dst = cvCreateImage( cvGetSize(image_src), 8, 1 );

CvPoint2D32f P1, P2, P3, P4;

P1.x = 449.549;
P1.y = 362.789;

P2.x = 591.633;
P2.y = 374.233;

P3.x = 588.618;
P3.y = 60.3181;

P4.x = 513.554;
P4.y = 53.2928;

printf("P1.x = %f\n",P1.x);
printf("P1.y = %f\n",P1.y);

printf("P2.x = %f\n",P2.x);
printf("P2.y = %f\n",P2.y);

printf("P3.x = %f\n",P3.x);
printf("P3.y = %f\n",P3.y);

printf("P4.x = %f\n",P4.x);
printf("P4.y = %f\n",P4.y);

CvPoint2D32f P1T, P2T, P3T, P4T;
P1T.x = 400.0;
P1T.y = 253.0;
P2T.x = 660.0;
P2T.y = 253.0;
P3T.x = 510.0;
P3T.y = 106.0;
P4T.x = 555.0;
P4T.y = 106.0;

cvCircle( image_src, cvPointFrom32f(P1T), 9, CV_RGB(0,0,255), 3); //bleu
cvCircle( image_src, cvPointFrom32f(P2T), 9, CV_RGB(0,255,0), 3); //vert
cvCircle( image_src, cvPointFrom32f(P3T), 9, CV_RGB(255,0,0), 3); // rouge
cvCircle( image_src, cvPointFrom32f(P4T), 9, CV_RGB(255,255,0), 3); // jaune

double* system = mygetPerspectiveTransform(P1T,P2T,P3T,P4T);

printf("system[0] : %g\n",system[0]);
printf("system[1] : %g\n",system[1]);
printf("system[2] : %g\n",system[2]);
printf("system[3] : %g\n",system[3]);
printf("system[4] : %g\n",system[4]);
printf("system[5] : %g\n",system[5]);
printf("system[6] : %g\n",system[6]);
printf("system[7] : %g\n",system[7]);

        double width=image_src->width;
        double height=image_src->height;

        // pour chaque pixel de l'image cible   
        for(int y=0 ; y<image_src->width ; y++) {
                for(int x=0; x < image_src->height;x++) {
 
                        // conversion dans le repère orthonormé (u,v) [0,1]x[0,1]
                        double u = (double)x/width;
                        double v = (double)y/height;
 
			double* P = invert(u,v,system);

                        int sx=(int)round(P[0]);
                        int sy=(int)round(P[1]);
 
			//printf("SY = %d - SX = %d Y = %d - X = %d \n",sy,sx,y,x);
                       // pixel source 
			if(sy < image_src->width && sy >= 0 && sx < image_src->height && sx >= 0){
			// on évite que le get cherche un pixel qui est en dehors de l'image
                        CvScalar rgb=cvGet2D(image_src,sx,sy);
 			int dx=x;
                        int dy=y;
			cvSet2D(image_dst,dx,dy,rgb);
   			}

                }
        }

 cvNamedWindow("Source Image", 1 );
 cvNamedWindow("Destination Image", 1 );
	cvShowImage( "Source Image", image_src );
      cvShowImage( "Destination Image",image_dst );
 cvWaitKey();


return 0;
}

merci de votre aide
car le but étant qu'à la fin du projet je puisse connaitre la distance entre la caméra et un damier.
la calibration marche mais il me faut juste trouver la hauteur du damier redressé pour y appliquer thalès et trouver la distance que je cherche

[zaziechris]

Bonjour,
j'essai de passer de la perspective en projection orthogonale (vue de haut) sur matlab. savez vous s'il ya un moyen de le faire? car pour l'instant je n'arrive à rien.
Merci