2009-02-12 3 views
7

Dato due segmenti di linea 2D, A e B, come si calcola la lunghezza del segmento di linea 2D più breve, C, che collega A e B?Collegamento di due segmenti di linea

+0

Definire connettersi. Intendi connettere le loro estremità o trovare il segmento di linea più breve tra qualsiasi punto sulle linee? – strager

+0

@strager, nella geometria euclidea, sono paralleli o i punti finali sono più vicini, quindi si controllano i vettori A1-B1, A1-B2, A2-B1 e A2-B2. – paxdiablo

+0

2D o 3D? Il caso 2D è quasi banale, mentre il caso 3D richiede algoritmi più complessi. – fbonnet

risposta

6

consideri vostra linea due segmenti a e B di essere rappresentato da due punti ciascuno (in Fortran!):

linea a rappresentata dalla A1 (x, y), A2 (x, y)

Linea B rappresentata da B1 (x, y) B2 (x, y)

Prima verifica se le due linee si intersecano utilizzando questo algoritmo.

Se intersecano, la distanza tra le due linee è zero e il segmento di linea che le unisce è il punto di intersezione.

Se non si intersecano, utilizzare questo metodo: http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/ per calcolare la distanza più breve tra:

  1. punto A1 e la linea B
  2. Point A2 e la linea B
  3. Point B1 e la linea Un B2
  4. punto e linea A

La più breve di questi quattro segmenti di linea è la tua risposta.

+0

Inoltre: primo controllo per A e B che si incrociano (A1 + vettore (A1-> A2) * a = B1 + vettore (B1-> B2) * b con numeri reali a e b! = 0). Secondo controllo per A e B che sono paralleli (cioè vettore (A1-> A2) * a = vettore (B1-> B2) con un numero reale! = 0). –

+0

Non credo che questa risposta sia corretta. per uno come sopra se le linee attraversano la distanza più breve tra le linee è zero. tuttavia la distanza più breve da qualsiasi punto di arrivo verso l'altra linea> 0. –

+0

Hai ragione entrambi: è necessario verificare l'intersezione. Tuttavia, penso che non sia necessario verificare se le linee sono parallele. – Alterlife

4

This page ha informazioni che si sta cercando.

+0

La mia testa è asplode! – spoulson

+0

una risposta per 3d funzionerà per 2d, 2d solo un caso speciale per 3d dove z sempre == 0. quindi nel codice sudo nella parte inferiore z (i) == z (j) == 0 –

2

punta rapida: se si desidera confrontare le distanze sulla base di punti, non è necessario fare le radici quadrate.

E.g. per vedere se P-to-Q è una distanza minore di Q-a-R, basta controllare (pseudocodice):

square(P.x-Q.x) + square(P.y-Q.y) < square(Q.x-R.x) + square(Q.y-R.y) 
0

This page ha una bella descrizione breve per trovare la distanza minima tra due linee, anche se @ Strager di collegamento comprende un codice

+0

Questo è 3D, e si riferisce a linee, non a segmenti di linea. Non rilevante qui. –

+0

La domanda originale non specificava 2D. –

+0

