Najkrótsza odległość między punktem a segmentem linii

358

Potrzebuję podstawowej funkcji, aby znaleźć najkrótszą odległość między punktem a segmentem linii. Napisz rozwiązanie w dowolnym języku; Potrafię przetłumaczyć to na to, czego używam (JavaScript).

EDYCJA: Mój segment linii jest zdefiniowany przez dwa punkty końcowe. Więc mój segment linii ABjest zdefiniowany przez dwa punkty A (x1,y1)i B (x2,y2). Próbuję znaleźć odległość między tym odcinkiem linii a punktem C (x3,y3). Moje umiejętności geometrii są zardzewiałe, więc przykłady, które widziałem, są mylące, przykro mi to przyznać.

Eli Courtwright
źródło
Nie wiem, jak reprezentujesz linie i punkty, ale oto cała matematyka, której potrzebujesz, aby zacząć. Nie powinno być zbyt trudne, aby dowiedzieć się, co musisz zrobić.
mandaleeka
4
@ArthurKalliokoski: ten link jest martwy, ale znalazłem kopię: paulbourke.net/geometry/pointline
Gunther Struyf
11
@GuntherStruyf: ten link też jest martwy, ale ten podobny link działa: paulbourke.net/geometry/pointlineplane
Michael
1
Jeśli ktoś szuka odległości między punktem a linią, a nie punktu i linii SEGMENT, sprawdź ten link: gist.github.com/rhyolight/2846020
Nick Budden
1
Powyższy link jest martwy. Oto pseudo-kod i próbka c ++, wyjaśniona i wyprowadzona tak szczegółowo, jak podręcznik, geomalgorithms.com/a02-_lines.html
Eric

Odpowiedzi:

447

Eli, ustalony kod jest niepoprawny. Punkt w pobliżu linii, na której leży segment, ale daleko od jednego końca segmentu, zostałby nieprawidłowo oceniony w pobliżu segmentu. Aktualizacja: Wspomniana nieprawidłowa odpowiedź nie jest już zaakceptowana.

Oto poprawny kod w C ++. Zakłada klasę wektora 2D class vec2 {float x,y;}, w zasadzie, z operatorami do dodawania, odejmowania, skalowania itp. Oraz funkcją iloczynu odległości i kropki (tj x1 x2 + y1 y2.).

float minimum_distance(vec2 v, vec2 w, vec2 p) {
  // Return minimum distance between line segment vw and point p
  const float l2 = length_squared(v, w);  // i.e. |w-v|^2 -  avoid a sqrt
  if (l2 == 0.0) return distance(p, v);   // v == w case
  // Consider the line extending the segment, parameterized as v + t (w - v).
  // We find projection of point p onto the line. 
  // It falls where t = [(p-v) . (w-v)] / |w-v|^2
  // We clamp t from [0,1] to handle points outside the segment vw.
  const float t = max(0, min(1, dot(p - v, w - v) / l2));
  const vec2 projection = v + t * (w - v);  // Projection falls on the segment
  return distance(p, projection);
}

EDYCJA: Potrzebowałem implementacji Javascript, więc tutaj jest bez zależności (ani komentarzy, ale jest to bezpośredni port powyższego). Punkty są reprezentowane jako obiekty xi yatrybuty.

function sqr(x) { return x * x }
function dist2(v, w) { return sqr(v.x - w.x) + sqr(v.y - w.y) }
function distToSegmentSquared(p, v, w) {
  var l2 = dist2(v, w);
  if (l2 == 0) return dist2(p, v);
  var t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
  t = Math.max(0, Math.min(1, t));
  return dist2(p, { x: v.x + t * (w.x - v.x),
                    y: v.y + t * (w.y - v.y) });
}
function distToSegment(p, v, w) { return Math.sqrt(distToSegmentSquared(p, v, w)); }

EDYCJA 2: Potrzebowałem wersji Java, ale co ważniejsze, potrzebowałem jej w 3d zamiast 2d.

float dist_to_segment_squared(float px, float py, float pz, float lx1, float ly1, float lz1, float lx2, float ly2, float lz2) {
  float line_dist = dist_sq(lx1, ly1, lz1, lx2, ly2, lz2);
  if (line_dist == 0) return dist_sq(px, py, pz, lx1, ly1, lz1);
  float t = ((px - lx1) * (lx2 - lx1) + (py - ly1) * (ly2 - ly1) + (pz - lz1) * (lz2 - lz1)) / line_dist;
  t = constrain(t, 0, 1);
  return dist_sq(px, py, pz, lx1 + t * (lx2 - lx1), ly1 + t * (ly2 - ly1), lz1 + t * (lz2 - lz1));
}
Grumdrig
źródło
1
Jako osobną odpowiedź dodałem rozwiniętą wersję tego.
M Katz,
4
Dzięki @Grumdrig Twoje rozwiązanie w javascript było na miejscu i pozwoliło zaoszczędzić sporo czasu. Przeniesiłem twoje rozwiązanie do Celu-C i dodałem je poniżej.
awolf
1
Naprawdę staramy się uniknąć podziału przez zero.
Grumdrig
9
Rzut punktu pna linię to punkt na linii najbliższej p. (I prostopadłe do linii na występ przechodzi przez p). Liczba t, jak daleko wzdłuż odcinka linii od punktu vna wtym, że występ spada. Więc jeśli twynosi 0, rzutowanie zaczyna się od razu v; jeśli to 1, to jest włączone w; jeśli wynosi na przykład 0,5, to jest w połowie drogi między. Jeśli tjest mniejsze niż 0 lub większe niż 1, spada na linię za jednym końcem lub drugim segmentem. W takim przypadku odległość do segmentu będzie odległością do bliższego końca.
Grumdrig
1
Ups - nie zauważyłem, że ktoś dostarczył wersję 3D. @RogiSolorzano, musisz najpierw przekonwertować współrzędne łacińskie, długie na współrzędne x, y, z w 3 spacji.
Grumdrig
110

Oto najprostszy kompletny kod w JavaScript.

x, y to punkt docelowy, a x1, y1 do x2, y2 to segment linii.

ZAKTUALIZOWANO: naprawiono problem z linią długości 0 z komentarzy.

function pDistance(x, y, x1, y1, x2, y2) {

  var A = x - x1;
  var B = y - y1;
  var C = x2 - x1;
  var D = y2 - y1;

  var dot = A * C + B * D;
  var len_sq = C * C + D * D;
  var param = -1;
  if (len_sq != 0) //in case of 0 length line
      param = dot / len_sq;

  var xx, yy;

  if (param < 0) {
    xx = x1;
    yy = y1;
  }
  else if (param > 1) {
    xx = x2;
    yy = y2;
  }
  else {
    xx = x1 + param * C;
    yy = y1 + param * D;
  }

  var dx = x - xx;
  var dy = y - yy;
  return Math.sqrt(dx * dx + dy * dy);
}

Obraz, który pomaga zwizualizować rozwiązanie

Joshua
źródło
8
Ze wszystkich kodów, które widziałem, aby rozwiązać ten problem, najbardziej podoba mi się ten. Jest bardzo jasny i łatwy do odczytania. Jednak matematyka jest nieco mistyczna. Co na przykład naprawdę reprezentuje iloczyn skalarny podzielony przez kwadrat długości?
user1815201
2
Iloczyn iloczynu podzielony przez kwadrat do długości daje odległość projekcji od (x1, y1). Jest to ułamek linii, do którego punkt (x, y) jest najbliżej. Zwróć uwagę na ostatnią klauzulę else, w której obliczane jest (xx, yy) - jest to rzut punktu (x, y) na segment (x1, y1) - (x2, y2).
Logan Pickup
4
Sprawdzanie segmentów linii o długości 0 jest zbyt daleko w kodzie. „len_sq” będzie wynosić zero, a kod podzieli się przez 0, zanim przejdzie do kontroli bezpieczeństwa.
HostedMetrics.com
17
Metry Zwraca się go w metrach.
Joshua
1
@nevermind, nazwijmy nasz punkt p0 i punkty, które definiują linię jako p1 i p2. Otrzymujemy wtedy wektory A = p0 - p1 i B = p2 - p1. Param to wartość skalarna, która po pomnożeniu przez B daje punkt na linii najbliższej p0. Jeśli param <= 0, najbliższym punktem jest p1. Jeśli param> = 1, najbliższym punktem jest p1. Jeśli jest między 0 a 1, to jest gdzieś pomiędzy p1 i p2, więc interpolujemy. XX i YY to najbliższy punkt na odcinku linii, dx / dy to wektor od p0 do tego punktu, a na koniec zwracamy długość tego wektora.
Sean
70

