Jak Mogę określić, czy punkt 2D znajduje się w Wielokątu?

Myślę, że poniższy fragment kodu jest najlepszym rozwiązaniem (zaczerpniętym z Tutaj):

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ątu. To, czy powtórzyć pierwszy wierzchołek na końcu, zostało omówione w artykule, o którym mowa powyżej.
• 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.

To i to krótki i wydajny i działa zarówno dla wypukłych, jak i wklęsłych wielokątów. Jak sugerowano wcześniej, należy najpierw sprawdzić obwiedniowy prostokąt i traktować otwory wielokątów oddzielnie.

Uruchamiam pół-nieskończony promień poziomo (zwiększając x, ustalając y) z punktu testowego i policzę, ile krawędzi przecina. Na każdym przejściu promień przełącza się między wewnątrz i na zewnątrz. To się nazywa Krzywa Jordana twierdzenie.

Zmienna c zmienia się z 0 Na 1 i 1 na 0 za każdym razem, gdy promień poziomy przekracza dowolną krawędź. W zasadzie to śledzenie, czy liczba skrzyżowanych krawędzi jest parzysta, czy nieparzysta. 0 oznacza parzyste, a 1 Nieparzyste.

Oto Wersja C# odpowiedzi udzielonej przez nirg, która pochodzi odtego profesora RPI . Zauważ, że użycie kodu z tego źródła RPI wymaga atrybucji.

Na górze dodano zaznaczenie obwiedni. Jednak, jak zauważa James Brown, kod główny jest prawie tak szybki jak samo zaznaczenie obwiedni, więc może faktycznie spowolnić ogólną operację, w przypadku gdy większość zaznaczanych punktów znajduje się wewnątrz obwiedni. Więc możesz opuścić obwiednię sprawdzoną lub alternatywą byłoby wcześniejsze obliczenie obwiedni 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;
    }

    // http://www.ecse.rpi.edu/Homepages/wrf/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 bazujący na podejściu Nirg:

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 z wierzchołków wielokąta. Jeśli całkowity kąt zorientowany wynosi 360 stopni, punkt znajduje się wewnątrz. Jeśli suma jest równa 0, punkt jest na zewnątrz.

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

Metody obliczające parzystość liczby przecięć są ograniczone, ponieważ można "trafić" wierzchołek podczas obliczania liczby przecięć.

EDIT: By The Sposób, ta metoda działa z wklęsłymi i wypukłymi wielokątami.

EDIT: znalazłem Ostatnio cały Artykuł Wikipedii na ten temat.

Eric Haines Artykuł 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 optymalizacje, takie jak użycie siatki wyszukiwania do dalszego podziału wielokąta na sektory "wejście" i "wyjście", mogą sprawić, że test będzie niewiarygodnie szybki nawet na wielokątach o > 1000 stronach.

W każdym razie, jest wcześnie, ale mój głos idzie do " crossings" metoda, czyli chyba to, co Mecki opisuje. Jednak znalazłem go najbardziej zwięźle opisany i skodyfikowany przez Davida Bourke . Podoba mi się to, że nie jest wymagana prawdziwa Trygonometria, a to działa na wypukłe i wklęsłe, i działa dość dobrze, jak liczba boków wzrasta.

                       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

To pytanie jest bardzo interesujące. Mam inny, wykonalny pomysł inny niż inne odpowiedzi tego postu.

Zapoznaj się z tym obrazem, aby uzyskać podstawowe zrozumienie idei:

Mój algorytm zakłada, że zgodnie z ruchem wskazówek zegara jest 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 Pythona, 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

Pracowałem nad tym, kiedy byłem badaczem Pod Michael Stonebraker - wiesz, profesor, który wymyślił Ingres, PostgreSQL , itd.

Zdaliśmy sobie sprawę, że najszybszym sposobem jest wykonanie obwiedni, ponieważ jest ona SUPER szybka. Jeśli jest poza obwiednią, to jest na zewnątrz. Inaczej wykonasz cięższą pracę...

Jeśli chcesz mieć świetny algorytm, poszukaj kodu źródłowego open source projektu PostgreSQL do pracy geo...

