Discussion :

[ijk-ref]

Si ton but est principalement de trouver l'équation permettant de corriger la perspective d'une image, tu devrais lire un peu ceci.

[ijk-ref]

Bonjour jp_developpeur, je préfère te répondre sur le forum

Désolé mais je me suis pas trop attarder sur ton code. C'est un langage que je n'ai jamais pratiquer et puis c'est pas très propre...

Par contre j'ai l'impression que tu ne travailles pas en 2D projective jusqu'au bout.

Je vois tes points fuites avec leurs coordonnées x,y Mais où elle la 3ème !?

Faut mieux tout garder ça simplifie grandement la représentation.

Par exemple l'un des 2 points de fuite M du rectangle A B C D se calcule simplement par :

M = (A^B)^(C^D) ( '^' : produit vectoriel )

[jp_developpeur]

En fait ma troisième coordonnée correpond à 1 si on a point fixe et 0 si on est à l'infini.

Mais tu l'as dejà codé en quel langage? ca pourrait m'intéresser si tu l'as fait en C.

En tout cas merci pour ta rep!

JP

[ijk-ref]

C#

Mon p'tit programme prend en entrée 4 points à placé aux quatre coins d'un rectangle dans images et l'affiche bien plate.

Mais bon j'ai juste fais ça pour m'amuser il y a quelque temps et c'est pas très propre.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
A----B
|    |
|    |
C----C

En tout cas le principe est :
- trouver le poids des 4 points A,B,C et D en s'aidant des points de fuite

puis après suffit de résoudre la matrice M qui permet d'avoir :
A.M = (0,0,1)
B.M = (1,0,1)
C.M = (0,1,1)

Tu bloques où ?

[jp_developpeur]

En fait c ce que j'ai fait mais en prenant un quadrilatère quelconque!

Dans mon programme (désolé de la propreté ms je dois me dépéché et je garde en commentaire ce que je me suis servie ou ce que je vais me servir), je prend 4 points. ensuite je détermine 2 points de fuite ((xPF1,yPF1) et (xPF2,yPF2)) gràce à l'annulation du déterminant (3 points alignés). De cela je détermine la matrice M_direct (matrice de passage d'un rectangle à une projection) j'en déduit la matrice inverse (par calcul: formule que j'ai trouvé sur net et que j'ai vérifié).
De cela, je parcours avec les indices (i,j) les pixels d'une image ('fictif') j'en déduis les coordonnées sur un carré ou rectangle (peu importe). Avec la matrice M de passage, j'en déduit les coordonnées sur l'image réel. Et de cette dernière j'en déduit (par équation linéaire) les indices correspondantes (i1,j1). J'obtiens donc la valeur du pixel de la 'vraie' image que je fais correpondre sur le carré ou rectangle, que je fais correspondre avc les indice (i,j).


Mon pb c'est que j'obtiens une image mais en fait je m'apperçoit que l'image réel n'est pas parcourue entièrement( avec les indice (i,j) (image redressée) j'obtiens pas tous les indices (i1,j1) de la 'vrai' image.


Pourtant les équations me paraissent juste.

J'ai vraiment besoin d'un coup de main,

Merci

JP

[ijk-ref]

Bon allons par étape :

1/ vérifier si le calcul de la matrice est exacte avec des cas simples.

Pour ma part il suffit juste de rentrer 3 points sources et 3 points destinations.

tu choisis de les placer où tu veux sur ton images sources :
[(x1,y1,1), (x2,y2,1), (x3,y3,1)]

la destination doit ressembler à cela
[(0, 0, 1), (I_MAX, 0, 1), (0, J_MAX, 1)]

La matrice trouvée doit vérifier qu'avec le calcul de

r = (0,0,1).M;

on ait bien x1 égal à r.X/r.W et y1 égal à r.Y/r.W

vérifier de même avec les 2 autres points.

2/ si tu arrives à cela, il suffira juste d'appliquer le bon poids aux points sources comme je l'explique dans mon 1er message posté ici.

En tout cas ta façon de faire ne peut pas être exacte dans tous les cas. Car un point de fuite peut très bien être en infini et toi tu négliges complètement la 3ème composante de la "2D projective"

[luxigo]

Je peux le faire en Java, mais en Matlab: non.

Cool ca marche pour passer du système de coordonées du plan de référence vers celui du plan de projection, merci 1000x

Par contre obtenir la fonction inverse, cad passer des coordonées de la projection vers celles du plan de référence;

je peux le faire en passant par une table,

mais en modifiant l'algorithme mathématique, non...

ben en fait si j'y suis arrivé par logique, voir pages suivantes du forum le code source pour opencv