Jest to implementacja stworzona dla SEGMENTÓW LINII SZYBKICH, a nie nieskończonych linii, jak się wydaje w większości innych funkcji (dlatego to zrobiłem).

Implementacja teorii przez Paula Bourke .

Pyton:

def dist(x1, y1, x2, y2, x3, y3): # x3,y3 is the point
    px = x2-x1
    py = y2-y1

    norm = px*px + py*py

    u =  ((x3 - x1) * px + (y3 - y1) * py) / float(norm)

    if u > 1:
        u = 1
    elif u < 0:
        u = 0

    x = x1 + u * px
    y = y1 + u * py

    dx = x - x3
    dy = y - y3

    # Note: If the actual distance does not matter,
    # if you only want to compare what this function
    # returns to other results of this function, you
    # can just return the squared distance instead
    # (i.e. remove the sqrt) to gain a little performance

    dist = (dx*dx + dy*dy)**.5

    return dist

AS3:

public static function segmentDistToPoint(segA:Point, segB:Point, p:Point):Number
{
    var p2:Point = new Point(segB.x - segA.x, segB.y - segA.y);
    var something:Number = p2.x*p2.x + p2.y*p2.y;
    var u:Number = ((p.x - segA.x) * p2.x + (p.y - segA.y) * p2.y) / something;

    if (u > 1)
        u = 1;
    else if (u < 0)
        u = 0;

    var x:Number = segA.x + u * p2.x;
    var y:Number = segA.y + u * p2.y;

    var dx:Number = x - p.x;
    var dy:Number = y - p.y;

    var dist:Number = Math.sqrt(dx*dx + dy*dy);

    return dist;
}

Jawa

private double shortestDistance(float x1,float y1,float x2,float y2,float x3,float y3)
    {
        float px=x2-x1;
        float py=y2-y1;
        float temp=(px*px)+(py*py);
        float u=((x3 - x1) * px + (y3 - y1) * py) / (temp);
        if(u>1){
            u=1;
        }
        else if(u<0){
            u=0;
        }
        float x = x1 + u * px;
        float y = y1 + u * py;

        float dx = x - x3;
        float dy = y - y3;
        double dist = Math.sqrt(dx*dx + dy*dy);
        return dist;

    }
quano
źródło
2
Przepraszam, ale próbowałem tego i nadal daje mi to wyniki, jakby linia rozciągała się w nieskończoność. Znalazłem jednak odpowiedź Grumdiga na pracę.
Frederik
1
W takim przypadku używasz go źle lub oznacza coś innego z nieskończonym. Zobacz przykład tego kodu tutaj: boomie.se/upload/Drawdebug.swf
quano
Wygląda na błąd w kodzie lub coś takiego, otrzymuję taki sam wynik jak Frederik /
Kromster
30
Wybór nazw zmiennych jest daleki od dobrego (p2, coś, u, ...)
miguelSantirso
2
Wypróbowałem wersję funkcji w języku Python i okazało się, że pokazuje niepoprawne wyniki, jeśli parametry są liczbami całkowitymi. distAnother(0, 0, 4, 0, 2, 2)daje 2.8284271247461903 (niepoprawnie). distAnother(0., 0., 4., 0., 2., 2.)daje 2.0 (poprawnie). Pamiętaj o tym. Myślę, że kod można ulepszyć, aby mieć gdzieś konwersję float.
Vladimir Obrizan
22

W moim wątku pytania, jak obliczyć najkrótszą odległość 2D między punktem a segmentem linii we wszystkich przypadkach w C, C # / .NET 2.0 lub Java? Poproszono mnie o umieszczenie tutaj odpowiedzi w języku C #, więc tutaj jest zmodyfikowana z http://www.topcoder.com/tc?d1=tutorials&d2=geometry1&module=Static :

//Compute the dot product AB . BC
private double DotProduct(double[] pointA, double[] pointB, double[] pointC)
{
    double[] AB = new double[2];
    double[] BC = new double[2];
    AB[0] = pointB[0] - pointA[0];
    AB[1] = pointB[1] - pointA[1];
    BC[0] = pointC[0] - pointB[0];
    BC[1] = pointC[1] - pointB[1];
    double dot = AB[0] * BC[0] + AB[1] * BC[1];

    return dot;
}

//Compute the cross product AB x AC
private double CrossProduct(double[] pointA, double[] pointB, double[] pointC)
{
    double[] AB = new double[2];
    double[] AC = new double[2];
    AB[0] = pointB[0] - pointA[0];
    AB[1] = pointB[1] - pointA[1];
    AC[0] = pointC[0] - pointA[0];
    AC[1] = pointC[1] - pointA[1];
    double cross = AB[0] * AC[1] - AB[1] * AC[0];

    return cross;
}

//Compute the distance from A to B
double Distance(double[] pointA, double[] pointB)
{
    double d1 = pointA[0] - pointB[0];
    double d2 = pointA[1] - pointB[1];

    return Math.Sqrt(d1 * d1 + d2 * d2);
}

//Compute the distance from AB to C
//if isSegment is true, AB is a segment, not a line.
double LineToPointDistance2D(double[] pointA, double[] pointB, double[] pointC, 
    bool isSegment)
{
    double dist = CrossProduct(pointA, pointB, pointC) / Distance(pointA, pointB);
    if (isSegment)
    {
        double dot1 = DotProduct(pointA, pointB, pointC);
        if (dot1 > 0) 
            return Distance(pointB, pointC);

        double dot2 = DotProduct(pointB, pointA, pointC);
        if (dot2 > 0) 
            return Distance(pointA, pointC);
    }
    return Math.Abs(dist);
} 

Nie jestem w stanie odpowiedzieć, ale zadawać pytania, więc mam nadzieję, że z pewnych powodów nie otrzymam miliona głosów, ale skonstruuję krytykę. Chciałem tylko (i zachęciłem) do dzielenia się czyimiś pomysłami, ponieważ rozwiązania w tym wątku są albo z jakimś egzotycznym językiem (Fortran, Mathematica) lub oznaczone przez kogoś jako wadliwe. Jedyny użyteczny (dla mnie Grumdrig) jest napisany w C ++ i nikt nie oznaczył go jako wadliwy. Brakuje jednak wywoływanych metod (kropka itp.).

char m
źródło
1
Dzięki za opublikowanie tego. Ale wygląda na to, że w ostatniej metodzie możliwa jest oczywista optymalizacja: nie obliczaj dist, dopóki nie zostanie ustalone, że jest potrzebne.
RenniePet
2
Komentarz do DotProduct mówi, że oblicza AB.AC, ale oblicza AB.BC.
Metal450
Produkt krzyżowy z definicji zwraca wektor, ale zwraca tutaj skalar.
SteakOverflow
21

W F # odległość od punktu cdo odcinka linii pomiędzy ai bjest określona przez:

let pointToLineSegmentDistance (a: Vector, b: Vector) (c: Vector) =
  let d = b - a
  let s = d.Length
  let lambda = (c - a) * d / s
  let p = (lambda |> max 0.0 |> min s) * d / s
  (a + p - c).Length

Wektor dwskazuje od ado bwzdłuż odcinka linii. Iloczyn iloczynu d/sz c-adaje parametr punktu najbliższego zbliżenia między linią nieskończoną a punktem c. Funkcje mini maxsłużą do zaciskania tego parametru do zakresu, 0..stak aby punkt znajdował się między aa b. Wreszcie długość a+p-cto odległość odc najbliższego punktu na odcinku linii.

Przykładowe zastosowanie:

pointToLineSegmentDistance (Vector(0.0, 0.0), Vector(1.0, 0.0)) (Vector(-1.0, 1.0))
Jon Harrop
źródło
1
Myślę, że ostatni wiersz jest niepoprawny i powinien brzmieć:(a + p - c).Length
Blair Holloway
To wciąż nie rozwiązuje w pełni problemu. Jednym ze sposobów, aby poprawić funkcję będzie przedefiniować lambdai pjak let lambda = (c - a) * d / (s * s)i let p = a + (lambda |> max 0.0 |> min 1.0) * d, odpowiednio. Następnie funkcja zwraca prawidłową odległość, np. Dla przypadku a = (0,1), gdy b = (1,0)i c = (1,1).
mikkoma
20

