Distance d'un point à un segment de droite

**defluc** · 26/04/2008, 23h42

Bonjour,

Pour trouver la distance d'un point à un segment de droite, j'utilise la fonction suivante :

Code :

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

1
2
3
4
5
6
7
8
9
Function DistPntSegm(XPnt, YPnt, Xb, Yb, Xe, Ye : Double) : Double;
Var
 TestX : Double;
begin
   Result := 1.7 * Power(10,308);
   TestX := PerpendFoot(XPnt, YPnt, Xb, Yb, Xe, Ye)[1];
   if PointOnSegment(Xb, Xe, TestX) then
      Result := DistPntLine(XPnt, YPnt, Xb, Yb, Xe, Ye);
end;

Dans laquelle je cherche les coordonnées du pied de la perpendiculaire abaissée de ce point sur la droite :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Function PerpendFoot(XPnt, YPnt, Xb, Yb, Xe, Ye : Double) : TPntArray;
Var
  A, B, C : Double;
begin
  Result[1] := 0;
  Result[2] := 0;
  A := Xb - Xe;
  if A = 0 then exit;
  A := (Yb - Ye)/A; // b = begin e = ending
  B := (Xb*Ye-Xe*Yb)/(Xb-Xe);
  C := (YPnt-A*XPnt-B)/(1+Sqr(A));
  Result[1] := XPnt + A*C;  //  Coordonnée X
  Result[2] := YPnt - C;      //  Coordonnée Y 
end;

. Puis, je vérifie si ce point est bien sur la droite :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
Function PointOnSegment(Xb, Xe, XPnt : Double) : Boolean;
Var
  A: Double;
begin
  Result := False;
  if A = 0 then exit;
  A:= (Xe-Xb);
  if A > 0 then A := (XPnt-Xb)/A;
  if Not((A<0) or (A > 1)) then Result := True
                           else result := False;
end;

. Et enfin, si c'est le cas, je calcule la distance du point à la droite :

Code :

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

1
2
3
4
5
6
7
8
9
10
Function DistPntLine(XPnt, YPnt, Xb, Yb, Xe, Ye : Double) : Double;
Var
  A, B, C, L : Double;
begin
  A := Yb - Ye;
  B := Xe - Xb;
  C := (Xb*Ye)-(Xe*Yb);
  L := Sqrt(Sqr(Xb-Xe)+Sqr(Yb-Ye));
  Result := Abs(((A*XPnt)+(B*YPnt)+ C)/L);
end;

Les résultats ne sont pas satisfaisants pour des raisons que j'attribue aux comparaisons de valeurs en virgule flottante.
Cela m'aiderait si quelqu'un pouvait me proposer un algo donnant satisfaction à 100%.

Merci d'avance

**pseudocode** · 27/04/2008, 00h40

Légère adaptation du code de la classe QuadEdge de l'implémentation de la Triangulation de Delaunay

Code java :

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
 
/**
 * Return the distance from a point to a segment
 *
 * @param ps,pe the start/end of the segment
 * @param p the given point
 * @return the distance from the given point to the segment
 */
private static double distanceToSegment(Point ps, Point pe, Point p) {
 
	if (ps.x==pe.x && ps.y==pe.y) return distance(ps,p);
 
	int sx=pe.x-ps.x;
	int sy=pe.y-ps.y;
 
	int ux=p.x-ps.x;
	int uy=p.y-ps.y;
 
	int dp=sx*ux+sy*uy;
	if (dp<0) return distance(ps,p);
 
	int sn2 = sx*sx+sy*sy;
	if (dp>sn2) return distance(pe,p);
 
	double ah2 = dp*dp / sn2;
	int un2=ux*ux+uy*uy;
	return Math.sqrt(un2-ah2);
}
 
/**
 * return the distance between two points
 *
 * @param p1,p2 the two points
 * @return dist the distance
 */