J'imagine que vous devriez savoir

Merci d'avance !

Code c :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
#include <stdlib.h>
#include <math.h>
#include "ppm.h"
 
typedef struct Point {
        double x;
        double y;
 
} Point;
 
double *getSystem(Point *P);
void invert(double *u, double *v, double *system);
void image_rectify(struct ppm_struct *outppm,struct ppm_struct *ppm,Point P[4],double W, double H);
 
 
// vous pouvez zoomer in et out en multipliant W et H par le facteur de zoom souhaité
//  WxH > image de sortie = zoom in.
//  WxH < image de sortie = zoom out
 
void image_rectify(struct ppm_struct *outppm,struct ppm_struct *ppm,Point P[4],double W, double H) {
 
        double x;
        double y;
        double *S;
        struct color rgb;
 
        double zoom_factor=(outppm->width-W)/W;
        double factor=zoom_factor/2;
 
        // resolution du systeme
        S = getSystem(P);
 
        // pour chaque pixel de l'image cible
        for(y=-H*factor;y<H+H*factor;y++) {
                for(x=-W*factor;x<W+W*factor;x++) {
 
                        // conversion dans le repère orthonormé (u,v) [0,1]x[0,1]
                        double u = x/(W-1);
                        double v = y/(H-1);
 
                        // passage dans le repère perspective
                        invert(&u, &v, S);
 
                        // pixel (px,py) correspondant de l'image source
                        int px=(int)round(u);
                        int py=(int)round(v);
 
                        ppm_peek(ppm,px,py,&rgb);
                        ppm_plot(outppm,x+W*factor,y+H*factor,&rgb);
                }
        }
 
        free(S);
}
 
double *getSystem(Point *P) {
 
        double *system=calloc(8,sizeof(double));
        if (!system) {
                exit(1);
        }
 
        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;
}
 
void invert(double *u, double *v, double *system) {
        double x = (system[0]*u[0]+system[1]*v[0]+system[2])/(system[6]*u[0]+system[7]*v[0]+1);
        double y = (system[3]*u[0]+system[4]*v[0]+system[5])/(system[6]*u[0]+system[7]*v[0]+1);
        u[0]=x;
        v[0]=y;
}

on m'a demandé ppm.h, mais vous n'avez qu'à remplacer les types et les noms de procédures par ceux qui vous conviennent...

Vous trouverez aussi dans la page suivante du forum le code source pour opencv, qui gère aussi la transformation inverse, au cas où ..

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
#ifndef __ppm__ 
#define __ppm__ 
 
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/types.h>
#include <unistd.h>
#include <arpa/inet.h>
 
typedef int boolean;
 
typedef struct color {
        unsigned long red;
        unsigned long green;
        unsigned long blue;
} color;
 
typedef struct ppm_struct {
        unsigned char   *data;
        size_t          datasize;
        int             width;
        int             height;
        int             maxpixelvalue;
        int             pixelsize; // in bytes
        int             depth; // in bits
        unsigned char   *pixmap;
} ppm_struct;   
 
struct ppm_struct *ppm_new();
 
#ifdef __ppm
#define __extern
#else
#define __extern extern
#endif
 
__extern void ppm_dispose(struct ppm_struct *ppm);
 
__extern void  ppm_load(struct ppm_struct *ppm, char *ppm_filename);
__extern int ppm_parse(ppm_struct *p);
__extern void  ppm_save(struct ppm_struct *ppm,char *ppm_filename);
__extern void ppm_plot(struct ppm_struct *ppm, int x, int y, struct color *rgb);
__extern void ppm_peek(struct ppm_struct *ppm, int x, int y, struct color *rgb);
__extern int ppm_copy(struct ppm_struct *out, struct ppm_struct *in,boolean copydata);
 
 
#undef __extern 
 
#endif

et tant qu'à faire le code source de ppm.o .. sans garanties blabla gpl blabla free

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
 
#define __ppm
#include "ppm.h"
 
struct ppm_struct *ppm_new() {
        struct ppm_struct *ppm=(struct ppm_struct*)calloc(1,sizeof(struct ppm_struct));
        if (!ppm) {
                exit(1);
        }
        return ppm;
}
 
void ppm_dispose(struct ppm_struct *ppm) {
        if (!ppm) return;
        if (ppm->data) free(ppm->data);
        free(ppm);
}
 