Dla wszystkich zainteresowanych, oto trywialna konwersja kodu JavaScript Joshua do Objective-C:

- (double)distanceToPoint:(CGPoint)p fromLineSegmentBetween:(CGPoint)l1 and:(CGPoint)l2
{
    double A = p.x - l1.x;
    double B = p.y - l1.y;
    double C = l2.x - l1.x;
    double D = l2.y - l1.y;

    double dot = A * C + B * D;
    double len_sq = C * C + D * D;
    double param = dot / len_sq;

    double xx, yy;

    if (param < 0 || (l1.x == l2.x && l1.y == l2.y)) {
        xx = l1.x;
        yy = l1.y;
    }
    else if (param > 1) {
        xx = l2.x;
        yy = l2.y;
    }
    else {
        xx = l1.x + param * C;
        yy = l1.y + param * D;
    }

    double dx = p.x - xx;
    double dy = p.y - yy;

    return sqrtf(dx * dx + dy * dy);
}

Potrzebowałem tego rozwiązania do pracy, MKMapPointwięc podzielę się nim na wypadek, gdyby ktoś go potrzebował. Tylko niewielka zmiana, która zwróci odległość w metrach:

- (double)distanceToPoint:(MKMapPoint)p fromLineSegmentBetween:(MKMapPoint)l1 and:(MKMapPoint)l2
{
    double A = p.x - l1.x;
    double B = p.y - l1.y;
    double C = l2.x - l1.x;
    double D = l2.y - l1.y;

    double dot = A * C + B * D;
    double len_sq = C * C + D * D;
    double param = dot / len_sq;

    double xx, yy;

    if (param < 0 || (l1.x == l2.x && l1.y == l2.y)) {
        xx = l1.x;
        yy = l1.y;
    }
    else if (param > 1) {
        xx = l2.x;
        yy = l2.y;
    }
    else {
        xx = l1.x + param * C;
        yy = l1.y + param * D;
    }

    return MKMetersBetweenMapPoints(p, MKMapPointMake(xx, yy));
}
Ben Gotow
źródło
To wydaje się działać dobrze dla mnie. Dzięki za konwersję.
Gregir
Warto zauważyć, że (xx, yy) to lokalizacja najbliższego punktu. Zmodyfikowałem trochę twój kod, więc zwraca zarówno punkt, jak i odległość, zrefaktoryzowane nazwy, więc opisują to, co jest i podany przykład: stackoverflow.com/a/28028023/849616 .
Vive
20

W Mathematica

Wykorzystuje parametryczny opis segmentu i rzutuje punkt na linię zdefiniowaną przez segment. Ponieważ parametr zmienia się od 0 do 1 w segmencie, jeśli rzut znajduje się poza tymi granicami, obliczamy odległość do odpowiedniego punktu, zamiast linii prostej prostopadłej do segmentu.

Clear["Global`*"];
 distance[{start_, end_}, pt_] := 
   Module[{param},
   param = ((pt - start).(end - start))/Norm[end - start]^2; (*parameter. the "."
                                                       here means vector product*)

   Which[
    param < 0, EuclideanDistance[start, pt],                 (*If outside bounds*)
    param > 1, EuclideanDistance[end, pt],
    True, EuclideanDistance[pt, start + param (end - start)] (*Normal distance*)
    ]
   ];  

Wynik wydruku:

Plot3D[distance[{{0, 0}, {1, 0}}, {xp, yp}], {xp, -1, 2}, {yp, -1, 2}]

alternatywny tekst

Wykreśl te punkty bliżej niż odległość odcięcia :

alternatywny tekst

Wykres konturowy:

wprowadź opis zdjęcia tutaj

Dr. Belisarius
źródło
11

Hej, właśnie to napisałem wczoraj. Jest w ActionScript 3.0, który jest w zasadzie Javascript, chociaż możesz nie mieć tej samej klasy Point.

//st = start of line segment
//b = the line segment (as in: st + b = end of line segment)
//pt = point to test
//Returns distance from point to line segment.  
//Note: nearest point on the segment to the test point is right there if we ever need it
public static function linePointDist( st:Point, b:Point, pt:Point ):Number
{
    var nearestPt:Point; //closest point on seqment to pt

    var keyDot:Number = dot( b, pt.subtract( st ) ); //key dot product
    var bLenSq:Number = dot( b, b ); //Segment length squared

    if( keyDot <= 0 )  //pt is "behind" st, use st
    {
        nearestPt = st  
    }
    else if( keyDot >= bLenSq ) //pt is "past" end of segment, use end (notice we are saving twin sqrts here cuz)
    {
        nearestPt = st.add(b);
    }
    else //pt is inside segment, reuse keyDot and bLenSq to get percent of seqment to move in to find closest point
    {
        var keyDotToPctOfB:Number = keyDot/bLenSq; //REM dot product comes squared
        var partOfB:Point = new Point( b.x * keyDotToPctOfB, b.y * keyDotToPctOfB );
        nearestPt = st.add(partOfB);
    }

    var dist:Number = (pt.subtract(nearestPt)).length;

    return dist;
}

Jest też dość kompletna i czytelna dyskusja na temat problemu: notejot.com

Matt W.
źródło
Dzięki - to jest dokładnie ten rodzaj kodu, którego szukałem. Poniżej zamieściłem własną odpowiedź, ponieważ udało mi się stworzyć coś, co działa w obecnej przeglądarce JavaScript, ale oznaczyłem twoją odpowiedź jako zaakceptowaną, ponieważ jest to prosta, dobrze napisana, łatwa do zrozumienia, i bardzo mile widziane.
Eli Courtwright
Czy nie brakuje metody kropkowej? W każdym razie łatwo jest obliczyć: vec1.x * vec2.x + vec1.y * vec2.y
quano
11

Dla leniwych, oto mój port Objective-C rozwiązania @ Grumdrig powyżej:

CGFloat sqr(CGFloat x) { return x*x; }
CGFloat dist2(CGPoint v, CGPoint w) { return sqr(v.x - w.x) + sqr(v.y - w.y); }
CGFloat distanceToSegmentSquared(CGPoint p, CGPoint v, CGPoint w)
{
    CGFloat l2 = dist2(v, w);
    if (l2 == 0.0f) return dist2(p, v);

    CGFloat t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
    if (t < 0.0f) return dist2(p, v);
    if (t > 1.0f) return dist2(p, w);
    return dist2(p, CGPointMake(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y)));
}
CGFloat distanceToSegment(CGPoint point, CGPoint segmentPointV, CGPoint segmentPointW)
{
    return sqrtf(distanceToSegmentSquared(point, segmentPointV, segmentPointW));
}
awolf
źródło
Dostaję „nan” z tej linii. Masz pomysł, dlaczego? ( return dist2(p, CGPointMake(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y)))
Nawiasem mówiąc,
sqrtf () kwadrat x, nie uzyskuje pierwiastka kwadratowego
Senseful
@Senseful Nie jestem pewien, co masz na myśli. sqrtf jest pierwiastkiem kwadratowym. developer.apple.com/library/mac/documentation/Darwin/Reference/…
awolf
@awolf: Spójrz na pierwszą linię kodu powyżej. Określa metodę sqrtf(x) = x*x.
Rozsądny
@ Dzięki wdzięczny, został źle nazwany, a nie wykonał niewłaściwej operacji.
awolf
10

Nie mogłem się oprzeć zakodowaniu go w pythonie :)

from math import sqrt, fabs
def pdis(a, b, c):
    t = b[0]-a[0], b[1]-a[1]           # Vector ab
    dd = sqrt(t[0]**2+t[1]**2)         # Length of ab
    t = t[0]/dd, t[1]/dd               # unit vector of ab
    n = -t[1], t[0]                    # normal unit vector to ab
    ac = c[0]-a[0], c[1]-a[1]          # vector ac
    return fabs(ac[0]*n[0]+ac[1]*n[1]) # Projection of ac to n (the minimum distance)

print pdis((1,1), (2,2), (2,0))        # Example (answer is 1.414)


Ditto for fortran :)

real function pdis(a, b, c)
    real, dimension(0:1), intent(in) :: a, b, c
    real, dimension(0:1) :: t, n, ac
    real :: dd
    t = b - a                          ! Vector ab
    dd = sqrt(t(0)**2+t(1)**2)         ! Length of ab
    t = t/dd                           ! unit vector of ab
    n = (/-t(1), t(0)/)                ! normal unit vector to ab
    ac = c - a                         ! vector ac
    pdis = abs(ac(0)*n(0)+ac(1)*n(1))  ! Projection of ac to n (the minimum distance)
end function pdis


program test
    print *, pdis((/1.0,1.0/), (/2.0,2.0/), (/2.0,0.0/))   ! Example (answer is 1.414)
end program test
cyberthanasis
źródło
10
czy to nie oblicza odległości punktu od linii zamiast segmentu?
balint.miklos
6
Jest to rzeczywiście odległość do linii, na której znajduje się segment, a nie do segmentu.
Grumdrig
12
To nie działa. Jeśli masz segment (0,0) i (5,0) i spróbuj przeciw punktowi (7,0), zwróci 0, co nie jest prawdą. Odległość powinna wynosić 2.
quano
8
Nie wziął pod uwagę przypadku, w którym rzut punktu na segment jest poza przedziałem od A do B. To może być to, czego chciał pytający, ale nie to, o co prosił.
phkahler
5
Nie o to pierwotnie pytano.
Sambatyon
10

Oto pełniejsza pisownia rozwiązania Grumdrig. Ta wersja zwraca również najbliższy punkt.

#include "stdio.h"
#include "math.h"

class Vec2
{
public:
    float _x;
    float _y;

    Vec2()
    {
        _x = 0;
        _y = 0;
    }

    Vec2( const float x, const float y )
    {
        _x = x;
        _y = y;
    }

    Vec2 operator+( const Vec2 &v ) const
    {
        return Vec2( this->_x + v._x, this->_y + v._y );
    }

    Vec2 operator-( const Vec2 &v ) const
    {
        return Vec2( this->_x - v._x, this->_y - v._y );
    }

    Vec2 operator*( const float f ) const
    {
        return Vec2( this->_x * f, this->_y * f );
    }

    float DistanceToSquared( const Vec2 p ) const
    {
        const float dX = p._x - this->_x;
        const float dY = p._y - this->_y;

        return dX * dX + dY * dY;
    }

    float DistanceTo( const Vec2 p ) const
    {
        return sqrt( this->DistanceToSquared( p ) );
    }

    float DotProduct( const Vec2 p ) const
    {
        return this->_x * p._x + this->_y * p._y;
    }
};

// return minimum distance between line segment vw and point p, and the closest point on the line segment, q
float DistanceFromLineSegmentToPoint( const Vec2 v, const Vec2 w, const Vec2 p, Vec2 * const q )
{
    const float distSq = v.DistanceToSquared( w ); // i.e. |w-v|^2 ... avoid a sqrt
    if ( distSq == 0.0 )
    {
        // v == w case
        (*q) = v;

        return v.DistanceTo( p );
    }

    // consider the line extending the segment, parameterized as v + t (w - v)
    // we find projection of point p onto the line
    // it falls where t = [(p-v) . (w-v)] / |w-v|^2

    const float t = ( p - v ).DotProduct( w - v ) / distSq;
    if ( t < 0.0 )
    {
        // beyond the v end of the segment
        (*q) = v;

        return v.DistanceTo( p );
    }
    else if ( t > 1.0 )
    {
        // beyond the w end of the segment
        (*q) = w;

        return w.DistanceTo( p );
    }

    // projection falls on the segment
    const Vec2 projection = v + ( ( w - v ) * t );

    (*q) = projection;

    return p.DistanceTo( projection );
}

float DistanceFromLineSegmentToPoint( float segmentX1, float segmentY1, float segmentX2, float segmentY2, float pX, float pY, float *qX, float *qY )
{
    Vec2 q;

    float distance = DistanceFromLineSegmentToPoint( Vec2( segmentX1, segmentY1 ), Vec2( segmentX2, segmentY2 ), Vec2( pX, pY ), &q );

    (*qX) = q._x;
    (*qY) = q._y;

    return distance;
}

void TestDistanceFromLineSegmentToPoint( float segmentX1, float segmentY1, float segmentX2, float segmentY2, float pX, float pY )
{
    float qX;
    float qY;
    float d = DistanceFromLineSegmentToPoint( segmentX1, segmentY1, segmentX2, segmentY2, pX, pY, &qX, &qY );
    printf( "line segment = ( ( %f, %f ), ( %f, %f ) ), p = ( %f, %f ), distance = %f, q = ( %f, %f )\n",
            segmentX1, segmentY1, segmentX2, segmentY2, pX, pY, d, qX, qY );
}

void TestDistanceFromLineSegmentToPoint()
{
    TestDistanceFromLineSegmentToPoint( 0, 0, 1, 1, 1, 0 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 20, 10, 5, 4 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 20, 10, 30, 15 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 20, 10, -30, 15 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 10, 0, 5, 1 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 0, 10, 1, 5 );
}
M. Katz
źródło
Dzięki za opublikowanie tego. Bardzo dobrze skonstruowany, skomentowany i sformatowany - prawie sprawił, że zapomniałem, jak bardzo nie lubię C ++. Użyłem tego, aby utworzyć odpowiednią wersję C #, którą teraz opublikowałem tutaj.
RenniePet
10

Jedno liniowe rozwiązanie wykorzystujące arcus tangens:

Chodzi o to, aby przesunąć A do (0, 0) i obrócić trójkąt zgodnie z ruchem wskazówek zegara, aby C leżał na osi X, kiedy to się stanie, By będzie odległością.

  1. kąt = Atan (Cy - Ay, Cx - Axe);
  2. kąt b = Atan (By - Ay, Bx - Axe);
  3. Długość AB = Sqrt ((Bx - Axe) ^ 2 + (By - Ay) ^ 2)
  4. By = Sin (bAngle - aAngle) * ABLength

DO#

public double Distance(Point a, Point b, Point c)
{
    // normalize points
    Point cn = new Point(c.X - a.X, c.Y - a.Y);
    Point bn = new Point(b.X - a.X, b.Y - a.Y);

    double angle = Math.Atan2(bn.Y, bn.X) - Math.Atan2(cn.Y, cn.X);
    double abLength = Math.Sqrt(bn.X*bn.X + bn.Y*bn.Y);

    return Math.Sin(angle)*abLength;
}

Jedna linia C # (do konwersji na SQL)

double distance = Math.Sin(Math.Atan2(b.Y - a.Y, b.X - a.X) - Math.Atan2(c.Y - a.Y, c.X - a.X)) * Math.Sqrt((b.X - a.X) * (b.X - a.X) + (b.Y - a.Y) * (b.Y - a.Y))
ADOConnection
źródło
7

Zastanów się nad modyfikacją powyższej odpowiedzi Grumdriga. Wiele razy przekonasz się, że niedokładność zmiennoprzecinkowa może powodować problemy. Używam podwójnych w poniższej wersji, ale możesz łatwo zmienić na zmiennoprzecinkowe. Ważną częścią jest to, że używa epsilon do obsługi „wpadki”. Ponadto wiele razy będziesz chciał wiedzieć, GDZIE nastąpiło skrzyżowanie, lub jeśli w ogóle. Jeśli zwrócone t wynosi <0,0 lub> 1,0, nie doszło do kolizji. Jednak nawet jeśli nie doszło do kolizji, wiele razy będziesz chciał wiedzieć, gdzie jest najbliższy punkt w segmencie do P, dlatego używam qx i qy, aby zwrócić tę lokalizację.

double PointSegmentDistanceSquared( double px, double py,
                                    double p1x, double p1y,
                                    double p2x, double p2y,
                                    double& t,
                                    double& qx, double& qy)
{
    static const double kMinSegmentLenSquared = 0.00000001;  // adjust to suit.  If you use float, you'll probably want something like 0.000001f
    static const double kEpsilon = 1.0E-14;  // adjust to suit.  If you use floats, you'll probably want something like 1E-7f
    double dx = p2x - p1x;
    double dy = p2y - p1y;
    double dp1x = px - p1x;
    double dp1y = py - p1y;
    const double segLenSquared = (dx * dx) + (dy * dy);
    if (segLenSquared >= -kMinSegmentLenSquared && segLenSquared <= kMinSegmentLenSquared)
    {
        // segment is a point.
        qx = p1x;
        qy = p1y;
        t = 0.0;
        return ((dp1x * dp1x) + (dp1y * dp1y));
    }
    else
    {
        // Project a line from p to the segment [p1,p2].  By considering the line
        // extending the segment, parameterized as p1 + (t * (p2 - p1)),
        // we find projection of point p onto the line. 
        // It falls where t = [(p - p1) . (p2 - p1)] / |p2 - p1|^2
        t = ((dp1x * dx) + (dp1y * dy)) / segLenSquared;
        if (t < kEpsilon)
        {
            // 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 > -kEpsilon)
            {
                // 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 - kEpsilon))
        {
            // 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 + kEpsilon))
            {
                // intersects at 2nd segment vertex
                t = 1.0;
            }
            // set our 'intersection' point to p2.
            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 = p1x + (t * dx);
            qy = p1y + (t * dy);
        }
        // return the squared distance from p to the intersection point.  Note that we return the squared distance
        // as an optimization because many times you just need to compare relative distances and the squared values
        // works fine for that.  If you want the ACTUAL distance, just take the square root of this value.
        double dpqx = px - qx;
        double dpqy = py - qy;
        return ((dpqx * dpqx) + (dpqy * dpqy));
    }
}
devnullicus
źródło
6

Zakładam, że chcesz znaleźć najkrótszyodległość między punktem a odcinkiem linii; aby to zrobić, musisz znaleźć linię (linię A), która jest prostopadła do segmentu linii (linia B), która przechodzi przez twój punkt, określić przecięcie między tą linią (linia A) a linią, która przechodzi przez twój segment linii (linia B) ; jeśli ten punkt znajduje się między dwoma punktami segmentu linii, to odległość jest odległością między twoim punktem a właśnie znalezionym punktem, który jest przecięciem linii A i linii B; jeśli punkt nie znajduje się między dwoma punktami segmentu linii, musisz uzyskać odległość między swoim punktem a bliższą dwoma końcami segmentu linii; można to łatwo zrobić, biorąc odległość kwadratową (aby uniknąć pierwiastka kwadratowego) między punktem a dwoma punktami segmentu linii; cokolwiek jest bliżej, weź pierwiastek kwadratowy z tego.

Paul Sonier
źródło
6

Implementacja C ++ / JavaScript Grumdriga była dla mnie bardzo przydatna, więc podałem bezpośredni port Python, którego używam. Pełny kod jest tutaj .

class Point(object):
  def __init__(self, x, y):
    self.x = float(x)
    self.y = float(y)

def square(x):
  return x * x

def distance_squared(v, w):
  return square(v.x - w.x) + square(v.y - w.y)

def distance_point_segment_squared(p, v, w):
  # Segment length squared, |w-v|^2
  d2 = distance_squared(v, w) 
  if d2 == 0: 
    # v == w, return distance to v
    return distance_squared(p, v)
  # Consider the line extending the segment, parameterized as v + t (w - v).
  # We find projection of point p onto the line.
  # It falls where t = [(p-v) . (w-v)] / |w-v|^2
  t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / d2;
  if t < 0:
    # Beyond v end of the segment
    return distance_squared(p, v)
  elif t > 1.0:
    # Beyond w end of the segment
    return distance_squared(p, w)
  else:
    # Projection falls on the segment.
    proj = Point(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y))
    # print proj.x, proj.y
    return distance_squared(p, proj)
wełna
źródło
5

Kod Matlab z wbudowanym „autotestem”, jeśli wywołują funkcję bez argumentów:

function r = distPointToLineSegment( xy0, xy1, xyP )
% r = distPointToLineSegment( xy0, xy1, xyP )

if( nargin < 3 )
    selfTest();
    r=0;
else
    vx = xy0(1)-xyP(1);
    vy = xy0(2)-xyP(2);
    ux = xy1(1)-xy0(1);
    uy = xy1(2)-xy0(2);
    lenSqr= (ux*ux+uy*uy);
    detP= -vx*ux + -vy*uy;

    if( detP < 0 )
        r = norm(xy0-xyP,2);
    elseif( detP > lenSqr )
        r = norm(xy1-xyP,2);
    else
        r = abs(ux*vy-uy*vx)/sqrt(lenSqr);
    end
end


    function selfTest()
        %#ok<*NASGU>
        disp(['invalid args, distPointToLineSegment running (recursive)  self-test...']);

        ptA = [1;1]; ptB = [-1;-1];
        ptC = [1/2;1/2];  % on the line
        ptD = [-2;-1.5];  % too far from line segment
        ptE = [1/2;0];    % should be same as perpendicular distance to line
        ptF = [1.5;1.5];      % along the A-B but outside of the segment

        distCtoAB = distPointToLineSegment(ptA,ptB,ptC)
        distDtoAB = distPointToLineSegment(ptA,ptB,ptD)
        distEtoAB = distPointToLineSegment(ptA,ptB,ptE)
        distFtoAB = distPointToLineSegment(ptA,ptB,ptF)
        figure(1); clf;
        circle = @(x, y, r, c) rectangle('Position', [x-r, y-r, 2*r, 2*r], ...
            'Curvature', [1 1], 'EdgeColor', c);
        plot([ptA(1) ptB(1)],[ptA(2) ptB(2)],'r-x'); hold on;
        plot(ptC(1),ptC(2),'b+'); circle(ptC(1),ptC(2), 0.5e-1, 'b');
        plot(ptD(1),ptD(2),'g+'); circle(ptD(1),ptD(2), distDtoAB, 'g');
        plot(ptE(1),ptE(2),'k+'); circle(ptE(1),ptE(2), distEtoAB, 'k');
        plot(ptF(1),ptF(2),'m+'); circle(ptF(1),ptF(2), distFtoAB, 'm');
        hold off;
        axis([-3 3 -3 3]); axis equal;
    end

end
Peter Karasev
źródło
Dzięki, ten kod Matlaba faktycznie oblicza najkrótszą odległość do linii SEGMENT, a nie odległość do nieskończonej linii, na której leży segment.
Rudolf Meijering
4

A teraz moje rozwiązanie również ...... (JavaScript)

Jest bardzo szybki, ponieważ staram się unikać jakichkolwiek funkcji Math.pow.

Jak widać na końcu funkcji mam odległość linii.

kod pochodzi z lib http://www.draw2d.org/graphiti/jsdoc/#!/example

/**
 * Static util function to determine is a point(px,py) on the line(x1,y1,x2,y2)
 * A simple hit test.
 * 
 * @return {boolean}
 * @static
 * @private
 * @param {Number} coronaWidth the accepted corona for the hit test
 * @param {Number} X1 x coordinate of the start point of the line
 * @param {Number} Y1 y coordinate of the start point of the line
 * @param {Number} X2 x coordinate of the end point of the line
 * @param {Number} Y2 y coordinate of the end point of the line
 * @param {Number} px x coordinate of the point to test
 * @param {Number} py y coordinate of the point to test
 **/
graphiti.shape.basic.Line.hit= function( coronaWidth, X1, Y1,  X2,  Y2, px, py)
{
  // Adjust vectors relative to X1,Y1
  // X2,Y2 becomes relative vector from X1,Y1 to end of segment
  X2 -= X1;
  Y2 -= Y1;
  // px,py becomes relative vector from X1,Y1 to test point
  px -= X1;
  py -= Y1;
  var dotprod = px * X2 + py * Y2;
  var projlenSq;
  if (dotprod <= 0.0) {
      // px,py is on the side of X1,Y1 away from X2,Y2
      // distance to segment is length of px,py vector
      // "length of its (clipped) projection" is now 0.0
      projlenSq = 0.0;
  } else {
      // switch to backwards vectors relative to X2,Y2
      // X2,Y2 are already the negative of X1,Y1=>X2,Y2
      // to get px,py to be the negative of px,py=>X2,Y2
      // the dot product of two negated vectors is the same
      // as the dot product of the two normal vectors
      px = X2 - px;
      py = Y2 - py;
      dotprod = px * X2 + py * Y2;
      if (dotprod <= 0.0) {
          // px,py is on the side of X2,Y2 away from X1,Y1
          // distance to segment is length of (backwards) px,py vector
          // "length of its (clipped) projection" is now 0.0
          projlenSq = 0.0;
      } else {
          // px,py is between X1,Y1 and X2,Y2
          // dotprod is the length of the px,py vector
          // projected on the X2,Y2=>X1,Y1 vector times the
          // length of the X2,Y2=>X1,Y1 vector
          projlenSq = dotprod * dotprod / (X2 * X2 + Y2 * Y2);
      }
  }
    // Distance to line is now the length of the relative point
    // vector minus the length of its projection onto the line
    // (which is zero if the projection falls outside the range
    //  of the line segment).
    var lenSq = px * px + py * py - projlenSq;
    if (lenSq < 0) {
        lenSq = 0;
    }
    return Math.sqrt(lenSq)<coronaWidth;
};
użytkownik1397870
źródło
4

zakodowane w t-sql

punkt to (@px, @py), a odcinek linii biegnie od (@ax, @ay) do (@bx, @by)

create function fn_sqr (@NumberToSquare decimal(18,10)) 
returns decimal(18,10)
as 
begin
    declare @Result decimal(18,10)
    set @Result = @NumberToSquare * @NumberToSquare
    return @Result
end
go

create function fn_Distance(@ax decimal (18,10) , @ay decimal (18,10), @bx decimal(18,10),  @by decimal(18,10)) 
returns decimal(18,10)
as
begin
    declare @Result decimal(18,10)
    set @Result = (select dbo.fn_sqr(@ax - @bx) + dbo.fn_sqr(@ay - @by) )
    return @Result
end
go

create function fn_DistanceToSegmentSquared(@px decimal(18,10), @py decimal(18,10), @ax decimal(18,10), @ay decimal(18,10), @bx decimal(18,10), @by decimal(18,10)) 
returns decimal(18,10)
as 
begin
    declare @l2 decimal(18,10)
    set @l2 = (select dbo.fn_Distance(@ax, @ay, @bx, @by))
    if @l2 = 0
        return dbo.fn_Distance(@px, @py, @ax, @ay)
    declare @t decimal(18,10)
    set @t = ((@px - @ax) * (@bx - @ax) + (@py - @ay) * (@by - @ay)) / @l2
    if (@t < 0) 
        return dbo.fn_Distance(@px, @py, @ax, @ay);
    if (@t > 1) 
        return dbo.fn_Distance(@px, @py, @bx, @by);
    return dbo.fn_Distance(@px, @py,  @ax + @t * (@bx - @ax),  @ay + @t * (@by - @ay))
end
go

create function fn_DistanceToSegment(@px decimal(18,10), @py decimal(18,10), @ax decimal(18,10), @ay decimal(18,10), @bx decimal(18,10), @by decimal(18,10)) 
returns decimal(18,10)
as 
begin
    return sqrt(dbo.fn_DistanceToSegmentSquared(@px, @py , @ax , @ay , @bx , @by ))
end
go

--example execution for distance from a point at (6,1) to line segment that runs from (4,2) to (2,1)
select dbo.fn_DistanceToSegment(6, 1, 4, 2, 2, 1) 
--result = 2.2360679775

--example execution for distance from a point at (-3,-2) to line segment that runs from (0,-2) to (-2,1)
select dbo.fn_DistanceToSegment(-3, -2, 0, -2, -2, 1) 
--result = 2.4961508830

--example execution for distance from a point at (0,-2) to line segment that runs from (0,-2) to (-2,1)
select dbo.fn_DistanceToSegment(0,-2, 0, -2, -2, 1) 
--result = 0.0000000000
rob mcnicol
źródło
4

Wygląda na to, że prawie wszyscy inni na StackOverflow udzielił odpowiedzi (do tej pory 23 odpowiedzi), więc oto mój wkład w C #. Jest to głównie oparte na odpowiedzi M. Katza, która z kolei oparta jest na odpowiedzi Grumdriga.

   public struct MyVector
   {
      private readonly double _x, _y;


      // Constructor
      public MyVector(double x, double y)
      {
         _x = x;
         _y = y;
      }


      // Distance from this point to another point, squared
      private double DistanceSquared(MyVector otherPoint)
      {
         double dx = otherPoint._x - this._x;
         double dy = otherPoint._y - this._y;
         return dx * dx + dy * dy;
      }


      // Find the distance from this point to a line segment (which is not the same as from this 
      //  point to anywhere on an infinite line). Also returns the closest point.
      public double DistanceToLineSegment(MyVector lineSegmentPoint1, MyVector lineSegmentPoint2,
                                          out MyVector closestPoint)
      {
         return Math.Sqrt(DistanceToLineSegmentSquared(lineSegmentPoint1, lineSegmentPoint2, 
                          out closestPoint));
      }


      // Same as above, but avoid using Sqrt(), saves a new nanoseconds in cases where you only want 
      //  to compare several distances to find the smallest or largest, but don't need the distance
      public double DistanceToLineSegmentSquared(MyVector lineSegmentPoint1, 
                                              MyVector lineSegmentPoint2, out MyVector closestPoint)
      {
         // Compute length of line segment (squared) and handle special case of coincident points
         double segmentLengthSquared = lineSegmentPoint1.DistanceSquared(lineSegmentPoint2);
         if (segmentLengthSquared < 1E-7f)  // Arbitrary "close enough for government work" value
         {
            closestPoint = lineSegmentPoint1;
            return this.DistanceSquared(closestPoint);
         }

         // Use the magic formula to compute the "projection" of this point on the infinite line
         MyVector lineSegment = lineSegmentPoint2 - lineSegmentPoint1;
         double t = (this - lineSegmentPoint1).DotProduct(lineSegment) / segmentLengthSquared;

         // Handle the two cases where the projection is not on the line segment, and the case where 
         //  the projection is on the segment
         if (t <= 0)
            closestPoint = lineSegmentPoint1;
         else if (t >= 1)
            closestPoint = lineSegmentPoint2;
         else 
            closestPoint = lineSegmentPoint1 + (lineSegment * t);
         return this.DistanceSquared(closestPoint);
      }


      public double DotProduct(MyVector otherVector)
      {
         return this._x * otherVector._x + this._y * otherVector._y;
      }

      public static MyVector operator +(MyVector leftVector, MyVector rightVector)
      {
         return new MyVector(leftVector._x + rightVector._x, leftVector._y + rightVector._y);
      }

      public static MyVector operator -(MyVector leftVector, MyVector rightVector)
      {
         return new MyVector(leftVector._x - rightVector._x, leftVector._y - rightVector._y);
      }

      public static MyVector operator *(MyVector aVector, double aScalar)
      {
         return new MyVector(aVector._x * aScalar, aVector._y * aScalar);
      }

      // Added using ReSharper due to CodeAnalysis nagging

      public bool Equals(MyVector other)
      {
         return _x.Equals(other._x) && _y.Equals(other._y);
      }

      public override bool Equals(object obj)
      {
         if (ReferenceEquals(null, obj)) return false;
         return obj is MyVector && Equals((MyVector) obj);
      }

      public override int GetHashCode()
      {
         unchecked
         {
            return (_x.GetHashCode()*397) ^ _y.GetHashCode();
         }
      }

      public static bool operator ==(MyVector left, MyVector right)
      {
         return left.Equals(right);
      }

      public static bool operator !=(MyVector left, MyVector right)
      {
         return !left.Equals(right);
      }
   }

A oto mały program testowy.

   public static class JustTesting
   {
      public static void Main()
      {
         Stopwatch stopwatch = new Stopwatch();
         stopwatch.Start();

         for (int i = 0; i < 10000000; i++)
         {
            TestIt(1, 0, 0, 0, 1, 1, 0.70710678118654757);
            TestIt(5, 4, 0, 0, 20, 10, 1.3416407864998738);
            TestIt(30, 15, 0, 0, 20, 10, 11.180339887498949);
            TestIt(-30, 15, 0, 0, 20, 10, 33.541019662496844);
            TestIt(5, 1, 0, 0, 10, 0, 1.0);
            TestIt(1, 5, 0, 0, 0, 10, 1.0);
         }

         stopwatch.Stop();
         TimeSpan timeSpan = stopwatch.Elapsed;
      }


      private static void TestIt(float aPointX, float aPointY, 
                                 float lineSegmentPoint1X, float lineSegmentPoint1Y, 
                                 float lineSegmentPoint2X, float lineSegmentPoint2Y, 
                                 double expectedAnswer)
      {
         // Katz
         double d1 = DistanceFromPointToLineSegment(new MyVector(aPointX, aPointY), 
                                              new MyVector(lineSegmentPoint1X, lineSegmentPoint1Y), 
                                              new MyVector(lineSegmentPoint2X, lineSegmentPoint2Y));
         Debug.Assert(d1 == expectedAnswer);

         /*
         // Katz using squared distance
         double d2 = DistanceFromPointToLineSegmentSquared(new MyVector(aPointX, aPointY), 
                                              new MyVector(lineSegmentPoint1X, lineSegmentPoint1Y), 
                                              new MyVector(lineSegmentPoint2X, lineSegmentPoint2Y));
         Debug.Assert(Math.Abs(d2 - expectedAnswer * expectedAnswer) < 1E-7f);
          */

         /*
         // Matti (optimized)
         double d3 = FloatVector.DistanceToLineSegment(new PointF(aPointX, aPointY), 
                                                new PointF(lineSegmentPoint1X, lineSegmentPoint1Y), 
                                                new PointF(lineSegmentPoint2X, lineSegmentPoint2Y));
         Debug.Assert(Math.Abs(d3 - expectedAnswer) < 1E-7f);
          */
      }

      private static double DistanceFromPointToLineSegment(MyVector aPoint, 
                                             MyVector lineSegmentPoint1, MyVector lineSegmentPoint2)
      {
         MyVector closestPoint;  // Not used
         return aPoint.DistanceToLineSegment(lineSegmentPoint1, lineSegmentPoint2, 
                                             out closestPoint);
      }

      private static double DistanceFromPointToLineSegmentSquared(MyVector aPoint, 
                                             MyVector lineSegmentPoint1, MyVector lineSegmentPoint2)
      {
         MyVector closestPoint;  // Not used
         return aPoint.DistanceToLineSegmentSquared(lineSegmentPoint1, lineSegmentPoint2, 
                                                    out closestPoint);
      }
   }

Jak widać, próbowałem zmierzyć różnicę między użyciem wersji, która omija metodę Sqrt (), a wersją normalną. Moje testy wskazują, że możesz zaoszczędzić około 2,5%, ale nawet nie jestem tego pewien - różnice w różnych przebiegach testowych były tego samego rzędu wielkości. Próbowałem także zmierzyć wersję opublikowaną przez Matti (plus oczywistą optymalizację), i ta wersja wydaje się być o około 4% wolniejsza niż wersja oparta na kodzie Katz / Grumdrig.

Edycja: Nawiasem mówiąc, próbowałem również zmierzyć metodę, która znajduje odległość do nieskończonej linii (nie segmentu linii) przy użyciu iloczynu krzyżowego (i Sqrt ()), i jest o około 32% szybsza.

RenniePet
źródło
3

Oto wersja C ++ devnullicus przekonwertowana na C #. Do mojej realizacji musiałem znać punkt przecięcia i znaleźć jego rozwiązanie, które działa dobrze.

public static bool PointSegmentDistanceSquared(PointF point, PointF lineStart, PointF lineEnd, out double distance, out PointF intersectPoint)
{
    const double kMinSegmentLenSquared = 0.00000001; // adjust to suit.  If you use float, you'll probably want something like 0.000001f
    const double kEpsilon = 1.0E-14; // adjust to suit.  If you use floats, you'll probably want something like 1E-7f
    double dX = lineEnd.X - lineStart.X;
    double dY = lineEnd.Y - lineStart.Y;
    double dp1X = point.X - lineStart.X;
    double dp1Y = point.Y - lineStart.Y;
    double segLenSquared = (dX * dX) + (dY * dY);
    double t = 0.0;

    if (segLenSquared >= -kMinSegmentLenSquared && segLenSquared <= kMinSegmentLenSquared)
    {
        // segment is a point.
        intersectPoint = lineStart;
        t = 0.0;
        distance = ((dp1X * dp1X) + (dp1Y * dp1Y));
    }
    else
    {
        // Project a line from p to the segment [p1,p2].  By considering the line
        // extending the segment, parameterized as p1 + (t * (p2 - p1)),
        // we find projection of point p onto the line. 
        // It falls where t = [(p - p1) . (p2 - p1)] / |p2 - p1|^2
        t = ((dp1X * dX) + (dp1Y * dY)) / segLenSquared;
        if (t < kEpsilon)
        {
            // intersects at or to the "left" of first segment vertex (lineStart.X, lineStart.Y).  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 > -kEpsilon)
            {
                // intersects at 1st segment vertex
                t = 0.0;
            }
            // set our 'intersection' point to p1.
            intersectPoint = lineStart;
            // Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
            // we were doing PointLineDistanceSquared, then intersectPoint.X would be (lineStart.X + (t * dx)) and intersectPoint.Y would be (lineStart.Y + (t * dy)).
        }
        else if (t > (1.0 - kEpsilon))
        {
            // intersects at or to the "right" of second segment vertex (lineEnd.X, lineEnd.Y).  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 + kEpsilon))
            {
                // intersects at 2nd segment vertex
                t = 1.0;
            }
            // set our 'intersection' point to p2.
            intersectPoint = lineEnd;
            // Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
            // we were doing PointLineDistanceSquared, then intersectPoint.X would be (lineStart.X + (t * dx)) and intersectPoint.Y would be (lineStart.Y + (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.
            intersectPoint = new PointF((float)(lineStart.X + (t * dX)), (float)(lineStart.Y + (t * dY)));
        }
        // return the squared distance from p to the intersection point.  Note that we return the squared distance
        // as an optimization because many times you just need to compare relative distances and the squared values
        // works fine for that.  If you want the ACTUAL distance, just take the square root of this value.
        double dpqX = point.X - intersectPoint.X;
        double dpqY = point.Y - intersectPoint.Y;

        distance = ((dpqX * dpqX) + (dpqY * dpqY));
    }

    return true;
}
headkaze
źródło
Działa jak marzenie!! Uratował mi niezliczone godziny. Dzięki wielkie!!
Steve Johnson
3

Tutaj używa Swift

    /* Distance from a point (p1) to line l1 l2 */
func distanceFromPoint(p: CGPoint, toLineSegment l1: CGPoint, and l2: CGPoint) -> CGFloat {
    let A = p.x - l1.x
    let B = p.y - l1.y
    let C = l2.x - l1.x
    let D = l2.y - l1.y

    let dot = A * C + B * D
    let len_sq = C * C + D * D
    let param = dot / len_sq

    var xx, yy: CGFloat

    if param < 0 || (l1.x == l2.x && l1.y == l2.y) {
        xx = l1.x
        yy = l1.y
    } else if param > 1 {
        xx = l2.x
        yy = l2.y
    } else {
        xx = l1.x + param * C
        yy = l1.y + param * D
    }

    let dx = p.x - xx
    let dy = p.y - yy

    return sqrt(dx * dx + dy * dy)
}
OzRunways
źródło
3

DO#

Na podstawie @Grumdrig

public static double MinimumDistanceToLineSegment(this Point p,
    Line line)
{
    var v = line.StartPoint;
    var w = line.EndPoint;

    double lengthSquared = DistanceSquared(v, w);

    if (lengthSquared == 0.0)
        return Distance(p, v);

    double t = Math.Max(0, Math.Min(1, DotProduct(p - v, w - v) / lengthSquared));
    var projection = v + t * (w - v);

    return Distance(p, projection);
}

public static double Distance(Point a, Point b)
{
    return Math.Sqrt(DistanceSquared(a, b));
}

public static double DistanceSquared(Point a, Point b)
{
    var d = a - b;
    return DotProduct(d, d);
}

public static double DotProduct(Point a, Point b)
{
    return (a.X * b.X) + (a.Y * b.Y);
}
Mateen Ulhaq
źródło
Próbowałem tego kodu, nie działa całkiem poprawnie. Czasami wydaje się, że pokonuje niewłaściwy dystans.
WDUK
3

Rozwiązanie 2D i 3D

Rozważ zmianę podstawy, tak aby segment linii stał (0, 0, 0)-(d, 0, 0)się punktem (u, v, 0). Najkrótsza odległość występuje w tej płaszczyźnie i jest podana przez

    u ≤ 0 -> d(A, C)
0 ≤ u ≤ d -> |v|
d ≤ u     -> d(B, C)

(odległość do jednego z punktów końcowych lub do linii pomocniczej, w zależności od rzutu na linię. Lokus izo-odległości składa się z dwóch półokręgów i dwóch odcinków linii).

wprowadź opis zdjęcia tutaj

W powyższym wyrażeniu, d jest długością odcinka AB, zaś u, v oznaczają odpowiednio iloczyn skalarny i (moduł) iloczynu AB / d (wektor jednostkowy w kierunku AB) i AC. Stąd wektorowo

AB.AC ≤ 0             -> |AC|
    0 ≤ AB.AC ≤ AB²   -> |ABxAC|/|AB|
          AB² ≤ AB.AC -> |BC|
Yves Daoust
źródło
2

zobacz zestaw narzędzi Matlab GEOMETRY w następującej witrynie: http://people.sc.fsu.edu/~jburkardt/m_src/geometry/geometry.html

ctrl + f i wpisz „segment”, aby znaleźć funkcje związane z segmentem linii. funkcje „segment_point_dist_2d.m” i „segment_point_dist_3d.m” są tym, czego potrzebujesz.

Kody GEOMETRII są dostępne w wersji C i wersji C ++ oraz wersji FORTRAN77 oraz wersji FORTRAN90 i wersji MATLAB.

Lilia
źródło
2

Wersja AutoHotkeys oparta na JavaScript Joshua's:

plDist(x, y, x1, y1, x2, y2) {
    A:= x - x1
    B:= y - y1
    C:= x2 - x1
    D:= y2 - y1

    dot:= A*C + B*D
    sqLen:= C*C + D*D
    param:= dot / sqLen

    if (param < 0 || ((x1 = x2) && (y1 = y2))) {
        xx:= x1
        yy:= y1
    } else if (param > 1) {
        xx:= x2
        yy:= y2
    } else {
        xx:= x1 + param*C
        yy:= y1 + param*D
    }

    dx:= x - xx
    dy:= y - yy

    return sqrt(dx*dx + dy*dy)
}
Chronocide
źródło
2

Nie widziałem tutaj implementacji Java, więc przetłumaczyłem funkcję JavaScript z zaakceptowanej odpowiedzi na kod Java:

static double sqr(double x) {
    return x * x;
}
static double dist2(DoublePoint v, DoublePoint w) {
    return sqr(v.x - w.x) + sqr(v.y - w.y);
}
static double distToSegmentSquared(DoublePoint p, DoublePoint v, DoublePoint w) {
    double l2 = dist2(v, w);
    if (l2 == 0) return dist2(p, v);
    double t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
    if (t < 0) return dist2(p, v);
    if (t > 1) return dist2(p, w);
    return dist2(p, new DoublePoint(
            v.x + t * (w.x - v.x),
            v.y + t * (w.y - v.y)
    ));
}
static double distToSegment(DoublePoint p, DoublePoint v, DoublePoint w) {
    return Math.sqrt(distToSegmentSquared(p, v, w));
}
static class DoublePoint {
    public double x;
    public double y;

    public DoublePoint(double x, double y) {
        this.x = x;
        this.y = y;
    }
}
Jurij Fiodorow
źródło
2

Wersja WPF:

public class LineSegment
{
    private readonly Vector _offset;
    private readonly Vector _vector;

    public LineSegment(Point start, Point end)
    {
        _offset = (Vector)start;
        _vector = (Vector)(end - _offset);
    }

    public double DistanceTo(Point pt)
    {
        var v = (Vector)pt - _offset;

        // first, find a projection point on the segment in parametric form (0..1)
        var p = (v * _vector) / _vector.LengthSquared;

        // and limit it so it lays inside the segment
        p = Math.Min(Math.Max(p, 0), 1);

        // now, find the distance from that point to our point
        return (_vector * p - v).Length;
    }
}
torvina
źródło
1

Oto kod, który skończyłem pisać. Ten kod zakłada, że ​​punkt jest zdefiniowany w postaci {x:5, y:7}. Zauważ, że nie jest to absolutnie najbardziej wydajny sposób, ale jest to najprostszy i najłatwiejszy do zrozumienia kod, jaki mogłem wymyślić.

// a, b, and c in the code below are all points

function distance(a, b)
{
    var dx = a.x - b.x;
    var dy = a.y - b.y;
    return Math.sqrt(dx*dx + dy*dy);
}

function Segment(a, b)
{
    var ab = {
        x: b.x - a.x,
        y: b.y - a.y
    };
    var length = distance(a, b);

    function cross(c) {
        return ab.x * (c.y-a.y) - ab.y * (c.x-a.x);
    };

    this.distanceFrom = function(c) {
        return Math.min(distance(a,c),
                        distance(b,c),
                        Math.abs(cross(c) / length));
    };
}
Eli Courtwright
źródło
1
Ten kod zawiera błąd. Punkt w pobliżu linii, na której leży odcinek, ale daleko od jednego końca odcinka, zostałby błędnie oceniony jako znajdujący się w pobliżu odcinka.
Grumdrig
Interesujące, przyjrzę się temu następnym razem, gdy będę pracować nad tą bazą kodów, aby potwierdzić twoje twierdzenie. Dzięki za wskazówkę.
Eli Courtwright
1

Powyższa funkcja nie działa na liniach pionowych. Oto funkcja, która działa dobrze! Linia z punktami p1, p2. a CheckPoint to p;

public float DistanceOfPointToLine2(PointF p1, PointF p2, PointF p)
{
  //          (y1-y2)x + (x2-x1)y + (x1y2-x2y1)
  //d(P,L) = --------------------------------
  //         sqrt( (x2-x1)pow2 + (y2-y1)pow2 )

  double ch = (p1.Y - p2.Y) * p.X + (p2.X - p1.X) * p.Y + (p1.X * p2.Y - p2.X * p1.Y);
  double del = Math.Sqrt(Math.Pow(p2.X - p1.X, 2) + Math.Pow(p2.Y - p1.Y, 2));
  double d = ch / del;
  return (float)d;
}
Dmitry
źródło
Nie odpowiada na pytanie. Działa to tylko w przypadku linii (tych, które rozciągają się nieskończenie w przestrzeni), a nie segmentów linii (które mają skończoną długość).
Trynidad,
„powyżej funkcji” jest niejednoznacznym odniesieniem. (Mnie irytuje bo czasami ta odpowiedź jest pokazany pod moją odpowiedź.)
RenniePet
1

Oto to samo, co odpowiedź C ++, ale została przeniesiona do pascal. Kolejność parametru punktu zmieniła się, aby pasować do mojego kodu, ale jest taka sama.

function Dot(const p1, p2: PointF): double;
begin
  Result := p1.x * p2.x + p1.y * p2.y;
end;
function SubPoint(const p1, p2: PointF): PointF;
begin
  result.x := p1.x - p2.x;
  result.y := p1.y - p2.y;
end;

function ShortestDistance2(const p,v,w : PointF) : double;
var
  l2,t : double;
  projection,tt: PointF;
begin
  // Return minimum distance between line segment vw and point p
  //l2 := length_squared(v, w);  // i.e. |w-v|^2 -  avoid a sqrt
  l2 := Distance(v,w);
  l2 := MPower(l2,2);
  if (l2 = 0.0) then begin
    result:= Distance(p, v);   // v == w case
    exit;
  end;
  // Consider the line extending the segment, parameterized as v + t (w - v).
  // We find projection of point p onto the line.
  // It falls where t = [(p-v) . (w-v)] / |w-v|^2
  t := Dot(SubPoint(p,v),SubPoint(w,v)) / l2;
  if (t < 0.0) then begin
    result := Distance(p, v);       // Beyond the 'v' end of the segment
    exit;
  end
  else if (t > 1.0) then begin
    result := Distance(p, w);  // Beyond the 'w' end of the segment
    exit;
  end;
  //projection := v + t * (w - v);  // Projection falls on the segment
  tt.x := v.x + t * (w.x - v.x);
  tt.y := v.y + t * (w.y - v.y);
  result := Distance(p, tt);
end;
użytkownik1401452
źródło
Istnieje kilka problemów z tą odpowiedzią: Typ PointF nie jest deklarowany (być może jest to typ standardowy w niektórych implementacjach Pascal). Prawdopodobnie jest to rekord x, y: podwójny; koniec; 2. funkcje Odległość i MPower nie są deklarowane i nie ma wyjaśnienia, co robią (możemy się domyślić, tak). 3. Projekcja zmiennej jest zadeklarowana, ale nigdy nie jest używana. Ogólnie rzecz biorąc, jest to raczej słaba odpowiedź.
dummzeuch