private static double distance(Point p1, Point p2) {
	int d2 = (p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y);
	return Math.sqrt(d2);
}

**Zavonen** · 27/04/2008, 09h08

On peut aussi faire ça ainsi:
On utilise une notion: celle de puissance d'un point par rapport à une droite:
Si ux+vy+w=0 est l'équation de la droite D
et a et b sont les coordonnées du point M
Puissance(M,D) = ua+vb+w par définition.
Cette notion a deux intérêts:
Elle permet de distinguer les deux demi-plans séparé par la droite en:
P(M,D) >0 et P(M,D)<0
La distance d'un point à une droite s'exprime aussi au moyen de la puissance:
d(M,D)=p(M,D)/sqrt(u^2+v^2)
La méthode est alors la suivante.
Si [BE] est le segment, on calcule les équations des deux perpendiculaires à [BE]
passant par B et E, c'est immédiat, elle sont du type ux+vy+w=0 où u et v sont les coordonnées du vecteur BE, le w s'obtenant en traduisant l'appartenance de B (resp E) aux droites en question.
Appelons D1 et D2 ces perpendiculaires.
Tout se résume maintenant à savoir la position de M par rapport à ces deux droites, pour cela on calcule:
P(M,D1) P(M,D2) P(B,D2) P(E,D1)
Il suffit ensuite de comparer les signes de ces quantités pour savoir si M est à l'intérieur de la bande délimitée par D1 et D2 où bien à l'expérieur du côté D1 ou à l'extérieur côté D2.
PS: Ceci est plus facile à FAIRE qu'à DIRE
Avantages: Le calcul des puissances ne nécessite ni division ni racines carrées.
Il n'est pas nécessaire de connaitre la projection de M sur (BE), autrement dit la discrimination se fait uniquement sur des sommes et des produits. Bien entendu quand le cas est reconnu il faut calculer une racine et faire une division pas moyen d'y échapper quelle que soit la méthode.

**Alikendarfen** · 27/04/2008, 11h20

Pratiquement la même que Pseudocode, en fait... sauf pour la fin (je n'ai pas compris la projection ?

)

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
 
        /// <summary>
        /// Retourne la distance du point P au segment AB.
        /// Si A et B sont confondus, retourne la distance de A a P
        /// Si la projection perpendiculaire du point sur la droite contenant AB est hors de AB, retourne double.PositiveInfinity.
        /// Sinon retourne la distance du point à la droite contenant AB.
        /// </summary>
        /// <param name="P"></param>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <returns></returns>
        static double DistancePaAB( Point P, Point A, Point B )
        {
            Point AB = new Point() { X = B.X - A.X, Y = B.Y - A.Y }; // Vecteur AB
            double scalaireABAB = Scalaire( AB, AB );
 
            Point AP = new Point() { X = P.X - A.X, Y = P.Y - A.Y }; // Vecteur AP
 
            // A et B sont confondus
            if ( scalaireABAB <= double.Epsilon )
            {
                return Norme( AP );
            }
            else
            {
                double scalaireAPAB = Scalaire( AP, AB );
 
                // Si le produit scalaire est négatif : P est 'avant' A dans le sens AB.
                // S'il est plus grand que le scalaire AB AB : P est 'apres' B dans le sens AB.
                // Dans les deux cas, la projection de P sur la droite contenant AB est hors de AB
                if ( (scalaireAPAB < 0d) || (scalaireAPAB > scalaireABAB) )
                {
                    return double.PositiveInfinity;
                }
                // Sinon la distance est la composante scalaire de AP sur le vecteur perpendiculaire à AB
                else
                {
                    Point perpAB = new Point() { X = AB.Y, Y = -AB.X }; // Un vecteur perpendiculaire à AB
 
                    return Math.Abs( Scalaire( AP, perpAB ) / Math.Sqrt( scalaireABAB ) );
                }
            }
        }
 
        static double Norme( Point V )
        {
            return Math.Sqrt( (V.X * V.X) + (V.Y * V.Y) );
        }
 
        static double Scalaire( Point V1, Point V2 )
        {
            return (V1.X * V2.X) + (V1.Y * V2.Y);
        }