void  ppm_load(struct ppm_struct *ppm,char *ppm_filename) {
 
        FILE *f;
 
        f=fopen(ppm_filename,"r");
        if (!f) {
                fprintf(stderr,"cannot open file: %s\n",ppm_filename);
                exit(1);
        }
 
        fseek(f,0,SEEK_END);
        ppm->datasize=ftell(f);
 
        fseek(f,0,SEEK_SET);
 
        ppm->data=(unsigned char*)malloc(ppm->datasize);
        if (!ppm->data) {
                exit(1);
        }
 
        if (fread(ppm->data,1,ppm->datasize,f)!=ppm->datasize) {
                fprintf(stderr,"ppm_load: read error\n");
                exit(1);
        }
 
        ppm_parse(ppm);
}
 
#define ppm_skip_comments(x) \
        while(1) { \
                if (!x) { \
                        fprintf(stderr,"file format not recognized\n"); \
                        exit(1);\
                } \
                if (x[1]!='#') break; \
                x=strchr(x+1,0xa); \
        }
 
int ppm_parse(ppm_struct *p) {
 
        char *c;
        char *d;
 
// file signature
        if (p->data[0]!='P' || p->data[1]!='6') {
                fprintf(stderr,"ppm signature not found\n");
                exit(1);
        }
 
        c=(unsigned char*)strchr(p->data,0xa);
        ppm_skip_comments(c);
 
// width
        d=strchr(c+1,' ');
        if (!d) {
                fprintf(stderr,"ppm file format not recognized\n");
                exit(1);
        }
        *d=0;
        p->width=atoi(c+1);
 
// height
        c=strchr(d+1,0xa);
        if ((!c)) {
                fprintf(stderr,"ppm file format not recognized\n");
                exit(1);
        }
        *c=0;
        p->height=atoi(d+1);
 
// max pixel value
        ppm_skip_comments(c);
        d=strchr(c+1,0xa);
        if ((!d)) {
                fprintf(stderr,"ppm file format not recognized\n");
                exit(1);
        }
        *d=0;
        p->maxpixelvalue=atoi(c+1);
 
        p->pixmap=d+1;
//
 
        p->depth=log2(p->maxpixelvalue+1);
        p->pixelsize=p->depth*3/8.0+0.5;
 
 
        return 0;       
}
 
void  ppm_save(struct ppm_struct *ppm,char *ppm_filename) {
 
        FILE *f;
        size_t total=0;
 
        f=fopen(ppm_filename,"w");
        if (!f) {
                fprintf(stderr,"cannot open file for writing: %s\n",ppm_filename);
                exit(1);
        }
 
        fprintf(f,"P6\xa%d %d\xa%d\xa", ppm->width, ppm->height, ppm->maxpixelvalue);
        while(1) {
                size_t count=fwrite(ppm->pixmap, ppm->pixelsize, ppm->width * ppm->height,f);
                if (!count) {
                        break;
                }
                total+=count;
                if (total==ppm->width*ppm->height*ppm->pixelsize) {
                        break;
                }
        }
 
        fclose(f);
 
}
 
 
void ppm_plot(struct ppm_struct *ppm, int x, int y, struct color *rgb) {
 
        long offset;
        offset=y*ppm->pixelsize*ppm->width+x*ppm->pixelsize;    
        if ( (offset<0) || (offset >= ppm->datasize-1) ) return;
 
        switch (ppm->depth) {
 
                case 8:
                        ((unsigned char*)ppm->pixmap)[offset]=(unsigned char)(rgb->red);        
                        ((unsigned char*)ppm->pixmap)[offset+1]=(unsigned char)(rgb->green);    
                        ((unsigned char*)ppm->pixmap)[offset+2]=(unsigned char)(rgb->blue);     
                        break;
 
                case 16:
                        ((unsigned short*)ppm->pixmap)[offset]=htons((unsigned short)(rgb->red));       
                        ((unsigned short*)ppm->pixmap)[offset+1]=htons((unsigned short)(rgb->green));   
                        ((unsigned short*)ppm->pixmap)[offset+2]=htons((unsigned short)(rgb->blue));    
                        break;
 
                case 24:
                case 32:
                        ((unsigned long*)ppm->pixmap)[offset]=htonl(rgb->red);  
                        ((unsigned long*)ppm->pixmap)[offset+1]=htonl(rgb->green);      
                        ((unsigned long*)ppm->pixmap)[offset+2]=htonl(rgb->blue);       
                        break;
 
                default:
                        fprintf(stderr,"ppm_plot: image depth %d not supported\n",ppm->depth);
                        exit(1);
 
        }
}
 