OK. Scuse. Suggerisci a chiunque abbia downvoted di rimuovere il downvote. (Non l'ho fatto.) –

2

Afterlife ha detto: "Prima controlla se le due linee si intersecano usando questo algoritmo", ma non ha indicato quale algoritmo intendesse. Ovviamente, è l'intersezione della linea segmenti non le linee estese che contano; eventuali segmenti di linea non paralleli (esclusi i punti coincidenti che non definiscono una linea) si intersecheranno, ma la distanza tra i segmenti di linea non sarà necessariamente zero. Quindi presumo che intendesse "segmenti di linea" piuttosto che "linee" lì.

Il collegamento Afterlife fornito è un approccio molto elegante per trovare il punto più vicino su una linea (o segmento di linea o raggio) su un altro punto arbitrario. Questo funziona per trovare la distanza da ciascun endpoint all'altro segmento di linea (vincolando il parametro calcolato a non meno di 0 per un segmento di linea o raggio e non più di 1 per un segmento di linea), ma non gestire la possibilità che un punto interno su un segmento di linea sia più vicino di entrambi gli endpoint perché si intersecano effettivamente, quindi è richiesto il controllo extra sull'intersezione.

Per quanto riguarda l'algoritmo per la determinazione dell'intersezione del segmento di linea, un approccio sarebbe quello di trovare l'intersezione delle linee estese (se parallelo, quindi hai finito), e quindi determinare se quel punto si trova all'interno di entrambi i segmenti, tali come portando il punto-prodotto dei vettori dal punto di intersezione, T, ai due estremi:

((Tx - A1x) * (Tx - A2x)) + ((Ty - A1y) ​​* (Ty - A2y))

Se questo è negativo (o "zero"), T è compreso tra A1 e A2 (o in corrispondenza di un punto finale). Controllare in modo simile per l'altro segmento di linea. Se uno dei due era maggiore di "zero", i segmenti di linea non si intersecano. Naturalmente, questo dipende dal trovare prima l'intersezione delle linee estese, che può richiedere l'espressione di ogni linea come un'equazione e la risoluzione del sistema mediante riduzione gaussiana (ecc.).

Ma ci può essere un modo più diretto senza dover risolvere per il punto di intersezione, prendendo il prodotto incrociato dei vettori (B1-A1) e (B2-A1) e il prodotto incrociato dei vettori (B1 -A2) e (B2-A2). Se questi prodotti trasversali sono nella stessa direzione, A1 e A2 si trovano sullo stesso lato della linea B; se sono in direzioni opposte, si trovano sui lati opposti della linea B (e se 0, quindi uno o entrambi sono su linea B). Allo stesso modo controllare i prodotti incrociati di vettori (A1-B1) e (A2-B1) e di (A1-B2) e (A2-B2). Se uno qualsiasi di questi prodotti incrociati è "zero" o se i punti finali di entrambi i segmenti di linea cadono sui lati opposti dell'altra linea, i segmenti di linea stessi devono intersecarsi, altrimenti non si intersecano.

Ovviamente, è necessaria una comoda formula per calcolare un prodotto incrociato di due vettori dalle loro coordinate. O se tu potessi determinare gli angoli (essendo positivi o negativi), non avresti bisogno dell'attuale cross-product, poiché è la direzione degli angoli tra i vettori che ci interessano veramente (o il seno dell'angolo, in realtà) . Ma io che la formula per il cross-prodotto (in 2-D) è semplicemente:

trasversale (V1, V2) = (V1x * V2Y) - (V2x * V1Y)

Questa è l'asse z componente del vettore 3D del prodotto incrociato (dove i componenti x e y devono essere zero, poiché i vettori iniziali sono nel piano z = 0), quindi puoi semplicemente guardare il segno (o "zero").

Quindi, è possibile utilizzare uno di questi due metodi per verificare l'intersezione del segmento di linea nell'algoritmo Afterlife descrive (facendo riferimento al collegamento).

3

Utilizzando l'idea generale degli algoritmi Afterlife's e Rob Parker's qui sopra, ecco una versione C++ di un insieme di metodi per ottenere la distanza minima tra 2 segmenti 2D arbitrari. Questo gestirà segmenti sovrapposti, segmenti paralleli, segmenti intersecanti e non intersecanti. Inoltre, utilizza vari valori epsilon per proteggere dall'imprecisione in virgola mobile. Infine, oltre a restituire la distanza minima, questo algoritmo ti darà il punto sul segmento 1 più vicino al segmento 2 (che è anche il punto di intersezione se i segmenti si intersecano). Sarebbe piuttosto banale per tornare anche il punto su [P3, P4] più vicini a [P1, P2] se lo si desidera, ma lascio che come esercizio per il lettore :)

// minimum distance (squared) between vertices, i.e. minimum segment length (squared) 
#define EPSILON_MIN_VERTEX_DISTANCE_SQUARED 0.00000001 

// An arbitrary tiny epsilon. If you use float instead of double, you'll probably want to change this to something like 1E-7f 
#define EPSILON_TINY 1.0E-14 

// Arbitrary general epsilon. Useful for places where you need more "slop" than EPSILON_TINY (which is most places). 
// If you use float instead of double, you'll likely want to change this to something like 1.192092896E-04 
#define EPSILON_GENERAL 1.192092896E-07 

bool AreValuesEqual(double val1, double val2, double tolerance) 
{ 
    if (val1 >= (val2 - tolerance) && val1 <= (val2 + tolerance)) 
    { 
     return true; 
    } 

    return false; 
} 


double PointToPointDistanceSquared(double p1x, double p1y, double p2x, double p2y) 
{ 
    double dx = p2x - p1x; 
    double dy = p2y - p1y; 
    return (dx * dx) + (dy * dy); 
} 