**Zavonen** · 27/04/2008, 12h18

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
# -*- coding: cp1252 -*-
import math
 
# calcule la distance du point M de coordonnées xM,yM, au segment BE
# B coordonnées xB,yB  E cordonnes xE,yE
# en utilisant la puissance d'un point par rapport à une droite
 
def distance(xB,yB,xE,yE,xM,yM):
    # coordonnées a,b du vecteur EB
    a=xE-xB
    b=yE-yB
    # équation de la perpendiculaire D1 en B à (EB): ax+by+w1
    w1=-a*xB-b*yB
    # équation de la perpendiculaire D2 en E à (EB): ax+by+w2
    w2=-a*xE-b*yE
    # équation de la droite (EB) : bx-ay+w3
    w3= a*yB-b*xB
    #puissance de M par rapport à D1
    PMD1=a*xM+b*yM+w1
    # puissance de M par rapport à D2
    PMD2=a*xM+b*yM+w2
    #puissance de B par rapport à D2
    PBD2=a*xB+b*yB+w2
    #puissance de E par rapport à D1
    PED1=a*xE+b*yE+w1
    # A ce stade encore ni racine ni quotient
    if PMD1*PED1 <0: #M et E de part et d'autre de D1
        return math.sqrt((xM-xB)*(xM-xB)+(yM-yB)*(yM-yB)) # pas de quotient
    if PMD2*PBD2 <0: #M et B de part et d'autre de D2
        return math.sqrt((xM-xE)*(xM-xE)+(yM-yE)*(yM-yE)) # idem
    # sinon
    return abs(b*xM-a*yM+w3)/math.sqrt(a*a+b*b) # là rien à faire
 
 
def main():
    xB=yB=1
    xE=yE=2
    xM=yM=0
    print distance(xB,yB,xE,yE,xM,yM)
    xM=1
    print distance(xB,yB,xE,yE,xM,yM)
    xM=2
    print distance(xB,yB,xE,yE,xM,yM)
    xM=3
    print distance(xB,yB,xE,yE,xM,yM)
if __name__ == '__main__':
    main()

**pseudocode** · 27/04/2008, 14h35

Envoyé par Alikendarfen

Pratiquement la même que Pseudocode, en fait... sauf pour la fin (je n'ai pas compris la projection ?

)

Tu as raison, c'était pas très clair. J'ai modifié le code pour que ca soit plus compréhensible.

==> |AH|² = (AB.AC)² / |AB|²

Pythagore : AC²=AH²+HC²

==> distance(H,C) = Racine( AC² - AH² )

**souviron34** · 27/04/2008, 18h59

je re-reprends mes citations :

http://www.developpez.net/forums/sho...7&postcount=10

Cette fois c'est le point 1.02

**defluc** · 27/04/2008, 19h05

Je vous remercie pour vos sources.
J'ai essayé de les transposer en Delphi mais les résultats ne correspondent pas à mes attentes. Peut-être parce que ma "traduction" est erronée.

**pseudocode** · 27/04/2008, 20h41

Envoyé par defluc

Je vous remercie pour vos sources.
J'ai essayé de les transposer en Delphi mais les résultats ne correspondent pas à mes attentes. Peut-être parce que ma "traduction" est erronée.

L'autre solution c'est qu'on se soit tous plantés... mais franchement, ca ne serait pas de bol.

**defluc** · 28/04/2008, 17h43

Voici ma source en Delphi qui «laisse passer» des cas qui ne le devraient pas :

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
Function DistPntToSegment(Xp, Yp, Xb, Yb, Xe, Ye : Double) : Double;
Var
  ScalaireABAB, ScalaireAPAB, ABx, ABy, APx, APy : Double;