void ppm_peek(struct ppm_struct *ppm, int x, int y, struct color *rgb) {
 
        long offset;
 
        offset=y*ppm->pixelsize*ppm->width+x*ppm->pixelsize;    
        if ( (offset<0) || (offset >= ppm->datasize-1) ) return;
 
        switch (ppm->depth) {
 
                case 8:
                        rgb->red=((unsigned char*)ppm->pixmap)[offset];
                        rgb->green=((unsigned char*)ppm->pixmap)[offset+1];
                        rgb->blue=((unsigned char*)ppm->pixmap)[offset+2];
                        break;
 
                case 16:
                        rgb->red=ntohs(((unsigned short*)ppm->pixmap)[offset]);
                        rgb->green=ntohs(((unsigned short*)ppm->pixmap)[offset+1]);
                        rgb->blue=ntohs(((unsigned short*)ppm->pixmap)[offset+2]);
                        break;
 
                case 24:
                case 32:
                        rgb->red=ntohs(((unsigned long*)ppm->pixmap)[offset]);
                        rgb->green=ntohs(((unsigned long*)ppm->pixmap)[offset+1]);
                        rgb->blue=ntohs(((unsigned long*)ppm->pixmap)[offset+2]);
                        break;
 
                default:
                        fprintf(stderr,"ppm_plot: image depth %d not supported\n",ppm->depth);
                        exit(1);
 
        }
}
int ppm_copy(struct ppm_struct *out, struct ppm_struct *in,boolean copydata) {
 
        memcpy(out,in,sizeof(struct ppm_struct));
 
        if (copydata) {
                out->data=malloc(in->datasize);
        } else {
                out->data=calloc(1,in->datasize);
        }
 
        if (!out->data) {
                exit(1);
        }
 
        if (copydata) {
                memcpy(out->data,in->data,in->datasize);
        }
 
        out->pixmap=out->data+(in->pixmap-in->data);
        return 0;
}

[ijk-ref]

Par contre obtenir la fonction inverse, cad passer des coordonées de la projection vers celles du plan de référence; J'imagine que vous devriez savoir

Ca se résume à inverser une matrice 3x3

[luxigo]

Ca se résume à inverser une matrice 3x3

excellent merci j'ai deja trouvé le code ici pour ça

[luxigo]

Je peux le faire en Java, mais en Matlab: non.

encore merci à pseudocode

j'ai converti le code en C et j'ai implémenté un moyen de zoomer, comme si on changeait l'échelle du quadrilatère source par rapport à son centre.. ça peut servir..

[whiteboy]

bjr!si tu cherches a calculer l'homographie ,il éxiste une bibliotheque précompilé (c/c++):

http://www.ics.forth.gr/~lourakis/homest/

je pense que c est l'une des plus précise méthode pour éstimer l'homographie!

bon courage.

[luxigo]

oh bien sur il ya surement plein d'autres libs qui le font d'une façon ou d'une autre

Il me semble que l'algorithme posté par pseudocode donnait peut-être un meilleur résultat,

Il faudra que je teste pour en être sûr.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
/* 'from' et 'to' sont des pointeurs sur une liste de 4 points en float */
 
void warp(IplImage *src,CvPoint2D32F *from,IplImage *dst,CvPoint2D32F *to) {
        CvMat *H=cvCreateMat(3,3,CV_32FC1);
        cvGetPerspectiveTransform(from,to,H);
        cvWarpPerspective( src, dst, H, 0, cvScalarAll(0) );
        cvReleaseMat(&H);
}

[luxigo]

Ca se résume à inverser une matrice 3x3

Heu pas tout à fait, il manque une étape importante

En fait en relisant http://wcours.gel.ulaval.ca/2006/a/6...projective.pdf avec le code sous les yeux ca devient beaucoup plus clair

J'ai finalement réussi à deviner z= dans inverse_projective_mapping() pour pouvoir diviser x et y et faire fonctionner la transformation inverse

On peut utiliser la matrice adjointe au lieu de l'inverse de la matrice.

On peut aussi multiplier une matrice de transformation avec la matrice adjointe correspondant à un second quadrilatère pour obtenir le passage d'un quadrilatère à un autre via le carré unité en une seule opération.

Cette fois je crois que le cas est merci

Voilà pour opencv, c'est surement plus rapide mais y a pas d'interpolation.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
160
161
 
 
/*
 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 
 
 ( Pages 17-21 extraites de "Fundamentals of Texture Mapping and
   Image Warping", Paul Heckbert,
   Master’s thesis, UCB/CSD 89/516, CS Division, U.C. Berkeley, June 1989,
    <a href="http://www.cs.cmu.edu/~ph" target="_blank">http://www.cs.cmu.edu/~ph</a> )
 
 C'est Ephraim Cohen qui lui a montré le système 8x8
 "for projective mapping inference"
  au New York Institute of Technology Computer Graphics Lab.
  
 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;
                        }
                }
        }
}