calcul de la distance minimale entre un point et un ellipsoïde

**kamelcompte** · 09/04/2010, 17h50

Bonjour tout le monde
je cherche comment calculer la distance minimale séparant un point (x, y, z) et une ellipsoide((x0,y0,z0), a, b, c)
(x0,y0,zO) est le centre de l'ellipsoide et a, b et c sont ses 3 rayons.
dans le cas dune sphère on peut calculer cette distance comme suit:
d=sqrt[(x-x0)^2+(y-y0)^2+(z-z0)^2)] - rayon.
je cherche comment je peut calculer cette distance dans le cas d'une ellipsoide
merci d'avance

**ToTo13** · 10/04/2010, 12h05

Bonjour,

je pense que c'est une question pour les matheux du forum Algorithmes/Mathématques.

**Jean Dumoncel** · 10/04/2010, 12h47

Envoyé par ToTo13

je pense que c'est une question pour les matheux du forum Algorithmes/Mathématques.

Exact, du coup je déplace la discussion.

En ce qui concerne la question, c'est loin d'être aussi trivial que pour une sphère.

Après tu peux peut-être simplifier le problème (tout dépend de la précision recherchée) en prenant la distance minimum entre ton point M et les 2 foyers de l'ellipsoïde (F ou F') et en recherchant le point d'intersection entre le segment MF (ou MF') et l'ellipsoïde. Mais la distance obtenue ne sera pas la distance minimale entre ton point et l'ellipsoïde...

**plegat** · 10/04/2010, 13h59

Envoyé par magelan

En ce qui concerne la question, c'est loin d'être aussi trivial que pour une sphère.

Peut-être pas trivial, mais pas hyper compliqué non plus...