begin
  ABX := Xe - Xb;
  ABY := Ye - Yb;
  APX := Xp - Xb;
  APY := Yp - Yb;
  if (Xe = Xb) and (Ye = Yb) then
  begin   //  A et B son confondus
    Result := Abs(Sqrt(Sqr(Xp-Xb)+Sqr(Yp-Ye)));;
    exit;
  end;
  ScalaireABAB := (ABX * ABX) + (ABY * ABY);
  ScalaireAPAB := (APX * ABX) + (APY * ABY);
  if (ScalaireAPAB < 0) And (ScalaireAPAB > ScalaireABAB) then
  begin   // La projection de P sur la droite se fait en dehors du segment
    Result := 1.7 * Power(10,308);
    exit;
  end;
  Result := Abs(((APX * ABY) + (APY * -ABX))/Sqrt(ScalaireABAB));
end;

**Zavonen** · 28/04/2008, 18h24

Code :

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

1
2
3
4
  begin   // La projection de P sur la droite se fait en dehors du segment
    Result := 1.7 * Power(10,308);
    exit;
  end;

Et alors ? Est ce une raison pour retourner une valeur astronomique ???

Result := 1.7 * Power(10,308);

As tu seulement lu nos nombreuses réponses ?

**Alikendarfen** · 28/04/2008, 19h57

Envoyé par defluc

Code :

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

1
2
3
4
5
6
Function DistPntToSegment(Xp, Yp, Xb, Yb, Xe, Ye : Double) : Double;
Var
....
 
  if (ScalaireAPAB < 0) And (ScalaireAPAB > ScalaireABAB) then
...

Tu as mis "And" au lieu de "Or"

**defluc** · 29/04/2008, 14h54

Pour Alikendarfen : le résultat reste erroné même après la modification

if (ScalaireAPAB < 0) And (ScalaireAPAB > ScalaireABAB) then

Il doit donc y avoir une autre erreur dans ma source.

Pour Zavonen : Bon d'accord, «1.7 * Power(10,308)», c'était rapport à la notion d'infini que je pourrais effectivement limiter au maximum probable que mon application ne devrait pas rencontrer comme distance.

Pour ton deuxième commentaire, j'ai lu toutes les réponses mais les traductions en Delphi ne me donnent pas encore de résultats probants.

Concernant ta méthode, voici ce que j'ai codé pour ce qui est de

Si [BE] est le segment, on calcule les équations des deux perpendiculaires à [BE] passant par B et E

Code :

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

1
2
3
4
5
6
7
8
9
10
11
Function PntToSegmByPower(Xp, Yp, Xb, Yb, Xe, Ye : Double) : Double;
Var
  a, b, c, cp, Perp1, Perp2 : Double;
begin
  a := Ye - Xb;        // Coordonnées de l'équation paramétrique du vecteur sont (a;b)
  b := Xb - Xe;
  c := Xe*Yb - Xb-Ye;
  cp := b*Xb - a*Yb; 
  Perp1 := -b*Xb + a*Yb + c;  // Le vecteur perpendiculaire mené de (Xb;Yb) est (-b*Xb;a*Yb)
  Perp2 := -b*Xe + a*Ye + c;  // Le vecteur perpendiculaire mené de (Xe;Ye) est (-b*Xe;a*Ye)
end;

Mais pour la suite, je coince

**souviron34** · 29/04/2008, 15h23

on t'as dit qu'avec les équations paramétriques ce serait plus dur..

Dans la référence que j''ai pointé, tu as directement le paramètre r qui te donne la solution...

C'est pourtant pas compliqué

**pseudocode** · 29/04/2008, 17h00

Envoyé par defluc

Pour Zavonen : Bon d'accord, «1.7 * Power(10,308)», c'était rapport à la notion d'infini que je pourrais effectivement limiter au maximum probable que mon application ne devrait pas rencontrer comme distance.

