Jak mogę ustalić, czy punkt 2D znajduje się w wielokącie?

Próbuję utworzyć szybki punkt 2D wewnątrz algorytmu wielokąta, do użycia w testowaniu trafień (np Polygon.contains(p:Point).). Docenione zostaną sugestie dotyczące skutecznych technik.

W przypadku grafiki wolałbym nie preferować liczb całkowitych. Wiele systemów używa liczb całkowitych do malowania w interfejsie użytkownika (w końcu piksele są liczbami całkowitymi), ale na przykład macOS używa float do wszystkiego. macOS zna tylko punkty, a punkt może tłumaczyć się na jeden piksel, ale w zależności od rozdzielczości monitora może on przekładać się na coś innego. Na ekranach siatkówki pół punktu (0,5 / 0,5) to piksel. Mimo to nigdy nie zauważyłem, że interfejsy macOS są znacznie wolniejsze niż inne interfejsy. Po tym, jak wszystkie interfejsy API 3D (OpenGL lub Direct3D) współpracują również z pływakami, a nowoczesne biblioteki graficzne bardzo często korzystają z akceleracji GPU.

Teraz powiedziałeś, że prędkość jest twoim głównym zmartwieniem, dobra, chodźmy na szybkość. Zanim uruchomisz jakiś wyrafinowany algorytm, najpierw wykonaj prosty test. Utwórz obwiednię wyrównaną do osi wokół swojego wielokąta. Jest to bardzo łatwe, szybkie i może już zabezpieczyć wiele obliczeń. Jak to działa? Iteruj po wszystkich punktach wielokąta i znajdź wartości min / maks X i Y.

Np. Masz punkty (9/1), (4/3), (2/7), (8/2), (3/6). Oznacza to, że Xmin to 2, Xmax to 9, Ymin to 1, a Ymax to 7. Punkt poza prostokątem o dwóch krawędziach (2/1) i (9/7) nie może znajdować się w wielokącie.

// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
    // Definitely not within the polygon!
}

To pierwszy test, który można przeprowadzić w dowolnym momencie. Jak widać, ten test jest bardzo szybki, ale jest również bardzo szorstki. Aby obsłużyć punkty znajdujące się w obrębie prostokąta ograniczającego, potrzebujemy bardziej zaawansowanego algorytmu. Można to obliczyć na kilka sposobów. Która metoda działa, zależy również od tego, czy wielokąt może mieć otwory lub zawsze będzie solidny. Oto przykłady brył (jeden wypukły, jeden wklęsły):

A oto jeden z dziurą:

Zielony ma dziurę w środku!

Najłatwiejszy algorytm, który może obsłużyć wszystkie trzy powyższe przypadki i nadal jest dość szybki, nazywa się rzutowaniem promieniowym . Algorytm jest dość prosty: narysuj wirtualny promień z dowolnego miejsca poza wielokątem i policz, jak często uderza on w bok wielokąta. Jeśli liczba trafień jest parzysta, znajduje się poza wielokątem, a jeśli jest nieparzysty, to jest w środku.

Algorytm numer uzwojenie byłoby alternatywą, jest bardziej dokładna dla punktów będących bardzo blisko do linii wielokąta, ale jest to również znacznie wolniej. Rzutowanie promieniowe może się nie udać w przypadku punktów znajdujących się zbyt blisko boku wielokąta ze względu na ograniczoną precyzję zmiennoprzecinkową i problemy z zaokrąglaniem, ale w rzeczywistości nie stanowi to problemu, tak jakby punkt znajdował się tak blisko boku, często jest to wizualnie niemożliwe widz rozpoznaje, czy jest już w środku, czy na zewnątrz.

Nadal masz powyższą ramkę, pamiętasz? Wystarczy wybrać punkt poza obwiednią i użyć go jako punktu początkowego dla promienia. Np. Punkt (Xmin - e/p.y)na pewno znajduje się poza wielokątem.

Ale co to jest e? Cóż, e(właściwie epsilon) dodaje obrysowi trochę obicia . Jak powiedziałem, śledzenie promieni nie powiedzie się, jeśli zaczniemy zbyt blisko linii wielokąta. Ponieważ obwiednia może być równa wielobokowi (jeśli wielobok jest prostokątem wyrównanym do osi, obwiednia jest równa samemu wielobokowi!), Potrzebujemy trochę wypełnienia, aby zapewnić bezpieczeństwo, to wszystko. Jak duży powinieneś wybrać e? Nie za duży. To zależy od skali układu współrzędnych używanego do rysowania. Jeśli szerokość kroku w pikselach wynosi 1,0, wybierz po prostu 1,0 (jeszcze 0,1 też by działało)

Teraz, gdy mamy promień ze współrzędnymi początkowymi i końcowymi, problem przesuwa się z „ jest punktem wewnątrz wielokąta ” na „ jak często promień przecina stronę wielokąta ”. Dlatego nie możemy tak po prostu pracować z punktami wielokąta, jak teraz, teraz potrzebujemy rzeczywistych boków. Strona jest zawsze definiowana przez dwa punkty.

side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:

Musisz przetestować promień ze wszystkich stron. Traktuj promień jako wektor, a każdą stronę jako wektor. Promień musi trafić w każdą stronę dokładnie raz lub wcale. Nie może trafić dwukrotnie po tej samej stronie. Dwie linie w przestrzeni 2D zawsze przecinają się dokładnie raz, chyba że są równoległe, w takim przypadku nigdy się nie przecinają. Jednak ponieważ wektory mają ograniczoną długość, dwa wektory mogą nie być równoległe i nadal nigdy się nie przecinają, ponieważ są zbyt krótkie, aby się ze sobą spotkać.

// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
    // Test if current side intersects with ray.
    // If yes, intersections++;
}
if ((intersections & 1) == 1) {
    // Inside of polygon
} else {
    // Outside of polygon
}

Jak dotąd dobrze, ale jak sprawdzić, czy dwa wektory się przecinają? Oto kod C (nie testowany), który powinien załatwić sprawę:

#define NO 0
#define YES 1
#define COLLINEAR 2

int areIntersecting(
    float v1x1, float v1y1, float v1x2, float v1y2,
    float v2x1, float v2y1, float v2x2, float v2y2
) {
    float d1, d2;
    float a1, a2, b1, b2, c1, c2;

    // Convert vector 1 to a line (line 1) of infinite length.
    // We want the line in linear equation standard form: A*x + B*y + C = 0
    // See: http://en.wikipedia.org/wiki/Linear_equation
    a1 = v1y2 - v1y1;
    b1 = v1x1 - v1x2;
    c1 = (v1x2 * v1y1) - (v1x1 * v1y2);

    // Every point (x,y), that solves the equation above, is on the line,
    // every point that does not solve it, is not. The equation will have a
    // positive result if it is on one side of the line and a negative one 
    // if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
    // 2 into the equation above.
    d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
    d2 = (a1 * v2x2) + (b1 * v2y2) + c1;

    // If d1 and d2 both have the same sign, they are both on the same side
    // of our line 1 and in that case no intersection is possible. Careful, 
    // 0 is a special case, that's why we don't test ">=" and "<=", 
    // but "<" and ">".
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // The fact that vector 2 intersected the infinite line 1 above doesn't 
    // mean it also intersects the vector 1. Vector 1 is only a subset of that
    // infinite line 1, so it may have intersected that line before the vector
    // started or after it ended. To know for sure, we have to repeat the
    // the same test the other way round. We start by calculating the 
    // infinite line 2 in linear equation standard form.
    a2 = v2y2 - v2y1;
    b2 = v2x1 - v2x2;
    c2 = (v2x2 * v2y1) - (v2x1 * v2y2);

    // Calculate d1 and d2 again, this time using points of vector 1.
    d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
    d2 = (a2 * v1x2) + (b2 * v1y2) + c2;

    // Again, if both have the same sign (and neither one is 0),
    // no intersection is possible.
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // If we get here, only two possibilities are left. Either the two
    // vectors intersect in exactly one point or they are collinear, which
    // means they intersect in any number of points from zero to infinite.
    if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;

    // If they are not collinear, they must intersect in exactly one point.
    return YES;
}

Wartości wejściowe to dwa punkty końcowe wektora 1 ( v1x1/v1y1i v1x2/v1y2) i wektora 2 ( v2x1/v2y1i v2x2/v2y2). Masz więc 2 wektory, 4 punkty, 8 współrzędnych. YESi NOsą jasne. YESzwiększa skrzyżowania, NOnic nie robi.

Co z COLLINEAR? Oznacza to, że oba wektory leżą na tej samej nieskończonej linii, w zależności od położenia i długości, nie przecinają się wcale lub przecinają się w nieskończonej liczbie punktów. Nie jestem absolutnie pewien, jak poradzić sobie z tą sprawą, nie liczyłbym jej jako skrzyżowania. Cóż, ten przypadek i tak jest raczej rzadki w praktyce z powodu błędów zaokrąglania zmiennoprzecinkowego; Lepszy kod prawdopodobnie nie przetestowałby, == 0.0fale zamiast czegoś takiego < epsilon, gdzie epsilon jest raczej małą liczbą.

Jeśli chcesz przetestować większą liczbę punktów, z pewnością możesz trochę przyspieszyć, utrzymując w pamięci standardowe równania liniowe boków wielokąta, więc nie musisz ich ponownie obliczać za każdym razem. Pozwoli to zaoszczędzić dwa mnożenia zmiennoprzecinkowe i trzy odejmowania zmiennoprzecinkowe na każdym teście w zamian za przechowywanie w pamięci trzech wartości zmiennoprzecinkowych na stronę wielokąta. Jest to typowa kompromis między pamięcią a czasem obliczeniowym.

Na koniec: jeśli możesz użyć sprzętu 3D do rozwiązania problemu, istnieje interesująca alternatywa. Po prostu pozwól GPU wykonać całą pracę za Ciebie. Utwórz malowaną powierzchnię, która jest poza ekranem. Wypełnij go całkowicie kolorem czarnym. Teraz pozwól OpenGL lub Direct3D pomalować twój wielokąt (lub nawet wszystkie twoje wielokąty, jeśli chcesz tylko sprawdzić, czy punkt znajduje się w którymkolwiek z nich, ale nie zależy ci na którym) i wypełnij wielokąt (y) innym kolor, np. biały. Aby sprawdzić, czy punkt znajduje się w wielokącie, uzyskaj kolor tego punktu z powierzchni rysowania. To tylko pobieranie pamięci O (1).