Tu calcules k²=(x-x0)^2/a²+(y-y0)^2/b²+(z-z0)^2/c² (qui te donne le carré du rapport d'homothétie entre ton ellipsoïde et celui qui passerait par ton point), et tu calcules ta distance d(x,ellipsoïde)=d(x,x0)*(k-1)/k

Non? (bien sûr, en faisant attention au cas particulier k=0, soit (x,y,z)=(x0,y0,z0))

**kamelcompte** · 12/04/2010, 00h59

merci beaucoups pour les réponses
j'aimerai recevoir un lien d'un article ou d'un cours pour bien étudier le principe et la démonstration mathématique de la méthode:

d(x,ellipsoïde)=d(x,x0)*(k-1)/k
avec
k²=(x-x0)^2/a²+(y-y0)^2/b²+(z-z0)^2/c²

merci d'avance

**Zavonen** · 12/04/2010, 01h37

Un ellipsoïde E étant convexe, Si N est un point extérieur la distance de N à l'ellipsoïde, est la distance NH où H est la projection orthogonale de M sur E.
On part donc d'une représentation paramétrique de l'ellipsoïde M(s,t). En tout point M(s,t) on peut calculer, par dérivation partielle, les vecteurs directeurs u(s,t), v(s,t) du plan tangent à E en M. Il suffit d'ajuster s et t pour que les deux vecteurs u(s,t) et v(s,t) soient orthogonaux à NM(s,t). On tombe sur un système de deux équations à deux inconnues, mais il n'est pas du premier degré.
Voir par exemple:
http://translate.google.fr/translate...blog/%3Fp%3D86

**Jean Dumoncel** · 12/04/2010, 10h05

Envoyé par plegat

Tu calcules k²=(x-x0)^2/a²+(y-y0)^2/b²+(z-z0)^2/c² (qui te donne le carré du rapport d'homothétie entre ton ellipsoïde et celui qui passerait par ton point), et tu calcules ta distance d(x,ellipsoïde)=d(x,x0)*(k-1)/k

Justement ce n'est pas aussi simple... tes équations ne calculent pas la distance minimale. Pour l'illustrer, j'ai écrit un code matlab avec tes équations pour une ellipse :

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
 
clear
close all
% Paramètre de l'ellipse
a=4;
b=1;
xO = 1;
yO = 1;
ang = 0;
 
% Point du plan
xM = 5;
yM = 3;
 
% Affichage de l'ellipse et du point
figure,
t = 0:.1:2*pi+0.1;
plot(a*cos(t) + xO,b*sin(t) + yO);
hold on,
plot([xO xM],[yO yM],'-*');
text(xO,yO,'O');
text(xM,yM,'M');
 
% Calcul de la distance
k2=(xM-xO)^2/a^2+(yM-yO)^2/b^2;
d=sqrt((xM-xO)^2+(yM-yO)^2)*(sqrt(k2)-1)/sqrt(k2);
 
% Affichage du cercle de centre (xM,yM) et de rayon d
plot(d*cos(t) + xM,d*sin(t) + yM);
 
axis equal,axis([-4 8 -1 6]),grid;

Et voici ce que l'on obtiens :

Le résultat est le cercle de centre M et de rayon la distance obtenue. On vois bien que ce n'est pas la distance minimale.
Une meilleure approche serait de considérer les foyers plutôt que le centre comme je l'avais suggéré dans mon premier message, mais ce n'est toujours pas la distance minimale.

Pour obtenir une solution beaucoup plus précise, je suis d'accord avec Zanoven, il faut utiliser un algorithme de minimisation, mais ce n'est toujours pas une solution analytique.

**plegat** · 12/04/2010, 10h29

Envoyé par magelan

Le résultat est le cercle de centre M et de rayon la distance obtenue. On vois bien que ce n'est pas la distance minimale.

Plus que clair!

Donc ne pas suivre ce que j'ai dit!

**pseudocode** · 12/04/2010, 11h00

On doit pouvoir se ramener au cas du cercle/sphère par une transformation linéaire (changement de repère). non ?

**Zavonen** · 12/04/2010, 11h21

On doit pouvoir se ramener au cas du cercle/sphère par une transformation linéaire (changement de repère). non ?

Ben,.. peut-être mais en tout cas on peut se ramener à un problème plan.
Dans le lien que j'ai donné l'auteur travaille avec une représentation paramétrique issue d'une équation cartésienne. Mauvaise pioche! car il traîne des radicaux.
Commençons par traiter le problème dans le plan avec un paramétrage de l'ellipse, par des coordonnées polaires:
x=acos(t)
y=bsin(t)
Le vecteur tangent V(t) au point M(t) de paramètre t s'obtient tout se suite de même que le vecteur NM(t). On annule le produit scalaire NM(t).V(t) pour obtenir une équation en cos(t) sin(t). On la résoud comme d'hab. en passant par la tangente de l'arc-moitié.
Cela vu, si l'ellipsoïde est un ellipsoïde de révolution d'axe Oz, comme il est globalement inchangé par toute rotation d'axe Oz la distance de N à l'ellipsoïde sera la distance de r(N) à l'ellipsoïde où r désigne une rotation quelconque d'axe Oz. Il suffit donc d'amener le point N dans un plan de coordonnées par une telle rotation et d'appliquer la solution du problème plan.

**pseudocode** · 12/04/2010, 12h51

Envoyé par Zavonen

Commençons par traiter le problème dans le plan avec un paramétrage de l'ellipse, par des coordonnées polaires:
x=acos(t)
y=bsin(t)
Le vecteur tangent V(t) au point M(t) de paramètre t s'obtient tout se suite de même que le vecteur NM(t). On annule le produit scalaire NM(t).V(t) pour obtenir une équation en cos(t) sin(t). On la résoud comme d'hab. en passant par la tangente de l'arc-moitié.

C'est la méthode des multiplicateurs de Lagrange, non ?

**kamelcompte** · 12/04/2010, 14h25

Envoyé par Zavonen

Un ellipsoïde E étant convexe, Si N est un point extérieur la distance de N à l'ellipsoïde, est la distance NH où H est la projection orthogonale de M sur E.
On part donc d'une représentation paramétrique de l'ellipsoïde M(s,t). En tout point M(s,t) on peut calculer, par dérivation partielle, les vecteurs directeurs u(s,t), v(s,t) du plan tangent à E en M. Il suffit d'ajuster s et t pour que les deux vecteurs u(s,t) et v(s,t) soient orthogonaux à NM(s,t). On tombe sur un système de deux équations à deux inconnues, mais il n'est pas du premier degré.
Voir par exemple:
http://translate.google.fr/translate...blog/%3Fp%3D86

prière de revoiur le lein, parcequ'il ne fonctionne pas

**kamelcompte** · 12/04/2010, 14h29

Envoyé par magelan

Justement ce n'est pas aussi simple... tes équations ne calculent pas la distance minimale. Pour l'illustrer, j'ai écrit un code matlab avec tes équations pour une ellipse :

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
 
clear
close all
% Paramètre de l'ellipse
a=4;
b=1;
xO = 1;
yO = 1;
ang = 0;
 
% Point du plan
xM = 5;
yM = 3;
 
% Affichage de l'ellipse et du point
figure,
t = 0:.1:2*pi+0.1;
plot(a*cos(t) + xO,b*sin(t) + yO);
hold on,
plot([xO xM],[yO yM],'-*');
text(xO,yO,'O');
text(xM,yM,'M');
 
% Calcul de la distance
k2=(xM-xO)^2/a^2+(yM-yO)^2/b^2;
d=sqrt((xM-xO)^2+(yM-yO)^2)*(sqrt(k2)-1)/sqrt(k2);
 
% Affichage du cercle de centre (xM,yM) et de rayon d
plot(d*cos(t) + xM,d*sin(t) + yM);
 
axis equal,axis([-4 8 -1 6]),grid;

Et voici ce que l'on obtiens :

Le résultat est le cercle de centre M et de rayon la distance obtenue. On vois bien que ce n'est pas la distance minimale.
Une meilleure approche serait de considérer les foyers plutôt que le centre comme je l'avais suggéré dans mon premier message, mais ce n'est toujours pas la distance minimale.

Pour obtenir une solution beaucoup plus précise, je suis d'accord avec Zanoven, il faut utiliser un algorithme de minimisation, mais ce n'est toujours pas une solution analytique.

merci
le plus important c'est que la solution est unique pour chaque point de l'espace

**Jean Dumoncel** · 12/04/2010, 14h51

Envoyé par kamelcompte

prière de revoiur le lein, parcequ'il ne fonctionne pas

C'est peut-être google translate qui n'a pas réussi à charger la page. Le lien du document en anglais : http://cheshirekow.com/blog/?p=86

Envoyé par kamelcompte

le plus important c'est que la solution est unique pour chaque point de l'espace

Dans ce cas la solution de plegat (celle que j'ai codé) devrait te convenir.

Envoyé par Zavonen

Commençons par traiter le problème dans le plan avec un paramétrage de l'ellipse, par des coordonnées polaires:
x=acos(t)
y=bsin(t)
Le vecteur tangent V(t) au point M(t) de paramètre t s'obtient tout se suite de même que le vecteur NM(t). On annule le produit scalaire NM(t).V(t) pour obtenir une équation en cos(t) sin(t). On la résoud comme d'hab. en passant par la tangente de l'arc-moitié.

J'ai essayé de développer les équations et j'arrive avec une expression dont le numérateur est un polynôme de degré 4. J'ai cherché les solutions avec les outils de MATLAB, et il m'a trouvé 4 racines complexes sur mon exemple. Bon je suis peut-être aller un peu vite, je vais vérifier les coefficients.

**Nemerle** · 12/04/2010, 16h23

Je me trompe peut-être, mais je ne vois pas où est le problème.

Si l'ellipse a pour milieu de ces foyers O, a pour rayons a et b et pour excentricité e, alors son équation polaire (R,o) dans son repère naturel centré en O est donnée par R^2=b^2/(1-e^2*cos^2(o)).

Si ton point M a pour coordonnée (rho, theta) dans le repère naturel de l'ellipse, le point que tu cherches est (R, theta) où R^2=b^2/(1-e^2*cos^2(theta)).

C'est en somme l'intersection de l'ellipse avec la demi-droite OM.

**pseudocode** · 12/04/2010, 16h32

Envoyé par Nemerle

C'est en somme l'intersection de l'ellipse avec la demi-droite OM.

Comme l'a dit Zavonen, il faut que R soit le projeté orthogonal de M. et donc que RM soit perpendiculaire à la courbure ( = colineaire a la normale ).

Hors, je ne vois pas en quoi OM serait colineaire a la normale (sauf aux extremités de l'ellipse).

**Nemerle** · 12/04/2010, 16h50

tu vois, je ne voyais pas où était le problème, et j'avais tord, pensant trop "cercle".

Evidemmment mon point proposé n'est pas correct: un truc utile est que la normale à la tangente d'un point de l'ellipse est la bissectrice de son angle avec les focales. On devrait pouvoir proposer une solution polaire en utilisant ça, je vais réfléchir d'ici demain.

**Jean Dumoncel** · 12/04/2010, 17h00

A y est, j'ai corrigé mes coefficients

Dans le cas d'une ellipse avec un angle de 0 et centré sur l'origine, on trouve rapidement le point cherché (enfin il faut quand même un algo qui donne toutes les solutions d'une équation de degré 4, en l'occurence sous matlab, il suffit d'utiliser la fonction roots

).

Voici le code :

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
clear
close all
 
% Parametre de l'ellipse
a=4;
b=1;
xO = 0;
yO = 0;
ang = 0;
 
% Point du plan
xM = 5;
yM = 3;
 
% xM = -4;
% yM = 6;
 
% Affichage de l'ellipse et du point
figure,
t = 0:.1:2*pi+0.1;
plot(a*cos(t) + xO,b*sin(t) + yO);
hold on,
text(xO,yO,'O');
text(xM,yM,'M');
 
axis equal,axis([-5 9 -2 7]),grid;
 
% Calul des coefficients du polynôme :
c4 = -b*yM;
c3 = -2*a*xM-2*a^2+2*b^2;
c2 = 0;
c1 = -2*a*xM+2*a^2-2*b^2;
c0 = b*yM;
 
% Calcul des solutions :
p = roots([c4 c3 c2 c1 c0]);
for k=1:length(p)
    if isreal(p(k))
        ang = 2*atan2(p(k),1);
        x = a*cos(ang);
        y = b*sin(ang);
        plot(x,y,'*r');
        % Affichage des solutions :
        d = sqrt((xM-x)^2+(yM-y)^2);
        plot(d*cos(t) + xM,d*sin(t) + yM,'g');
        plot([x xM],[y yM],'-');
    end
end

Le résultat :

On obtient 2 solutions qui sont perpendiculaires par rapport à la tangente de l'ellipse.
Il suffit alors de garder la plus petite distance.

**Zavonen** · 12/04/2010, 17h09

J'ai essayé de développer les équations et j'arrive avec une expression dont le numérateur est un polynôme de degré 4.

C'est exact !
Pour une implémentation, voir:
http://www.google.fr/url?sa=t&source...6isw-dv504dTtA

C'est la méthode des multiplicateurs de Lagrange, non ?

Ah, ben si Giuseppe Lodovico Lagrangia se met à pomper mes idées maintenant, ...

**Zavonen** · 12/04/2010, 17h28

Cerises sur le gâteau:
Existence:
Si N est un point de l'espace la fonction M ---> d(M,N) où N est sur l'ellipsoïde est une fonction continue sur un compact elle atteint donc effectivement son minimum en (au moins) un point.
Unicité:
Résulte de la convexité.Si D est une droite avec deux points M1 et M2 distincts il existe un point M3 strictement compris entre M1 et M2 tel que d(N,M3) < Inf(d(N,M1),d(M,n2)).
Ramener 3D a 2D:
Si l'ellipsoïde E est d'axe Oz et si M est le point de l'ellipsoïde où le minimum de distance est atteint, si S est la sphère de centre N et de rayon MN, alors la figure formée par l'ellipsoïde et la sphère est symétrique par rapport au plan défini par Oz et le point M, le plan tangent en M à la sphère et le plan tangent en M à l'ellipsoïde coïncident. Le cas où M est sur Oz est trivial.
Les deux solutions exhibées par Magelan correspondent, comme on le voit au maxi et au mini, car il s'agit bien de la méthode de Lagrange.

calcul de la distance minimale entre un point et un ellipsoïde

Mathématiques

Discussions similaires

Partager

Partager