Ce n'était pas le sens de la remarque de Zavonen.

ce n'est pas parce que la projection du point ne tombe pas entre les extrémités du segment que la distance est forcément infinie !!

Code :

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

1
2
3
4
5
6
7
8
9
 
^
|
| P  A        B
| X  X--------X
|
+---------------->
 
Distance(P,[AB]) = 3

**Zavonen** · 29/04/2008, 19h02

Ce n'était pas le sens de la remarque de Zavonen.

Merci Pseudocode, il y en a au moins qui suit ...
Cela dit, je pense qu'on peut tout reprendre, je veux dire n'importe laquelle des méthodes présentées ici dans cette discussion pour résoudre l'autre problème:
Distance d'un point à un polygone CONVEXE; C'est à mon avis le minimum des distances aux côtés considérés comme des segments.
Preuve:
Soit T le point du polygone réalisant le minimum de distance avec M. Ce point existe (pour des raisons de compacité) mais n'est pas forcément unique.
Le segment [TM] contient forcément un point C du contour du polygone et la distance MC est plus petite que la distance MT donc la distance de M au polygone est égale à la distance au pourtour. Cette distance est atteinte en un point C toujours pourra raison de compacité, (non forcément unique à nouveau), C appartient à un côté MnMn+1 il est impossible de trouver sur ce côté un point plus proche de M, cela contredirait la définition même de T. Donc la distance de M au polygone est égale à la distance de M à ce côté.
D'où mon affirmation.
Maintenant il n'y a plus qu'à itérer une des méthodes données sur les côtés consécutifs du polygone.

**Alikendarfen** · 29/04/2008, 23h44

Envoyé par pseudocode

Ce n'était pas le sens de la remarque de Zavonen.

ce n'est pas parce que la projection du point ne tombe pas entre les extrémités du segment que la distance est forcément infinie !!

Code :

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

1
2
3
4
5
6
7
8
9
 
^
|
| P  A        B
| X  X--------X
|
+---------------->
 
Distance(P,[AB]) = 3

Ca ce n'est qu'une question de convention, non ? Tout dépend de l'énoncé et de ce qu'on considère comme valide comme distance point / segment.

**Zavonen** · 30/04/2008, 10h17

Non, il n'y a ici aucune question de convention.
Distance d'un point M à un ensemble K dans un espace métrique E c'est par définition:
Inf d(M,P) pour P dans K.
Ce Inf existe toujours car les nombres d(M,P) sont minorés par 0 (théorème de la borne inférieure). En outre, comme un segment est un compact (fermé borné) et que la fonction P --> d(M,P) est continue le inf est effectivement atteint ( Il existe un P dans K tel que d(M,P)=Inf), cependant ce point peut ne pas être unique.

**defluc** · 30/04/2008, 13h06

Et bien Zavonen, j'ai suivi ton conseil et j'ai encodé ce que le lien que tu m'as gentiment renseigné propose.

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
Function TMapForm.DistPntSegm(XPnt, YPnt, Xb, Yb, Xe, Ye : Double) : Double;
Var
 Len, ParamR, ParamS, X, Y : Double;
begin                               
  Result := 999999;
  Len := Dist2Pnts(Xb, Yb, Xe, Ye);
  ParamR := ((Xpnt - Xb)*(Xe - Xb) + (YPnt - Yb) + (Ye - Yb))/Sqr(Len);
  if (ParamR > 0) and (ParamR < 1) then
  begin
    X := Xb + ParamR * (Xe - Xb); // X point projection de Xpnt sur le segment
    Y := Yb + ParamR * (Ye - Yb); // Y point projection de Ypnt sur le segment
    ParamS := ((Yb - Ypnt)*(Xe - Xb) - (Xb - XPnt) + (Ye - Yb))/Sqr(Len);
    Result := Abs(ParamS) * Len;
  end;
end;