double PointSegmentDistanceSquared(double px, double py, 
            double p1x, double p1y, 
            double p2x, double p2y, 
            double& t, 
            double& qx, double& qy) 
{ 
    double dx = p2x - p1x; 
    double dy = p2y - p1y; 
    double dp1x = px - p1x; 
    double dp1y = py - p1y; 
    const double segLenSquared = (dx * dx) + (dy * dy); 
    if (AreValuesEqual(segLenSquared, 0.0, EPSILON_MIN_VERTEX_DISTANCE_SQUARED)) 
    { 
     // segment is a point. 
     qx = p1x; 
     qy = p1y; 
     t = 0.0; 
     return ((dp1x * dp1x) + (dp1y * dp1y)); 
    } 
    else 
    { 
     t = ((dp1x * dx) + (dp1y * dy))/segLenSquared; 
     if (t <= EPSILON_TINY) 
     { 
      // intersects at or to the "left" of first segment vertex (p1x, p1y). If t is approximately 0.0, then 
      // intersection is at p1. If t is less than that, then there is no intersection (i.e. p is not within 
      // the 'bounds' of the segment) 
      if (t >= -EPSILON_TINY) 
      { 
       // intersects at 1st segment vertex 
       t = 0.0; 
      } 
      // set our 'intersection' point to p1. 
      qx = p1x; 
      qy = p1y; 
      // Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if 
      // we were doing PointLineDistanceSquared, then qx would be (p1x + (t * dx)) and qy would be (p1y + (t * dy)). 
     } 
     else if (t >= (1.0 - EPSILON_TINY)) 
     { 
      // intersects at or to the "right" of second segment vertex (p2x, p2y). If t is approximately 1.0, then 
      // intersection is at p2. If t is greater than that, then there is no intersection (i.e. p is not within 
      // the 'bounds' of the segment) 
      if (t <= (1.0 + EPSILON_TINY)) 
      { 
       // intersects at 2nd segment vertex 
       t = 1.0; 
      } 
      qx = p2x; 
      qy = p2y; 
      // Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if 
      // we were doing PointLineDistanceSquared, then qx would be (p1x + (t * dx)) and qy would be (p1y + (t * dy)). 
     } 
     else 
     { 
      // The projection of the point to the point on the segment that is perpendicular succeeded and the point 
      // is 'within' the bounds of the segment. Set the intersection point as that projected point. 
      qx = ((1.0 - t) * p1x) + (t * p2x); 
      qy = ((1.0 - t) * p1y) + (t * p2y); 
      // for debugging 
      //ASSERT(AreValuesEqual(qx, p1x + (t * dx), EPSILON_TINY)); 
      //ASSERT(AreValuesEqual(qy, p1y + (t * dy), EPSILON_TINY)); 
     } 
     // return the squared distance from p to the intersection point. 
     double dpqx = px - qx; 
     double dpqy = py - qy; 
     return ((dpqx * dpqx) + (dpqy * dpqy)); 
    } 
} 


double SegmentSegmentDistanceSquared( double p1x, double p1y, 
             double p2x, double p2y, 
             double p3x, double p3y, 
             double p4x, double p4y, 
             double& qx, double& qy) 
{ 
    // check to make sure both segments are long enough (i.e. verts are farther apart than minimum allowed vert distance). 
    // If 1 or both segments are shorter than this min length, treat them as a single point. 
    double segLen12Squared = PointToPointDistanceSquared(p1x, p1y, p2x, p2y); 
    double segLen34Squared = PointToPointDistanceSquared(p3x, p3y, p4x, p4y); 
    double t = 0.0; 
    double minDist2 = 1E+38; 
    if (segLen12Squared <= EPSILON_MIN_VERTEX_DISTANCE_SQUARED) 
    { 
     qx = p1x; 
     qy = p1y; 
     if (segLen34Squared <= EPSILON_MIN_VERTEX_DISTANCE_SQUARED) 
     { 
      // point to point 
      minDist2 = PointToPointDistanceSquared(p1x, p1y, p3x, p3y); 
     } 
     else 
     { 
      // point - seg 
      minDist2 = PointSegmentDistanceSquared(p1x, p1y, p3x, p3y, p4x, p4y); 
     } 
     return minDist2; 
    } 
    else if (segLen34Squared <= EPSILON_MIN_VERTEX_DISTANCE_SQUARED) 
    { 
     // seg - point 
     minDist2 = PointSegmentDistanceSquared(p3x, p3y, p1x, p1y, p2x, p2y, t, qx, qy); 
     return minDist2; 
    } 

    // if you have a point class and/or methods to do cross products, you can use those here. 
    // This is what we're actually doing: 
    // Point2D delta43(p4x - p3x, p4y - p3y); // dir of p3 -> p4 
    // Point2D delta12(p1x - p2x, p1y - p2y); // dir of p2 -> p1 
    // double d = delta12.Cross2D(delta43); 
    double d = ((p4y - p3y) * (p1x - p2x)) - ((p1y - p2y) * (p4x - p3x)); 
    bool bParallel = AreValuesEqual(d, 0.0, EPSILON_GENERAL); 

    if (!bParallel) 
    { 
     // segments are not parallel. Check for intersection. 
     // Point2D delta42(p4x - p2x, p4y - p2y); // dir of p2 -> p4 
     // t = 1.0 - (delta42.Cross2D(delta43)/d); 
     t = 1.0 - ((((p4y - p3y) * (p4x - p2x)) - ((p4y - p2y) * (p4x - p3x)))/d); 
     double seg12TEps = sqrt(EPSILON_MIN_VERTEX_DISTANCE_SQUARED/segLen12Squared); 
     if (t >= -seg12TEps && t <= (1.0 + seg12TEps)) 
     { 
      // inside [p1,p2]. Segments may intersect. 
      // double s = 1.0 - (delta12.Cross2D(delta42)/d); 
      double s = 1.0 - ((((p4y - p2y) * (p1x - p2x)) - ((p1y - p2y) * (p4x - p2x)))/d); 
      double seg34TEps = sqrt(EPSILON_MIN_VERTEX_DISTANCE_SQUARED/segLen34Squared); 
      if (s >= -seg34TEps && s <= (1.0 + seg34TEps)) 
      { 
       // segments intersect! 
       minDist2 = 0.0; 
       qx = ((1.0 - t) * p1x) + (t * p2x); 
       qy = ((1.0 - t) * p1y) + (t * p2y); 
       // for debugging 
       //double qsx = ((1.0 - s) * p3x) + (s * p4x); 
       //double qsy = ((1.0 - s) * p3y) + (s * p4y); 
       //ASSERT(AreValuesEqual(qx, qsx, EPSILON_MIN_VERTEX_DISTANCE_SQUARED)); 
       //ASSERT(AreValuesEqual(qy, qsy, EPSILON_MIN_VERTEX_DISTANCE_SQUARED)); 
       return minDist2; 
      } 
     } 
    } 

    // Segments do not intersect. Find closest point and return dist. No other way at this 
    // point except to just brute-force check each segment end-point vs opposite segment. The 
    // minimum distance of those 4 tests is the closest point. 
    double tmpQx, tmpQy, tmpD2; 
    minDist2 = PointSegmentDistanceSquared(p3x, p3y, p1x, p1y, p2x, p2y, t, qx, qy); 
    tmpD2 = PointSegmentDistanceSquared(p4x, p4y, p1x, p1y, p2x, p2y, t, tmpQx, tmpQy); 
    if (tmpD2 < minDist2) 
    { 
     qx = tmpQx; 
     qy = tmpQy; 
     minDist2 = tmpD2; 
    } 
    tmpD2 = PointSegmentDistanceSquared(p1x, p1y, p3x, p3y, p4x, p4y, t, tmpQx, tmpQy); 
    if (tmpD2 < minDist2) 
    { 
     qx = p1x; 
     qy = p1y; 
     minDist2 = tmpD2; 
    } 
    tmpD2 = PointSegmentDistanceSquared(p2x, p2y, p3x, p3y, p4x, p4y, t, tmpQx, tmpQy); 
    if (tmpD2 < minDist2) 
    { 
     qx = p2x; 
     qy = p2y; 
     minDist2 = tmpD2; 
    } 

    return minDist2; 
} 
+0

Questo non compila, la riga 'minDist2 = PointSegmentDistanceSquared (p1x, p1y, p3x, p3y, p4x, p4y);' si aspetta argomenti aggiuntivi. – Richard

+0

Ah, mi dispiace per quello - apparentemente ho dimenticato di includere la versione del metodo che non ha bisogno di t, qx e qy. Dal momento che questi sono solo argomenti di ritorno comunque, puoi semplicemente aggiungere dei valori fittizi lì e ignorare i risultati restituiti. – devnullicus