I chcemy podkreślić, że nigdy nie mieliśmy wglądu w praworęczność i leworęczność (również jako problem "wewnątrz" i "na zewnątrz"...

UPDATE

Link BKB dostarczył sporej liczby rozsądnych algorytmów. Pracowałem nad problemami Nauki o Ziemi i dlatego potrzebowałem rozwiązania, które działa w szerokości/długości geograficznej, i ma specyficzny problem ręki - czy obszar wewnątrz mniejszy obszar czy większy obszar? Odpowiedź jest taka ,że" kierunek " wierzchołków ma znaczenie - jest leworęczny lub praworęczny i w ten sposób możesz wskazać dowolny obszar jako "wewnątrz" dowolnego wielokąta. W związku z tym w mojej pracy wykorzystano rozwiązanie trzecie wymienione na tej stronie.

DODATKOWO w mojej pracy wykorzystywałam osobne funkcje do testów "on the line".

...Ponieważ ktoś zapytał: odkryliśmy, że testy obwiedni są najlepsze, gdy liczba wierzchołków przekracza pewną liczbę - zrób bardzo szybki test przed wykonaniem dłuższego testu, jeśli to konieczne... Obwiednia jest tworzona przez po prostu weź największe x, najmniejsze x, największe y i najmniejsze y i złóż je razem, aby utworzyć cztery punkty pudełka...

Kolejna wskazówka dla tych, którzy podążają za nami: wykonaliśmy wszystkie nasze bardziej wyrafinowane i" przyciemniające światło "obliczenia w przestrzeni siatki, wszystkie w dodatnich punktach na płaszczyźnie, a następnie ponownie rzutowaliśmy z powrotem na" rzeczywistą " długość/szerokość geograficzną, unikając w ten sposób możliwych błędów owinięcia wokół, gdy jedna linia przekroczyła 180 długości geograficznej i podczas obsługi regionów polarnych. Świetnie!

Naprawdę podoba mi się rozwiązanie opublikowane przez Nirg i edytowane przez bobobobo. Po prostu uczyniłem go przyjaznym dla javascript i trochę bardziej czytelnym dla mojego użytku:

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;
}

Swift version of the answer by 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
    }
}

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

Jak wiele optymalizacji, są one oparte na konkretnych, a nie ogólnych przypadkach i dają korzyści w oparciu o zamortyzowany czas, a nie pojedyncze użycie.

Pracując w tej dziedzinie, znalazłem Joeseph o 'Rourkes' Geometria obliczeniowa w C ' ISBN 0-521-44034-3 być wielką pomocą.

Trywialnym rozwiązaniem byłoby podzielenie wielokąta na trójkąty i sprawdzenie trójkątów zgodnie z wyjaśnieniem tutaj

Jeśli twój wielokąt jest wypukły może być jednak lepsze podejście. Spójrz na wielokąt jako zbiór nieskończonych linii. Każda linia dzieląca przestrzeń na dwie. dla każdego punktu łatwo powiedzieć, czy jego po jednej lub 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 stare, ale oto algorytm rzucania promieni zaimplementowany w Cocoa, na wypadek gdyby ktoś był zainteresowany. Nie jestem pewien, czy jest to najskuteczniejszy sposób na robienie 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;
}

Obj-C Wersja odpowiedzi nirg z przykładową metodą punktów testowych. Odpowiedź Nirg zadziałała dobrze dla mnie.

- (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");
    }
}

Ja też myślałem, że sumowanie 360 działa tylko dla wielokątów wypukłych, ale to nieprawda.

Ta strona ma ładny diagram pokazujący dokładnie to, i dobre wyjaśnienie na testowanie trafień: Gamasutra-w Nowy Rok: Wykrywanie kolizji

Wersja C# odpowiedzi nirg jest tutaj: po prostu udostępnię kod. To może zaoszczędzić komuś 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;
    }

}

. NET port:

    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;
        }
    }

Nie ma nic piękniejszego niż indukcyjna definicja problemu. Dla kompletności tutaj masz wersję w prologu, która może również wyjaśnić myśli za Ray casting :

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

Niektóre predykaty pomocnicze:

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 podane 2 punkty A i B (linia(A,B)) wynosi:

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

Jest ważne, aby kierunek obrotu linii był ustawiony na zegar-mądry dla granic i anty-zegar-mądry dla otworów. Sprawdzimy czy punkt (X, Y), czyli testowany punkt znajduje się po lewej stronie pół płaszczyzny naszej linii (to kwestia gustu, może być też prawej stronie, ale także kierunek linii granicznych musi być zmieniony w tego przypadku), polega to na projekcji promienia z punktu w prawo (lub w lewo) i potwierdź skrzyżowanie z linią. Zdecydowaliśmy się na projekt 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 jest po lewej (lub prawej) stronie tylko odcinek linii, a nie całą płaszczyznę, więc musimy ogranicz wyszukiwanie tylko do tego segmentu, ale jest to łatwe, ponieważ być wewnątrz segmentu tylko jeden punkt w linii może być wyższy niż Y w osi pionowej. Ponieważ jest to silniejsze ograniczenie to potrzeby być pierwszym, który sprawdzi, więc najpierw bierzemy tylko te linie spełnienie tego wymogu, a następnie sprawdź jego pozycję. Przez Jordan Twierdzenie o krzywej każdy promień rzutowany na wielokąt musi przecinać się w parzysta liczba linii. Więc skończyliśmy, rzucimy Promień do prawo i za każdym razem, gdy przecina linię, przełącza jej stan. Jednak w naszej realizacji jesteśmy goint do sprawdzenia długości worek rozwiązań spełniających podane ograniczenia i decydują o / align = "left" / dla każdej linii w wielokąt to trzeba zrobić.

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 VBA:

Uwaga: Pamiętaj, że jeśli twój wielokąt jest obszarem na mapie, że Szerokość / Długość geograficzna są wartościami Y/X w przeciwieństwie do X / Y (Szerokość = Y, Długość = X) ze względu na to, co rozumiem, są implikacje historyczne z przeszłości, kiedy Długość geograficzna nie była pomiarem.

CLASS MODULE: 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

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

/* 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ń do atan2f, ale może być interesująca dla programistów czytających ten wątek(w moich testach jego ~23x wolniej niż przy użyciu metody przecięcia linii).

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

1. Jeśli punkt znajduje się wewnątrz obszaru wielokąta. (można to osiągnąć za pomocą algorytmu odlewania promieni)
2. Jeśli punkt znajduje się na granicy wielokąta(może być osiągnięty za pomocą tego samego algorytmu, który jest używany do wykrywania punktów na polilinii(linii)).

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

1. promień pokrywa się z bokiem 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 z kątów wielokąta.

Sprawdzanie Określanie, Czy Punkt Znajduje Się Wewnątrz Złożonego Wielokąta. Artykuł zapewnia łatwy sposób na ich rozwiązanie, więc będzie nie wymaga specjalnego traktowania w powyższych przypadkach.

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

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

Nie wiem co dzieje się tak, jeśli twój wielokąt ma w sobie dziurę, ale wydaje mi się, że główną ideę można dostosować do tej sytuacji

Równie dobrze możesz opublikować pytanie w społeczności matematycznej. Założę się, że mają na to milion sposobów

Jeśli szukasz biblioteki java-script, istnieje rozszerzenie javascript google maps v3 dla klasy Polygon do wykrywania, czy punkt znajduje się w niej.

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

Google Extention Github

Używając qt (Qt 4.3+), można użyć funkcji QPolygon containsPoint

Odpowiedź zależy od tego, czy masz wielokąty proste czy złożone. Wielokąty proste nie mogą mieć przecięć segmentów linii. Więc mogą mieć dziury, ale linie nie mogą się krzyżować. Złożone regiony mogą mieć przecięcia linii, więc mogą mieć nakładające się regiony lub regiony, które stykają się ze sobą tylko przez jeden punkt.

Dla prostych wielokątów najlepszym algorytmem jest algorytm rzutowania promieni (liczby krzyżowej). Dla wielokątów złożonych algorytm ten nie wykrywa punktów które znajdują się wewnątrz nakładających się regionów. Więc dla złożonych wielokątów musisz użyć algorytmu nawijania liczb.

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 Scala rozwiązania przez nirg (zakłada, że wstępne sprawdzenie prostokąta jest wykonywane 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)
}

Wykonałem implementację Pythona C++ kod:

Wejścia

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

• Bounding_box_positions: kandydat wskazuje na filtr. (W mojej realizacji stworzonej z obwiedni.

  (Dane wejściowe są listami krotek w formacie: [(xcord, ycord), ...])

Zwraca

• wszystkie punkty znajdujące się 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ł pochodzi z tutaj

To działa tylko dla wypukłych kształtów, ale Minkowski Portal Refinement, i GJK są również świetne opcje do testowania, czy punkt jest w wielokątu. Użyj odejmowania Minkowskiego, aby odjąć punkt od wielokąta, a następnie uruchom te algorytmy, aby sprawdzić, czy wielokąt zawiera początek.

Ponadto, co ciekawe, można opisać swoje kształty nieco bardziej pośrednio za pomocą funkcji pomocniczych, które przyjmują Wektor kierunku jako wejście i wypluwają najdalszy punkt wzdłuż tego wektora. Pozwala to na możesz opisać dowolny wypukły kształt.. zakrzywione, zbudowane z wielokątów lub mieszane. Można również wykonywać operacje, aby połączyć wyniki prostych funkcji wsparcia, aby tworzyć bardziej złożone kształty.

Więcej informacji: http://xenocollide.snethen.com/mpr2d.html

Również, game programming gems 7 mówi o tym, jak to zrobić w 3d (: