Projection d'une portion de sphère sur un rectangle

[vibi27]

Bonjour,

Voici mon problème:
-on prend 4 points A,B,C,D sur la terre dont on connait les latitudes/longitudes.
-on prend un 5ie point M à l'intérieur de ce "quadrilatère sphérique" (je ne sais pas si c'est le bon terme), dont on connait aussi la latitude/longitude
-ce quadrilatère est projeté sur un rectangle (disons mon écran d'ordinateur !), c'est à dire que dans ce repère, A est le coin bas gauche de mon écran, B le point bas droite de l'écran, C le point haut droite et D le point haut gauche. Je connais donc dans ce repère écran (R') les coordonnées des points A,B,C,D projetés, que j'appelle A',B',C',D'.
-on a, en considérant l'origine en bas à gauche, l'axe X horizontal vers la droite et l'axe Y vers le haut:
A'=(0,0), B'=(w, 0), C'=(w,h), D'=(0,h) (w et h sont des constantes connues)

Mon problème est de calculer les coordonnées M', dans le repère écran R', qui est le projeté du point M.

J'ai eu l'idée suivante: le rapport des distances orthodromiques dans R devrait être égal au rapport des distances cartésiennes dans R'
d(MA)/d(MB) = d(M'A')/d(M'B')
d(MB)/d(MC) = d(M'B')/d(M'C')
etc...

J'ai déjà la formule pour calculer la distance ortho = d(XY) quelque soit X et Y

Dans R' (le rectangle), on peut facilement trouver des triangles rectangles et donc des formules faisant intervenir les distances d(M'A'), d(M'B'), d(M'C') et d(M'D'), et ensuite calculer leur rapport ex: d(M'A')/d(M'B').
Mais j'obtiens ainsi des équations à résoudre de type y*y = ax*x + bx + c, et je ne m'en sors pas

Voici mes questions:
1) Pensez vous que la conservation du rapport des distances ortho soit une bonne idée ?
2) Y a t-il un calcul plus "formel"/mathématique qui pourrait convenir.

PS1: OK, la terre n'est pas une sphère, mais plutôt un ellipsoïde, mais on va simplifier un peu.
PS2: j'ai pensé aussi aux translations/rotations/homothéties (homographies), mais je pense que ce n'est pas approprié
PS3: si la formulation de mon problème n'est pas suffisante, n'hésitez pas à m'en faire part je tenterai d'être plus clair.
PS4: une approximation me convien
PS4: En finale, ma formule, si j'en trouve une, j'en ferai un algorithme dans un langage informatique (Java)

Merci d'avance pour vos réponses
-----

[Resartus]

On ne peut pas conserver le même rapport des distances entre la sphère (ou l'ellipsoide) et sa projection car cela fait un système de trois equations pour seulement deux variables. Il existe quantité de projections possibles. Elle sont en général conformes, c'est à dire qu'elles conservent les angles, mais aucune ne peut conserver toutes les distances. Il faut accepter des compromis

Pour une aire de dimension pas trop grande, le système de projection habituellement utilisé par les cartographes est celui de Lambert (utilisé par exemple pour les cartes de la France). Elle déforme toutes les distances (horizontales et verticales), mais c'est celle qui donne la déformation moyenne la plus faible

L'autre méthode assez standard est la projection UTM (Universal Transverse Mercator). Elle conserve les distances horizontales, conserve aussi les distances verticales près du méridien central, mais les déforme de plus en plus quand on s'en éloigne latéralement.

Vous trouverez sans problème sur internet toutes les équations correspondant à ces deux systèmes

[invite57a1e779]

on prend 4 points A,B,C,D sur la terre dont on connait les latitudes/longitudes.

par exemple :

A = Madrid
B = Le Caire
C = Fontainebleau
D = Paris

e quadrilatère est projeté sur un rectangle (disons mon écran d'ordinateur !), c'est à dire que dans ce repère, A est le coin bas gauche de mon écran, B le point bas droite de l'écran, C le point haut droite et D le point haut gauche

Quelles sont les projections capables de réaliser cet exploit ?

[Dynamix]

Salut

ce quadrilatère est projeté sur un rectangle

Tu projettes sur un plan , pas sur un rectangle .
Sinon ce n' est pas une projection et il faut définir ce que c' est .

[A voir en vidéo sur Futura]

[Dlzlogic]

Bonsoir,
On suppose que la dimension géographique à représenter est d'une dimension pas plus grand que la France. Sinon, il y a une petite opération supplémentaire à faire.
On a donc 4 points connus dans 2 systèmes différents, le système géographique et le système "écran".
La transformation qui permet de passe de l'un à l'autre est une transformation affine. Son équation s'écrit
X=TX + XX.x + XY.y
Y=TY + YX.x + YY.y
où TX, TY, XX, XY, YX, YY sont des paramètres à calculer.
On dispose de 4 points, soit 8 équations. donc des équations en sur-nombre. La solution la "plus probable", c'est à dire la meilleure est celle qui minimise la somme des carrés des écarts. On écrit cela. Cette somme sera minimum pour les valeurs qui annulent les dérivées partielles. On obtient donc un système linéaire de 6 équations à 6 inconnues, facile à résoudre.
Je peux vous donner le module qui calcule cela en C et/ou en PHP.
Ce calcul n'est pas difficile à faire. Il est vivement conseillé de résoudre le système par la méthode du pivot de Gauss (cad pas d'inversion de matrice).

[vibi27]

Merci à tous pour vos réponses constructives et très rapides.
En fait je ne cherchais pas la projection, c'est un partenaire qui l'a produite et je n'ai pas accès à son code source.
J'ai oublié de préciser que la zone A,B,C,D était de l'ordre de 1Km à 10km de "côté" environ.
La dernière réponse me plait bien, et si en plus j'avais un exemple de code en C pour la résolution ce serait génial , merci d'avance Dlzlogic.
Je vous envoie mon mail en MP (à moins qu'il ne soit déjà visible sur mon profil).

[vibi27]

Envoi de MP impossible pour des raisons techniques, aussi voici mon mail: *** Adresse mail ***

[Médiat]

Bonjour,

La charte de ce forum précise bien :

3. Ne sont pas autorisés les pseudos qui contiennent une adresse mail, messagerie instantanée ou web. Il vous est fortement déconseillé également d'utiliser ce que vous employez par ailleurs comme identifiant de sécurité (par exemple identifiant servant à vous connecter à votre FAI, à votre banque...). Il est interdit de mettre votre adresse email sur le forum.

Cordialement

Médiat

[leon1789]

Bonsoir

Bonsoir,
On suppose que la dimension géographique à représenter est d'une dimension pas plus grand que la France. Sinon, il y a une petite opération supplémentaire à faire.
On a donc 4 points connus dans 2 systèmes différents, le système géographique et le système "écran".
La transformation qui permet de passe de l'un à l'autre est une transformation affine.

Il y a plein d'autres transformations qui permettent d'envoyer 4 points sur 4 points. Mais on peut chercher une transformation affine... si elle existe !

Son équation s'écrit
X=TX + XX.x + XY.y
Y=TY + YX.x + YY.y
où TX, TY, XX, XY, YX, YY sont des paramètres à calculer.
On dispose de 4 points, soit 8 équations. donc des équations en sur-nombre.

...donc ta transformation affine n'existe pas dans l'absolu . CQFD

La solution la "plus probable",

pas de bol, il n'y a pas de probabilité dans l'histoire

c'est à dire la meilleure ...

la meilleure pour quel critère ?

...est celle qui minimise la somme des carrés des écarts. On écrit cela. Cette somme sera minimum pour les valeurs qui annulent les dérivées partielles. On obtient donc un système linéaire de 6 équations à 6 inconnues, facile à résoudre.

ok

Je peux vous donner le module qui calcule cela en C et/ou en PHP.
Ce calcul n'est pas difficile à faire. Il est vivement conseillé de résoudre le système par la méthode du pivot de Gauss (cad pas d'inversion de matrice).

Manque de chance, c'est que la méthode de Gauss ne donnera pas de solution, puisqu'il n'y a pas de solution exacte (en général).
La solution suivant les moindres carrés n'a rien à voir avec le pivot de Gauss

Bref, toujours autant de contre-vérités, ça ne change pas.

[leon1789]

Manque de chance, c'est que la méthode de Gauss ne donnera pas de solution, puisqu'il n'y a pas de solution exacte (en général).
La solution suivant les moindres carrés n'a rien à voir avec le pivot de Gauss

Hum... Rectification : la méthode de Gauss permet en effet de déterminer la solution pour les moindres carrés (en supposant qu'elle soit unique, ce qui est généralement le cas) car les dérivées partielles sont des polynômes de degré 1 (exceptionnellement) dans le contexte précisé par Dlzlogic.

[vibi27]

Bonsoir,
j'ai implémenté la méthode du pivot de Gauss (google: pivot de gauss java -> 1ere réponse) appliquée aux coefficients des fonctions affines suggérées par dlzlogic.
Et je n'ai pas utilisé le moindre carré (humour).
Cela fonctionne suffisamment pour me rendre compte que mes données d'entrée sont parfois farfelues (quadrilatère concave ou croisé alors que j'attendais un convexe, lat/lon fixes alors que l'image bouge...).

En tout cas merci pour vos réponses.

[leon1789]

C'est dommage si vos données sont parfois farfelues (comme vous dites).

Cela étant, la méthode proposée par Dlzlogic n'a rien d'être LA méthode, la plus probable, la meilleure, etc.
On pourrait tout bonnement changer la fonction recherchée en prenant autant de paramètres que de coordonnées de points, à savoir 8 (puisque 4 points dans un plan) :
f(x,y) = (a + b.x + c.y + d.xy, a' + b'.x + c'.y + d'.xy)
Pour envoyer 4 points sur 4 points de manière exacte (sauf cas particulier), il suffit alors de résoudre un système linéaire 8 x 8.

On pourrait aussi utiliser des homographies (liée à la géométrie projective) :
f(x,y) = ( (a + b.x + c.y)/(1 + d.x + e.y), (a' + b'.x + c'.y)/(1 + d'.x + e'.y) )
fonction qui est connue après résolution d'un système linéaire 10 x 10.

Tout cela, sans approximation, et autre détournement "frauduleux"

Etc. Il faut bien comprendre le contexte du problème pour en trouver la "vraie" solution (si elle existe).

[leon1789]

On pourrait aussi utiliser des homographies (liée à la géométrie projective) :
f(x,y) = ( (a + b.x + c.y)/(1 + d.x + e.y), (a' + b'.x + c'.y)/(1 + d'.x + e'.y) )

arf, écrit trop vite :
f(x,y) = ( (a + b.x + c.y)/(1 + d.x + e.y), (a' + b'.x + c'.y)/(1 + d.x + e.y) )
fonction qui est connue après résolution d'un système linéaire 8x8

Si vous cherchez à conserver le birapport, c'est peut-être ces fonctions dont vous avez besoin...
Voir : https://fr.wikipedia.org/wiki/Application_projective

[vibi27]

Excellent les formules, merci Léon.
Effectivement mon problème ressemble fortement à trouver une application projective mais je ne vois pas comment écrire les équations f(x,y)
avec 8 inconnues:
Si je reprends la formule analytique de wikipedia (concernant l'application projective) avec les notations suivantes:
n=2
a indice i,j = aij
x indice 1 = x c'est ma longitude
x indice 2 = y c'est ma latitude
x' indice 1 = x'
x' indice 2 = y'

J'obtiens:
x'=(a11.x + a12.y + a13)/(a31.x + a32.y + a33)
y'=(a21.x + a22.y + a23)/(a31.x + a32.y + a33)

Soit 9 inconnues a11,a12,a13, a21,a22,a23, a31,a32,a33.
J'ai dû me tromper quelquepart...?

[leon1789]

Soit 9 inconnues a11,a12,a13, a21,a22,a23, a31,a32,a33.
J'ai dû me tromper quelquepart...?

Non, la remarque est très juste, c'est bien 9 paramètres inconnus a priori. Mais regardons de plus près :

x'=(a11.x + a12.y + a13)/(a31.x + a32.y + a33)
y'=(a21.x + a22.y + a23)/(a31.x + a32.y + a33)

Si a33 est nul (cas exceptionnel) alors il reste 8 inconnues à déterminer ;
si a33 n'est pas nul, et alors on peut simplifier les fractions en divisant les numérateurs et dénominateurs par a33 et il vient alors
x'=(a'11.x + a'12.y + a'13)/(a'31.x + a'32.y + 1)
y'=(a'21.x + a'22.y + a'23)/(a'31.x + a'32.y + 1)
ce qui donne encore une fois 8 inconnues à déterminer.

[vibi27]

Merci,

OK j'ai pigé. Et comme je ne sais pas à l'avance si a33 vaut zéro je ferai les deux résolutions a33=0 et a33<>0 et je prendrai celle qui donne le moins d'erreurs de calculs.
Ces deux équations définissent une application non linéaire, à cause de la division, aussi je me vois mal faire le pivot de gauss dessus.
Avez vous une idée de la technique que je pourrais employer pour résoudre ces équations non linéaires,
je suis preneur de toute idée, y compris une méthode non formelle par exemple un nom d'algorithme informatique, car en final c'est dans une babasse que cela finira

[leon1789]

OK j'ai pigé. Et comme je ne sais pas à l'avance si a33 vaut zéro je ferai les deux résolutions a33=0 et a33<>0 et je prendrai celle qui donne le moins d'erreurs de calculs.

il n'y a aura pas d'erreur de calcul (en comparaison au approximation par les moindres carrés) : la fonction enverra exactement vos 4 premiers points sur les 4 points d'arrivée.
Par ailleurs, vous constaterez que a33<>0 (donc a33=1) est le cas le plus souvent rencontré (pour ne pas dire davantage... le fait de mettre a_33=0 provoquera la nullité des autres a_ij dans la résolution du système linéaire !)

Ces deux équations définissent une application non linéaire, à cause de la division, aussi je me vois mal faire le pivot de gauss dessus.

x'=(a11.x + a12.y + a13)/(a31.x + a32.y + a33)
en multipliant par le dénominateur :
(a31.x + a32.y + a33) . x'= (a11.x + a12.y + a13)
voilà, c'est linéaire en les inconnus.
(a31.x + a32.y + a33) . x' - (a11.x + a12.y + a13) = 0

[ansset]

bonjour,
je n'ai pas suivi vos calculs,
mais au vu des distances considérée ( 1 à 10 km ), on peut avec une erreur négligeable considérer la terre "plate" sur cette zone et en déduire l'équation du plan en question en (x,y,z). non?
Cdt

[jiherve]

Bonjour,
en effet sur de courtes distances on peut se passer de projection compliquées ce que font nombre de logiciel de génération d'images cartographiques.
La limite est cependant vite atteinte.
Ayant déjà eu a résoudre ce genre de problème il faut passer par la méthode des moindres carrés et Cholesky pour l'inversion , il faudra tout de même disposer d'assez de quadruplets (L,G,X,Y), nombre qui dépend du degré du polynôme, pour avoir un espoir de solution fiable .
Ceci dit inutile de réinventer l'eau tiède il faut s'accrocher à un système de projection connu.
JR

[Dlzlogic]

Bonjour Ansset et Jiherve,

*** Procès d'intention ***
il s'agit, non pas de calage de plans, mais de redressement de photo aérienne.
Dans le cas de calage de plan, on identifie un certain nombre de couples de points homogènes, et on choisi soit une méthode au hasard comme certains, soit la méthode des moindres carrés, comme font les autres.
Dans le cas de redressement, la transformation n'est naturellement pas une transformation affine qui conserve les parallèles, mais une transformation du type perspective.
Dans le cas général, c'est à dire photographie de surface non plane, le redressement nécessite 2 prises vue. La résolution est basée sur la stéréoscopie et cela fait partie de la photogrammétrie.
Dans le cas où l'objet est plan (façade ou terrain plan) les calculs sont plus simples et se ramènent à utiliser deux fois Thales.
Concernant le cas général, on peut utiliser une seule photo, si on connait l'altimétrie du terrain. Mais là on atteint un domaine que je domine mal. Mais dans tous les cas, il y aura toujours de points de calage en sur-nombre (disons, au moins 6) donc, toujours un calage par la méthode habituelle, après redressement.

[vibi27]

Merci pour les réponses.
Pour le moment je continue sur la piste de leon1789

Ci dessous j'essaye de construire la matrice:

Soit
(I) a'31.x.x' + a'32.y.x' + x' -a'11.x -a'12.y -a'13 = 0
(II) a'31.x.y' + a'32.y.y' + y' -a'21.x -a'22.y -a'23 = 0

J'insère un "1" dans le vecteur GaussCoef ci dessous, ce "1" sera multiplié par x' ou y' pour donner le 3ie terme de (I) ou (II).
GaussCoef est donc un vecteur 9x1 qui vaut: {a'31, a'32, 1, a'11, a'12, a'13, a'21, a'22, a'23}

MAT est une Matrice 9x9 remplie avec:
-les longitudes xi connues (x1, x2, x3, x4)
-les latitudes yi connues (y1, y2, y3, y4)
-les xi' pixel connus (x1', x2', x3', x4')
-les yi' pixel connus (y1', y2', y3', y4')

MAT = {
{ x1.x1', y1.x1', x1', -x1, -y1, -1, 0, 0, 0},
{ x1.x1', y1.x1', y1', 0, 0, 0, -x1, -y1, -1},
{ x2.x2', y2.x2', x2', -x2, -y2, -1, 0, 0, 0},
{ x2.x2', y2.x2', y2', 0, 0, 0, -x2, -y2, -1},
{ x3.x3', y3.x3', x3', -x3, -y3, -1, 0, 0, 0},
{ x3.x3', y3.x3', y3', 0, 0, 0, -x3, -y3, -1},
{ x4.x4', y4.x4', x4', -x4, -y4, -1, 0, 0, 0},
{ x4.x4', y4.x4', y4', 0, 0, 0, -x4, -y4, -1},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0},
};


// Cote droit des equations: vecteur 9x1
vecteurZero = { 0, 0, 0, 0, 0, 0, 0, 0, 0}


// Resolution par pivot de Gauss (ou Gauss-Jordan ?) de l'equation:
// MAT . gaussCoef = vecteurZero
GaussCoef = lsolve(MAT, vecteurZero);

// Application des coefficients trouvés pour trouver la valeur pixel de x'
// x'=(a'11.x + a'12.y + a'13)/(a'31.x + a'32.y + 1)
x'=(gaussCoef[3].x + gaussCoef[4].y + gaussCoef[5])/(gaussCoef[0].x + gaussCoef[1].y + 1)

// Application des coefficients trouvés pour trouver la valeur pixel de y'
// y'=(a'21.x + a'22.y + a'23)/(a'31.x + a'32.y + 1)
y'=(gaussCoef[6].x + gaussCoef[7].y + gaussCoef[8])/(gaussCoef[0].x + gaussCoef[1].y + 1)

Que pensez vous de cette façon de procéder ?

[leon1789]

Il faut répéter 4 fois les deux équations suivantes,

(I) a'31.x.x' + a'32.y.x' -a'11.x -a'12.y - a'13 = -x'
(II) a'31.x.y' + a'32.y.y' -a'21.x -a'22.y - a'23 = -y'

avec la constante (qui n'est pas en facteur d'un des a_ij) dans le second membre (qui n'est pas 0,0...,0).
On obtient une matrice 8x8 , pas 9x9.

MAT = {
{ x1.x1', y1.x1', -x1, -y1, -1, 0, 0, 0},
{ x1.x1', y1.x1', 0, 0, 0, -x1, -y1, -1},
{ x2.x2', y2.x2', -x2, -y2, -1, 0, 0, 0},
{ x2.x2', y2.x2', 0, 0, 0, -x2, -y2, -1},
{ x3.x3', y3.x3', -x3, -y3, -1, 0, 0, 0},
{ x3.x3', y3.x3', 0, 0, 0, -x3, -y3, -1},
{ x4.x4', y4.x4', -x4, -y4, -1, 0, 0, 0},
{ x4.x4', y4.x4', 0, 0, 0, -x4, -y4, -1}
};

vecteur = {-x1', -y1', -x2', -y2', -x3', -y3', -x4', -y4'}

GaussCoef = lsolve(MAT, vecteur);

[leon1789]

Dans le cas de calage de plan, on identifie un certain nombre de couples de points homogènes, et on choisi soit une méthode au hasard comme certains, soit la méthode des moindres carrés, comme font les autres.

...la méthode des moindres carrés, ok, mais appliquée à quel modèle mathématique ? Un modèle << au hasard >> ?
Certains savants disent qu'il faut une transformation affine, d'autres disent qu'il faut utiliser et ...

Dans le cas de redressement, la transformation n'est naturellement pas une transformation affine qui conserve les parallèles, mais une transformation du type perspective.
Dans le cas général, c'est à dire photographie de surface non plane, le redressement nécessite 2 prises vue. La résolution est basée sur la stéréoscopie et cela fait partie de la photogrammétrie.
Dans le cas où l'objet est plan (façade ou terrain plan) les calculs sont plus simples et se ramènent à utiliser deux fois Thales.

et pourtant Thalès conserve le parallélisme... alors ?

[vibi27]

Argh je n'étais pas très loin, mais c'est bien sur, enfin c'est très clair avec vôtre explication.
Thx Leon

[leon1789]

attention, erreur de coeff sur les lignes paires de la matrice MAT :

MAT = {
{ x1.x1', y1.x1', -x1, -y1, -1, 0, 0, 0},
{ x1.y1', y1.y1', 0, 0, 0, -x1, -y1, -1},
{ x2.x2', y2.x2', -x2, -y2, -1, 0, 0, 0},
{ x2.y2', y2.y2', 0, 0, 0, -x2, -y2, -1},
{ x3.x3', y3.x3', -x3, -y3, -1, 0, 0, 0},
{ x3.y3', y3.y3', 0, 0, 0, -x3, -y3, -1},
{ x4.x4', y4.x4', -x4, -y4, -1, 0, 0, 0},
{ x4.y4', y4.y4', 0, 0, 0, -x4, -y4, -1}
};

[jiherve]

Bonjour,
La methode dont j'ai fais état était à l'origine utilisée pour traiter des photo satellites , elle peut servir à d'autre chose car basiquement c'est une approximation polynomiale .Il y eu même des circuits intégrés dédiés à ces calculs TMC2301/TMC2302 de chez feu TRW
Tout le problème revient à avoir assez de quadruplets liant l'espace de départ à l'espace d'arrivée, en carto ils sont calculables et l'obtention des polynômes de déformation ne fait que simplifier les calculs, pour l'imagerie (anamorphose) il faut définir les points par d'autre méthode , s'il s’agit de photo aériennes on peut le faire par identification de points remarquables sur la photo et une carte.
JR

[vibi27]

Yes, c'est exact, j'avais reproduit mon erreur du post dans mon code et j'obtenais des y'=x', mes cercles étaient tous des segments
Maintenant cela fonctionne, je vois bien mes objets lat/lon projetés au bon endroit dans mon repère écran.
Attention l'axe Y peut être orienté vers le bas dans certaines librairies de dessin, et donc à l'envers du repère de math classique.
Merci pour tout le temps passé à m'aider