Oczywiście ta metoda jest użyteczna tylko wtedy, gdy twoja powierzchnia do rysowania nie musi być duża. Jeśli nie może zmieścić się w pamięci GPU, ta metoda jest wolniejsza niż w przypadku procesora. Gdyby musiała być ogromna, a Twój procesor graficzny obsługuje nowoczesne moduły cieniujące, nadal możesz używać procesora graficznego, implementując rzutowanie promieni pokazane powyżej jako moduł cieniujący GPU, co jest absolutnie możliwe. W przypadku większej liczby wielokątów lub dużej liczby punktów do przetestowania, opłaci się to, biorąc pod uwagę, że niektóre procesory graficzne będą w stanie przetestować od 64 do 256 punktów równolegle. Należy jednak pamiętać, że przesyłanie danych z procesora do GPU iz powrotem jest zawsze drogie, więc po prostu testowanie kilku punktów na kilku prostych wielokątach, w których albo punkty lub wielokąty są dynamiczne i często się zmieniają, podejście GPU rzadko płaci poza.

Myślę, że następujący fragment kodu jest najlepszym rozwiązaniem (wzięty stąd ):

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
     (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
  }
  return c;
}

Argumenty

• nvert : liczba wierzchołków w wielokącie. To, czy powtórzyć pierwszy wierzchołek na końcu, zostało omówione w powyższym artykule.
• vertx, verty : Tablice zawierające współrzędne x i y wierzchołków wielokąta.
• testx, testy : współrzędna X i Y punktu testowego.

Jest zarówno krótki, jak i wydajny i działa zarówno na wypukłe, jak i wklęsłe wielokąty. Jak zasugerowano wcześniej, należy najpierw sprawdzić prostokąt obwiedni i osobno potraktować otwory wielokąta.

Idea tego jest dość prosta. Autor opisuje to w następujący sposób:

Wypuszczam częściowo nieskończony promień w poziomie (zwiększając x, ustalone y) poza punkt testowy i liczę, ile krawędzi przecina. Na każdym skrzyżowaniu promień przełącza się między wnętrzem a zewnętrzem. Nazywa się to twierdzeniem krzywej Jordana.

Zmienna c zmienia się z 0 na 1 i 1 na 0 za każdym razem, gdy promień poziomy przecina dowolną krawędź. Zasadniczo śledzi więc, czy liczba skrzyżowanych krawędzi jest parzysta czy nieparzysta. 0 oznacza parzystą, a 1 nieparzystą.

Oto wersja C # odpowiedzi udzielonej przez nirg , która pochodzi od tego profesora RPI . Pamiętaj, że użycie kodu z tego źródła RPI wymaga przypisania.

U góry dodano obwiednię. Jednak, jak podkreśla James Brown, główny kod jest prawie tak szybki jak samo sprawdzenie pola granicznego, więc sprawdzenie pola granicznego może faktycznie spowolnić całą operację, w przypadku gdy większość punktów, które sprawdzasz, znajduje się wewnątrz ramki granicznej . Możesz więc pozostawić pole ograniczające wyewidencjonowane, lub alternatywnie byłoby wstępnie obliczyć obwiednie twoich wielokątów, jeśli nie zmieniają one kształtu zbyt często.

public bool IsPointInPolygon( Point p, Point[] polygon )
{
    double minX = polygon[ 0 ].X;
    double maxX = polygon[ 0 ].X;
    double minY = polygon[ 0 ].Y;
    double maxY = polygon[ 0 ].Y;
    for ( int i = 1 ; i < polygon.Length ; i++ )
    {
        Point q = polygon[ i ];
        minX = Math.Min( q.X, minX );
        maxX = Math.Max( q.X, maxX );
        minY = Math.Min( q.Y, minY );
        maxY = Math.Max( q.Y, maxY );
    }

    if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
    {
        return false;
    }

    // https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
    bool inside = false;
    for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
    {
        if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
             p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
        {
            inside = !inside;
        }
    }

    return inside;
}

Oto wariant JavaScript odpowiedzi M. Katza oparty na podejściu Nirga:

function pointIsInPoly(p, polygon) {
    var isInside = false;
    var minX = polygon[0].x, maxX = polygon[0].x;
    var minY = polygon[0].y, maxY = polygon[0].y;
    for (var n = 1; n < polygon.length; n++) {
        var q = polygon[n];
        minX = Math.min(q.x, minX);
        maxX = Math.max(q.x, maxX);
        minY = Math.min(q.y, minY);
        maxY = Math.max(q.y, maxY);
    }

    if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
        return false;
    }

    var i = 0, j = polygon.length - 1;
    for (i, j; i < polygon.length; j = i++) {
        if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
                p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
            isInside = !isInside;
        }
    }

    return isInside;
}

Oblicz zorientowaną sumę kątów między punktem p a każdym wierzchołkiem wielokąta. Jeśli całkowity zorientowany kąt wynosi 360 stopni, punkt znajduje się wewnątrz. Jeśli suma wynosi 0, punkt jest na zewnątrz.

Bardziej podoba mi się ta metoda, ponieważ jest bardziej niezawodna i mniej zależna od precyzji numerycznej.

Metody obliczania równości liczby skrzyżowań są ograniczone, ponieważ można „trafić” wierzchołek podczas obliczania liczby skrzyżowań.

EDYCJA: Przy okazji, ta metoda działa z wklęsłymi i wypukłymi wielokątami.

EDYCJA: Niedawno znalazłem cały artykuł w Wikipedii na ten temat.

To pytanie jest bardzo interesujące. Mam inny praktyczny pomysł, inny niż inne odpowiedzi na ten post. Chodzi o to, aby użyć sumy kątów, aby zdecydować, czy cel jest w środku, czy na zewnątrz. Lepiej znany jako liczba uzwojenia .

Niech x będzie punktem docelowym. Niech tablica [0, 1, .... n] będzie wszystkimi punktami obszaru. Połącz punkt docelowy z każdym punktem granicznym za pomocą linii. Jeśli punkt docelowy znajduje się w tym obszarze. Suma wszystkich kątów wyniesie 360 ​​stopni. Jeśli nie, kąty będą mniejsze niż 360.

Zapoznaj się z tym obrazem, aby uzyskać podstawowe zrozumienie tego pomysłu:

Mój algorytm zakłada, że ​​kierunek zgodny z ruchem wskazówek zegara to kierunek dodatni. Oto potencjalny wkład:

[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]

Poniżej znajduje się kod python, który implementuje ten pomysł:

def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
    a = border[i]
    b = border[i + 1]

    # calculate distance of vector
    A = getDistance(a[0], a[1], b[0], b[1]);
    B = getDistance(target[0], target[1], a[0], a[1])
    C = getDistance(target[0], target[1], b[0], b[1])

    # calculate direction of vector
    ta_x = a[0] - target[0]
    ta_y = a[1] - target[1]
    tb_x = b[0] - target[0]
    tb_y = b[1] - target[1]

    cross = tb_y * ta_x - tb_x * ta_y
    clockwise = cross < 0

    # calculate sum of angles
    if(clockwise):
        degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
    else:
        degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))

if(abs(round(degree) - 360) <= 3):
    return True
return False

Artykuł Erica Hainesa cytowany przez bobobobo jest naprawdę doskonały. Szczególnie interesujące są tabele porównujące wydajność algorytmów; metoda sumowania kątów jest naprawdę zła w porównaniu z innymi. Interesujące jest również to, że takie optymalizacje, jak użycie siatki odnośników do dalszego podziału wielokąta na sektory „wejściowy” i „wyjściowy”, mogą sprawić, że test będzie niesamowicie szybki nawet na wielokątach o> 1000 bokach.

W każdym razie, są wczesne dni, ale mój głos dotyczy metody „przejazdów”, co, jak sądzę, jest prawie tym, co opisuje Mecki. Jednak najbardziej trafnie go opisałem i skodyfikowałem przez Davida Bourke . Uwielbiam to, że nie jest wymagana prawdziwa trygonometria i działa ona na wypukłe i wklęsłe i działa dość dobrze, gdy rośnie liczba boków.

Nawiasem mówiąc, oto jedna z tabel wydajności z artykułu Erica Hainesa dla zainteresowania, testowania losowych wielokątów.

                       number of edges per polygon
                         3       4      10      100    1000
MacMartin               2.9     3.2     5.9     50.6    485
Crossings               3.1     3.4     6.8     60.0    624
Triangle Fan+edge sort  1.1     1.8     6.5     77.6    787
Triangle Fan            1.2     2.1     7.3     85.4    865
Barycentric             2.1     3.8    13.8    160.7   1665
Angle Summation        56.2    70.4   153.6   1403.8  14693

Grid (100x100)          1.5     1.5     1.6      2.1      9.8
Grid (20x20)            1.7     1.7     1.9      5.7     42.2
Bins (100)              1.8     1.9     2.7     15.1    117
Bins (20)               2.1     2.2     3.7     26.3    278

Szybka wersja odpowiedzi nirg :

extension CGPoint {
    func isInsidePolygon(vertices: [CGPoint]) -> Bool {
        guard !vertices.isEmpty else { return false }
        var j = vertices.last!, c = false
        for i in vertices {
            let a = (i.y > y) != (j.y > y)
            let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
            if a && b { c = !c }
            j = i
        }
        return c
    }
}

Naprawdę podoba się rozwiązanie opublikowane przez Nirga i zredagowane przez bobobobo. Właśnie uczyniłem go javascript przyjaznym i trochę bardziej czytelnym dla mojego zastosowania:

function insidePoly(poly, pointx, pointy) {
    var i, j;
    var inside = false;
    for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
        if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
    }
    return inside;
}

Pracowałem nad tym, kiedy byłem badaczem pod kierunkiem Michaela Stonebrakera - wiesz, profesor, który wymyślił Ingres , PostgreSQL itp.

Uświadomiliśmy sobie, że najszybszym sposobem było najpierw wykonać obwiednię, ponieważ jest ona superszybka. Jeśli jest poza obwiednią, jest na zewnątrz. W przeciwnym razie wykonasz cięższą pracę ...

Jeśli chcesz świetnego algorytmu, zajrzyj do kodu źródłowego PostgreSQL projektu open source do pracy geo ...

Chcę podkreślić, że nigdy nie mieliśmy wglądu w prawą i lewą rękę (wyrażalną również jako problem „wewnątrz” vs. „na zewnątrz” ...

AKTUALIZACJA

Link BKB dostarczył sporo rozsądnych algorytmów. Pracowałem nad problemami nauk o Ziemi i dlatego potrzebowałem rozwiązania, które działa na szerokości i długości geograficznej, i ma szczególny problem z poręcznością - czy obszar jest wewnątrz mniejszego lub większego obszaru? Odpowiedź jest taka, że ​​„kierunek” wierzchołków ma znaczenie - jest leworęczny lub praworęczny, w ten sposób możesz wskazać dowolny obszar jako „wewnątrz” dowolnego wielokąta. W związku z tym w mojej pracy wykorzystałem rozwiązanie trzy wymienione na tej stronie.

Ponadto w mojej pracy wykorzystałem osobne funkcje do testów „on-line”.

... Ponieważ ktoś zapytał: ustaliliśmy, że testy ramki ograniczającej były najlepsze, gdy liczba pionów przekroczyła pewną liczbę - wykonaj bardzo szybki test przed wykonaniem dłuższego testu, jeśli to konieczne ... Ramkę ograniczającą tworzy się po prostu biorąc największy x, najmniejszy x, największy y i najmniejszy y i zebranie ich razem, aby utworzyć cztery punkty pudełka ...

Kolejna wskazówka dla następnych: wykonaliśmy wszystkie nasze bardziej wyrafinowane i „ściemniające” obliczenia w przestrzeni siatki wszystko w dodatnich punktach na płaszczyźnie, a następnie ponownie rzutowaliśmy na „rzeczywistą” długość / szerokość geograficzną, unikając w ten sposób możliwych błędów owijanie się, gdy jedna linia przekroczy 180 długości geograficznej i podczas obchodzenia się z regionami polarnymi. Działa świetnie!

Odpowiedź Davida Segonda jest w zasadzie standardową odpowiedzią ogólną, a Richard T jest najczęstszą optymalizacją, choć są też inne. Inne silne optymalizacje oparte są na mniej ogólnych rozwiązaniach. Na przykład, jeśli zamierzasz sprawdzić ten sam wielokąt z dużą ilością punktów, triangulacja wielokąta może znacznie przyspieszyć sprawę, ponieważ istnieje wiele bardzo szybkich algorytmów wyszukiwania TIN. Innym jest, jeśli wielokąt i punkty znajdują się na ograniczonej płaszczyźnie w niskiej rozdzielczości, powiedzmy na ekranie, możesz pomalować wielokąt na buforze wyświetlania odwzorowanym w pamięci w danym kolorze i sprawdzić kolor danego piksela, aby zobaczyć, czy leży w wielokątach.

Podobnie jak wiele optymalizacji, opierają się one na konkretnych, a nie ogólnych przypadkach, i dają korzyści na podstawie zamortyzowanego czasu, a nie pojedynczego użycia.

Pracując w tej dziedzinie, znalazłem Joesepha O'Rourkesa „Geometria obliczeniowa w C” ISBN 0-521-44034-3, która była wielką pomocą.

Trywialnym rozwiązaniem byłoby podzielenie wielokąta na trójkąty i przetestowanie trójkątów, jak wyjaśniono tutaj

Jeśli twój wielokąt jest CONVEX , może być lepsze podejście. Spójrz na wielokąt jako zbiór nieskończonych linii. Każda linia dzieli przestrzeń na dwie części. dla każdego punktu łatwo jest powiedzieć, czy jest po jednej, czy po drugiej stronie linii. Jeśli punkt znajduje się po tej samej stronie wszystkich linii, to znajduje się wewnątrz wielokąta.

Zdaję sobie sprawę, że to jest stare, ale tutaj jest algorytm rzucania promieni zaimplementowany w Cocoa, na wypadek, gdyby ktoś był zainteresowany. Nie jestem pewien, czy jest to najskuteczniejszy sposób robienia rzeczy, ale może komuś pomóc.

- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
    NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
    BOOL result;
    float aggregateX = 0; //I use these to calculate the centroid of the shape
    float aggregateY = 0;
    NSPoint firstPoint[1];
    [currentPath elementAtIndex:0 associatedPoints:firstPoint];
    float olderX = firstPoint[0].x;
    float olderY = firstPoint[0].y;
    NSPoint interPoint;
    int noOfIntersections = 0;

    for (int n = 0; n < [currentPath elementCount]; n++) {
        NSPoint points[1];
        [currentPath elementAtIndex:n associatedPoints:points];
        aggregateX += points[0].x;
        aggregateY += points[0].y;
    }

    for (int n = 0; n < [currentPath elementCount]; n++) {
        NSPoint points[1];

        [currentPath elementAtIndex:n associatedPoints:points];
        //line equations in Ax + By = C form
        float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;  
        float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
        float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);

        float _A_BAR = olderY - points[0].y;
        float _B_BAR = points[0].x - olderX;
        float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);

        float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
        if (det != 0) {
            //intersection points with the edges
            float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
            float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
            interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
            if (olderX <= points[0].x) {
                //doesn't matter in which direction the ray goes, so I send it right-ward.
                if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {  
                    noOfIntersections++;
                }
            } else {
                if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
                     noOfIntersections++;
                } 
            }
        }
        olderX = points[0].x;
        olderY = points[0].y;
    }
    if (noOfIntersections % 2 == 0) {
        result = FALSE;
    } else {
        result = TRUE;
    }
    return result;
}

Wersja Obj-C odpowiedzi nirga z przykładową metodą testowania punktów. Odpowiedź Nirga działała dla mnie dobrze.

- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
    NSUInteger nvert = [vertices count];
    NSInteger i, j, c = 0;
    CGPoint verti, vertj;

    for (i = 0, j = nvert-1; i < nvert; j = i++) {
        verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
        vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
        if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
        ( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
            c = !c;
    }

    return (c ? YES : NO);
}

- (void)testPoint {

    NSArray *polygonVertices = [NSArray arrayWithObjects:
        [NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
        [NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
        [NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
        [NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
        [NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
        [NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
        nil
    ];

    CGPoint tappedPoint = CGPointMake(23.0, 70.0);

    if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
        NSLog(@"YES");
    } else {
        NSLog(@"NO");
    }
}

Nie ma nic piękniejszego niż indukcyjna definicja problemu. Dla kompletności tutaj masz wersję prologu, która może również wyjaśnić kwestie leżące u podstaw rzucania promieni :

Na podstawie symulacji algorytmu prostoty w http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Niektóre przewidywania pomocnika:

exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).

inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) +      X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).

get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).

Równanie linii z 2 punktami A i B (Linia (A, B)) jest następujące:

                    (YB-YA)
           Y - YA = ------- * (X - XA) 
                    (XB-YB) 

Ważne jest, aby kierunek obrotu linii był ustawiony zgodnie z ruchem wskazówek zegara dla granic i przeciwnie do zegara dla otworów. Sprawdzimy, czy punkt (X, Y), tj. Testowany punkt znajduje się w lewej półpłaszczyźnie naszej linii (jest to kwestia gustu, może to być również prawa strona, ale także kierunek granic) linie muszą zostać zmienione w tym przypadku), ma to na celu rzutowanie promienia z punktu w prawo (lub w lewo) i potwierdzenie przecięcia z linią. Zdecydowaliśmy się wyświetlać promień w kierunku poziomym (znowu jest to kwestia gustu, można to również zrobić w pionie z podobnymi ograniczeniami), więc mamy:

               (XB-XA)
           X < ------- * (Y - YA) + XA
               (YB-YA) 

Teraz musimy wiedzieć, czy punkt znajduje się tylko po lewej (lub po prawej) stronie segmentu linii, a nie na całej płaszczyźnie, więc musimy ograniczyć wyszukiwanie tylko do tego segmentu, ale jest to łatwe, ponieważ można znaleźć się w tym segmencie tylko jeden punkt na linii może znajdować się wyżej niż Y w osi pionowej. Ponieważ jest to silniejsze ograniczenie, najpierw musi to sprawdzić, dlatego najpierw bierzemy tylko te linie, które spełniają ten wymóg, a następnie sprawdzamy jego możliwości. Zgodnie z twierdzeniem Curve'a Jordan promień rzutowany na wielokąt musi przecinać się w parzystej liczbie linii. Więc skończymy, rzucimy promień w prawo, a następnie za każdym razem, gdy przecina on linię, przełącza swój stan. Jednak w naszym wdrożeniu sprawdzimy długość pakietu rozwiązań spełniających podane ograniczenia i zdecydujemy o tym. należy to zrobić dla każdej linii w wielokącie.

is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] =  [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA)); 
                                                        is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).

in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon),  in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line),    in_y_range_at_poly(Coordinate,Line,Polygon), Lines).

traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).

% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).

Wersja C # odpowiedzi nirg jest tutaj: po prostu podzielę się kodem. Może komuś zaoszczędzić trochę czasu.

public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) {
            bool result = false;
            int j = polygon.Count() - 1;
            for (int i = 0; i < polygon.Count(); i++) {
                if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) {
                    if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) {
                        result = !result;
                    }
                }
                j = i;
            }
            return result;
        }

Wersja Java:

public class Geocode {
    private float latitude;
    private float longitude;

    public Geocode() {
    }

    public Geocode(float latitude, float longitude) {
        this.latitude = latitude;
        this.longitude = longitude;
    }

    public float getLatitude() {
        return latitude;
    }

    public void setLatitude(float latitude) {
        this.latitude = latitude;
    }

    public float getLongitude() {
        return longitude;
    }

    public void setLongitude(float longitude) {
        this.longitude = longitude;
    }
}

public class GeoPolygon {
    private ArrayList<Geocode> points;

    public GeoPolygon() {
        this.points = new ArrayList<Geocode>();
    }

    public GeoPolygon(ArrayList<Geocode> points) {
        this.points = points;
    }

    public GeoPolygon add(Geocode geo) {
        points.add(geo);
        return this;
    }

    public boolean inside(Geocode geo) {
        int i, j;
        boolean c = false;
        for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
            if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
                    (geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
                c = !c;
        }
        return c;
    }

}

Port .Net:

    static void Main(string[] args)
    {

        Console.Write("Hola");
        List<double> vertx = new List<double>();
        List<double> verty = new List<double>();

        int i, j, c = 0;

        vertx.Add(1);
        vertx.Add(2);
        vertx.Add(1);
        vertx.Add(4);
        vertx.Add(4);
        vertx.Add(1);

        verty.Add(1);
        verty.Add(2);
        verty.Add(4);
        verty.Add(4);
        verty.Add(1);
        verty.Add(1);

        int nvert = 6;  //Vértices del poligono

        double testx = 2;
        double testy = 5;


        for (i = 0, j = nvert - 1; i < nvert; j = i++)
        {
            if (((verty[i] > testy) != (verty[j] > testy)) &&
             (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
                c = 1;
        }
    }

WERSJA VBA:

Uwaga: pamiętaj, że jeśli twój wielokąt jest obszarem na mapie, szerokość i długość geograficzna są wartościami Y / X w przeciwieństwie do X / Y (szerokość geograficzna = Y, długość geograficzna = X), ponieważ z tego, co rozumiem, są implikacje historyczne z czasów, gdy Długość geograficzna nie była miarą.

MODUŁ KLASOWY: CPoint

Private pXValue As Double
Private pYValue As Double

'''''X Value Property'''''

Public Property Get X() As Double
    X = pXValue
End Property

Public Property Let X(Value As Double)
    pXValue = Value
End Property

'''''Y Value Property'''''

Public Property Get Y() As Double
    Y = pYValue
End Property

Public Property Let Y(Value As Double)
    pYValue = Value
End Property

MODUŁ:

Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean

    Dim i As Integer
    Dim j As Integer
    Dim q As Object
    Dim minX As Double
    Dim maxX As Double
    Dim minY As Double
    Dim maxY As Double
    minX = polygon(0).X
    maxX = polygon(0).X
    minY = polygon(0).Y
    maxY = polygon(0).Y

    For i = 1 To UBound(polygon)
        Set q = polygon(i)
        minX = vbMin(q.X, minX)
        maxX = vbMax(q.X, maxX)
        minY = vbMin(q.Y, minY)
        maxY = vbMax(q.Y, maxY)
    Next i

    If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
        isPointInPolygon = False
        Exit Function
    End If


    ' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

    isPointInPolygon = False
    i = 0
    j = UBound(polygon)

    Do While i < UBound(polygon) + 1
        If (polygon(i).Y > p.Y) Then
            If (polygon(j).Y < p.Y) Then
                If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
                    isPointInPolygon = True
                    Exit Function
                End If
            End If
        ElseIf (polygon(i).Y < p.Y) Then
            If (polygon(j).Y > p.Y) Then
                If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
                    isPointInPolygon = True
                    Exit Function
                End If
            End If
        End If
        j = i
        i = i + 1
    Loop   
End Function

Function vbMax(n1, n2) As Double
    vbMax = IIf(n1 > n2, n1, n2)
End Function

Function vbMin(n1, n2) As Double
    vbMin = IIf(n1 > n2, n2, n1)
End Function


Sub TestPointInPolygon()

    Dim i As Integer
    Dim InPolygon As Boolean

'   MARKER Object
    Dim p As CPoint
    Set p = New CPoint
    p.X = <ENTER X VALUE HERE>
    p.Y = <ENTER Y VALUE HERE>

'   POLYGON OBJECT
    Dim polygon() As CPoint
    ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
    For i = 0 To <ENTER VALUE HERE> 'Same value as above
       Set polygon(i) = New CPoint
       polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
       polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
    Next i

    InPolygon = isPointInPolygon(p, polygon)
    MsgBox InPolygon

End Sub

Zrobiłem wdrożenie Pythona z nirg w C ++ kod :

Wejścia

• bounding_points: węzły tworzące wielokąt.

• bounding_box_positions: punkty kandydujące do filtrowania. (W mojej implementacji utworzonej z ramki granicznej.

  (Wejścia są wykazy krotek w formacie: [(xcord, ycord), ...])

Zwroty

• Wszystkie punkty wewnątrz wielokąta.

def polygon_ray_casting(self, bounding_points, bounding_box_positions):
    # Arrays containing the x- and y-coordinates of the polygon's vertices.
    vertx = [point[0] for point in bounding_points]
    verty = [point[1] for point in bounding_points]
    # Number of vertices in the polygon
    nvert = len(bounding_points)
    # Points that are inside
    points_inside = []

    # For every candidate position within the bounding box
    for idx, pos in enumerate(bounding_box_positions):
        testx, testy = (pos[0], pos[1])
        c = 0
        for i in range(0, nvert):
            j = i - 1 if i != 0 else nvert - 1
            if( ((verty[i] > testy ) != (verty[j] > testy))   and
                    (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
                c += 1
        # If odd, that means that we are inside the polygon
        if c % 2 == 1: 
            points_inside.append(pos)


    return points_inside

Ponownie pomysł zaczerpnięto stąd

Zaskoczony, że nikt wcześniej tego nie poruszał, ale dla pragmatyków wymagających bazy danych: MongoDB ma doskonałe wsparcie dla zapytań Geo, w tym tego.

To czego szukasz to:

db.neighborhoods.findOne ({geometria: {$ geoIntersects: {$ geometry: {type: „Point”, współrzędne: [„długość geograficzna”, „szerokość geograficzna”]}}}})

Neighborhoodsto kolekcja przechowująca jeden lub więcej wielokątów w standardowym formacie GeoJson. Jeśli zapytanie zwraca null, nie jest przecinane, inaczej jest.

Bardzo dobrze udokumentowane tutaj: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

Wydajność dla ponad 6000 punktów sklasyfikowanych w 330 nieregularnych siatkach wielokątów była krótsza niż jedna minuta bez żadnej optymalizacji i obejmowała czas na aktualizację dokumentów za pomocą odpowiedniego wielokąta.

Oto punkt w teście wielokąta w C, który nie używa rzutowania promieniami. I może działać dla nakładających się obszarów (przecięcia siebie), patrz use_holesargument.

/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);

/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
                         const bool use_holes)
{
    /* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
    float angletot = 0.0;
    float fp1[2], fp2[2];
    unsigned int i;
    const float *p1, *p2;

    p1 = verts[nr - 1];

    /* first vector */
    fp1[0] = p1[0] - pt[0];
    fp1[1] = p1[1] - pt[1];

    for (i = 0; i < nr; i++) {
        p2 = verts[i];

        /* second vector */
        fp2[0] = p2[0] - pt[0];
        fp2[1] = p2[1] - pt[1];

        /* dot and angle and cross */
        angletot += angle_signed_v2v2(fp1, fp2);

        /* circulate */
        copy_v2_v2(fp1, fp2);
        p1 = p2;
    }

    angletot = fabsf(angletot);
    if (use_holes) {
        const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
        angletot -= nested * (float)(M_PI * 2.0);
        return (angletot > 4.0f) != ((int)nested % 2);
    }
    else {
        return (angletot > 4.0f);
    }
}

/* math lib */

static float dot_v2v2(const float a[2], const float b[2])
{
    return a[0] * b[0] + a[1] * b[1];
}

static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
    const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
    return atan2f(perp_dot, dot_v2v2(v1, v2));
}

static void copy_v2_v2(float r[2], const float a[2])
{
    r[0] = a[0];
    r[1] = a[1];
}

Uwaga: jest to jedna z mniej optymalnych metod, ponieważ zawiera wiele wywołań atan2f, ale może zainteresować programistów czytających ten wątek (w moich testach jest ~ 23 razy wolniejszy niż metoda przecięcia linii).

Aby wykryć trafienie w wielokąt, musimy przetestować dwie rzeczy:

1. Jeśli punkt znajduje się w obszarze wielokąta. (można to osiągnąć za pomocą algorytmu Ray-Casting)
2. Jeśli punkt znajduje się na granicy wielokąta (można to zrobić za pomocą tego samego algorytmu, który jest używany do wykrywania punktu na polilinii (linii)).

Aby poradzić sobie z następującymi szczególnymi przypadkami w algorytmie rzutowania promieni :

1. Promień zachodzi na jedną ze stron wielokąta.
2. Punkt znajduje się wewnątrz wielokąta, a promień przechodzi przez wierzchołek wielokąta.
3. Punkt znajduje się poza wielokątem, a promień dotyka tylko jednego kąta wielokąta.

Sprawdź, czy punkt znajduje się w złożonym wielokącie . Artykuł zapewnia łatwy sposób ich rozwiązania, więc w powyższych przypadkach nie będzie wymagane specjalne leczenie.

Można to zrobić, sprawdzając, czy obszar utworzony przez połączenie pożądanego punktu z wierzchołkami wielokąta odpowiada obszarowi samego wielokąta.

Lub możesz sprawdzić, czy suma kątów wewnętrznych od twojego punktu do każdej pary dwóch kolejnych wierzchołków wielokąta do sumy punktów kontrolnych do 360, ale mam wrażenie, że pierwsza opcja jest szybsza, ponieważ nie wymaga podziałów ani obliczeń odwrotności funkcji trygonometrycznych.

Nie wiem, co się stanie, jeśli twój wielokąt ma dziurę w środku, ale wydaje mi się, że główny pomysł można dostosować do tej sytuacji

Możesz również opublikować pytanie w społeczności matematyki. Założę się, że mają na to milion sposobów

Jeśli szukasz biblioteki skryptów Java, istnieje rozszerzenie javascript google maps v3 dla klasy Polygon, aby wykryć, czy w nim znajduje się jakiś punkt.

var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);

Rozszerzenie Google Github

Podczas używania qt(Qt 4.3+), można użyć funkcji QPolygon za containsPoint

Odpowiedź zależy od tego, czy masz proste, czy złożone wielokąty. Proste wielokąty nie mogą mieć żadnych przecięć segmentów linii. Mogą mieć dziury, ale linie nie mogą się przecinać. Złożone regiony mogą mieć przecięcia linii - mogą więc mieć zachodzące na siebie regiony lub regiony, które stykają się ze sobą tylko jednym punktem.

W przypadku prostych wielokątów najlepszym algorytmem jest algorytm rzutowania promieniowego (liczba przecinająca). W przypadku złożonych wielokątów ten algorytm nie wykrywa punktów znajdujących się w nakładających się regionach. Tak więc w przypadku złożonych wielokątów należy użyć algorytmu liczby uzwojenia.

Oto doskonały artykuł z implementacją C obu algorytmów. Próbowałem ich i działają dobrze.

http://geomalgorithms.com/a03-_inclusion.html

Wersja rozwiązania Scala firmy nirg (zakłada, że ​​wstępne sprawdzenie prostokąta ograniczającego odbywa się osobno):

def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {

  val length = polygon.length

  @tailrec
  def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
    if (i == length)
      tracker
    else {
      val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
      oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
    }
  }

  oddIntersections(0, length - 1, tracker = false)
}

Oto wersja golang odpowiedzi @nirg (zainspirowana kodem C # przez @@ m-katz)

func isPointInPolygon(polygon []point, testp point) bool {
    minX := polygon[0].X
    maxX := polygon[0].X
    minY := polygon[0].Y
    maxY := polygon[0].Y

    for _, p := range polygon {
        minX = min(p.X, minX)
        maxX = max(p.X, maxX)
        minY = min(p.Y, minY)
        maxY = max(p.Y, maxY)
    }

    if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
        return false
    }

    inside := false
    j := len(polygon) - 1
    for i := 0; i < len(polygon); i++ {
        if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
            inside = !inside
        }
        j = i
    }

    return inside
}