Como posso determinar se um ponto 2D está dentro de um polígono?

Estou tentando criar um ponto 2D rápido dentro do algoritmo de polígono, para uso em testes de impacto (por exemplo Polygon.contains(p:Point)). Sugestões para técnicas eficazes seriam apreciadas.

Para gráficos, prefiro não preferir números inteiros. Muitos sistemas usam números inteiros para a pintura da interface do usuário (afinal, os pixels são ints), mas o macOS, por exemplo, usa float para tudo. O macOS conhece apenas pontos e um ponto pode ser traduzido para um pixel, mas, dependendo da resolução do monitor, pode ser traduzido para outra coisa. Nas telas de retina, meio ponto (0,5 / 0,5) é pixel. Ainda assim, nunca notei que as interfaces do macOS são significativamente mais lentas que as outras. Afinal, todas as APIs 3D (OpenGL ou Direct3D) também funcionam com carros alegóricos e as bibliotecas gráficas modernas costumam tirar proveito da aceleração da GPU.

Agora você disse que a velocidade é a sua principal preocupação, ok, vamos em velocidade. Antes de executar qualquer algoritmo sofisticado, primeiro faça um teste simples. Crie uma caixa delimitadora alinhada ao eixo em torno do seu polígono. Isso é muito fácil, rápido e já pode proteger muitos cálculos. Como isso funciona? Itere todos os pontos do polígono e encontre os valores mínimo / máximo de X e Y.

Por exemplo, você tem os pontos (9/1), (4/3), (2/7), (8/2), (3/6). Isso significa que Xmin é 2, Xmax é 9, Ymin é 1 e Ymax é 7. Um ponto fora do retângulo com as duas arestas (2/1) e (9/7) não pode estar dentro do polígono.

// 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!
}

Este é o primeiro teste a ser executado em qualquer ponto. Como você pode ver, este teste é ultra rápido, mas também é muito grosseiro. Para lidar com pontos que estão dentro do retângulo delimitador, precisamos de um algoritmo mais sofisticado. Existem algumas maneiras de como isso pode ser calculado. O método que funciona também depende do fato de o polígono poder ter furos ou sempre ser sólido. Aqui estão exemplos de sólidos (um convexo, um côncavo):

E aqui está um com um buraco:

O verde tem um buraco no meio!

O algoritmo mais fácil, que pode lidar com todos os três casos acima e ainda é bastante rápido, é chamado de conversão de raios . A idéia do algoritmo é bastante simples: desenhe um raio virtual de qualquer lugar fora do polígono até o seu ponto e conte quantas vezes ele atinge um lado do polígono. Se o número de acertos for par, está fora do polígono, se for ímpar, está dentro.

O algoritmo do número do enrolamento seria uma alternativa, é mais preciso para os pontos estarem muito próximos de uma linha poligonal, mas também é muito mais lento. A conversão de raios pode falhar em pontos muito próximos ao lado de um polígono devido a problemas limitados de precisão e arredondamento de pontos flutuantes, mas, na realidade, isso dificilmente é um problema, como se um ponto estivesse tão próximo de um lado, geralmente não é visualmente possível para um visualizador para reconhecer se ele já está dentro ou fora.

Você ainda tem a caixa delimitadora acima, lembra? Basta escolher um ponto fora da caixa delimitadora e usá-lo como ponto de partida para o seu raio. Por exemplo, o ponto (Xmin - e/p.y)está fora do polígono, com certeza.

Mas o que é isso e? Bem, e(na verdade epsilon), dá à caixa delimitadora um pouco de preenchimento . Como eu disse, o traçado de raios falha se começarmos muito perto de uma linha poligonal. Como a caixa delimitadora pode ser igual ao polígono (se o polígono é um retângulo alinhado ao eixo, a caixa delimitadora é igual ao próprio polígono!), Precisamos de algum preenchimento para tornar isso seguro, é tudo. Quão grande você deve escolher e? Não muito grande. Depende da escala do sistema de coordenadas que você usa para desenhar. Se a largura da etapa do pixel for 1.0, basta escolher 1,0 (ainda assim, 0,1 também funcionaria)

Agora que temos o raio com suas coordenadas de início e fim, o problema muda de " é o ponto dentro do polígono " para "com que frequência o raio cruza um lado do polígono ". Portanto, não podemos apenas trabalhar com os pontos poligonais como antes, agora precisamos dos lados reais. Um lado é sempre definido por dois pontos.

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

Você precisa testar o raio contra todos os lados. Considere o raio como um vetor e cada lado como um vetor. O raio tem que atingir cada lado exatamente uma vez ou nunca. Não pode atingir o mesmo lado duas vezes. Duas linhas no espaço 2D sempre se cruzam exatamente uma vez, a menos que sejam paralelas; nesse caso, nunca se cruzam. No entanto, como os vetores têm um comprimento limitado, dois vetores podem não ser paralelos e ainda nunca se cruzarem, porque são muito curtos para se encontrarem.

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

Até aqui tudo bem, mas como você testa se dois vetores se cruzam? Aqui está um código C (não testado), que deve fazer o truque:

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

Os valores de entrada são os dois pontos finais do vetor 1 ( v1x1/v1y1e v1x2/v1y2) e do vetor 2 ( v2x1/v2y1e v2x2/v2y2). Então você tem 2 vetores, 4 pontos, 8 coordenadas. YESe NOsão claros. YESaumenta interseções, NOnão faz nada.

E o COLLINEAR? Isso significa que ambos os vetores estão na mesma linha infinita, dependendo da posição e do comprimento, não se cruzam nem se interceptam em um número infinito de pontos. Não tenho certeza absoluta de como lidar com este caso, não o contaria como interseção de qualquer maneira. Bem, este caso é bastante raro na prática, devido a erros de arredondamento de ponto flutuante; código melhor provavelmente não testaria, == 0.0fmas sim algo como < epsilon, onde epsilon é um número bastante pequeno.

Se você precisar testar um número maior de pontos, certamente poderá acelerar um pouco a coisa toda, mantendo na memória as formas padrão das equações lineares dos lados dos polígonos, para que você não precise recalculá-las todas as vezes. Isso economizará duas multiplicações de ponto flutuante e três subtrações de ponto flutuante em cada teste, em troca de armazenar três valores de ponto flutuante por lado do polígono na memória. É uma troca típica de memória versus tempo de computação.

Por último, mas não menos importante: se você pode usar o hardware 3D para resolver o problema, existe uma alternativa interessante. Apenas deixe a GPU fazer todo o trabalho para você. Crie uma superfície de pintura que esteja fora da tela. Encha-o completamente com a cor preta. Agora deixe o OpenGL ou o Direct3D pintar seu polígono (ou mesmo todos os seus polígonos, se você quiser apenas testar se o ponto está em algum deles, mas você não se importa com qual deles) e preencher o (s) polígono (s) com um diferente cor, por exemplo, branco. Para verificar se um ponto está dentro do polígono, obtenha a cor desse ponto na superfície de desenho. Esta é apenas uma busca de memória O (1).

É claro que esse método só pode ser usado se a superfície do desenho não precisar ser grande. Se não puder caber na memória da GPU, esse método é mais lento do que na CPU. Se precisar ser enorme e sua GPU suportar shaders modernos, você ainda poderá usá-la implementando a transmissão de raios mostrada acima como shader de GPU, o que é absolutamente possível. Para um número maior de polígonos ou um grande número de pontos a serem testados, isso será recompensado, considere que algumas GPUs poderão testar 64 a 256 pontos em paralelo. Observe, no entanto, que transferir dados da CPU para a GPU e vice-versa é sempre caro, portanto, apenas para testar alguns pontos em alguns polígonos simples, em que os pontos ou os polígonos são dinâmicos e mudam com frequência, uma abordagem da GPU raramente paga fora.

Eu acho que o seguinte trecho de código é a melhor solução (retirada daqui ):

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

Argumentos

• nvert : Número de vértices no polígono. A repetição do primeiro vértice no final foi discutida no artigo mencionado acima.
• vertx, verty : matrizes contendo as coordenadas x e y dos vértices do polígono.
• testx, irritável : coordenadas X e Y do ponto de teste.

É curto e eficiente e funciona tanto para polígonos convexos quanto côncavos. Como sugerido anteriormente, você deve verificar primeiro o retângulo delimitador e tratar os orifícios dos polígonos separadamente.

A ideia por trás disso é bastante simples. O autor descreve da seguinte maneira:

Corro um raio semi-infinito horizontalmente (aumentando x, y fixo) a partir do ponto de teste e conto quantas arestas ele cruza. Em cada cruzamento, o raio alterna entre interior e exterior. Isso é chamado de teorema da curva de Jordan.

A variável c está alternando de 0 para 1 e 1 para 0 cada vez que o raio horizontal cruza qualquer aresta. Então, basicamente, é acompanhar se o número de arestas cruzadas é par ou ímpar. 0 significa par e 1 significa ímpar.

Aqui está uma versão em C # da resposta dada por nirg , que vem deste professor de RPI . Observe que o uso do código dessa fonte RPI requer atribuição.

Uma seleção de caixa delimitadora foi adicionada na parte superior. No entanto, como James Brown aponta, o código principal é quase tão rápido quanto a caixa delimitadora, então a verificação da caixa delimitadora pode realmente retardar a operação geral, no caso de a maioria dos pontos que você está verificando estarem dentro da caixa delimitadora . Portanto, você pode deixar a caixa delimitadora ou uma alternativa seria pré-calcular as caixas delimitadoras de seus polígonos se elas não mudarem de forma com muita frequência.

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

Aqui está uma variante JavaScript da resposta de M. Katz com base na abordagem de 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;
}

Calcule a soma orientada dos ângulos entre o ponto pe cada um dos ápices poligonais. Se o ângulo total orientado for 360 graus, o ponto estará dentro. Se o total for 0, o ponto está fora.

Eu gosto mais desse método porque é mais robusto e menos dependente da precisão numérica.

Os métodos que calculam a uniformidade do número de interseções são limitados porque você pode "atingir" um ápice durante o cálculo do número de interseções.

EDIT: By the Way, este método funciona com polígonos côncavos e convexos.

EDIT: Recentemente, encontrei um artigo completo da Wikipedia sobre o assunto.

Esta questão é tão interessante. Tenho outra ideia viável diferente de outras respostas a este post. A idéia é usar a soma dos ângulos para decidir se o alvo está dentro ou fora. Mais conhecido como número do enrolamento .

Seja x o ponto alvo. Seja o array [0, 1, .... n] todos os pontos da área. Conecte o ponto de destino a cada ponto de borda com uma linha. Se o ponto de destino estiver dentro desta área. A soma de todos os ângulos será 360 graus. Caso contrário, os ângulos serão inferiores a 360.

Consulte esta imagem para obter um entendimento básico da ideia:

Meu algoritmo assume que o sentido horário é a direção positiva. Aqui está uma entrada em potencial:

[[-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]]

A seguir está o código python que implementa a ideia:

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

O artigo de Eric Haines citado por bobobobo é realmente excelente. Particularmente interessantes são as tabelas que comparam o desempenho dos algoritmos; o método da soma dos ângulos é realmente ruim comparado aos outros. Também interessante é que otimizações como o uso de uma grade de pesquisa para subdividir ainda mais o polígono em setores "dentro" e "fora" podem tornar o teste incrivelmente rápido, mesmo em polígonos com> 1000 lados.

Enfim, é cedo, mas meu voto vai para o método "travessias", que é basicamente o que Mecki descreve, eu acho. No entanto, achei mais sucintamente descrito e codificado por David Bourke . Eu amo que não há trigonometria real necessária, e funciona para convexo e côncavo, e funciona razoavelmente bem à medida que o número de lados aumenta.

A propósito, aqui está uma das tabelas de desempenho do artigo de Eric Haines, que interessa, testando em polígonos aleatórios.

                       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

Versão rápida da resposta por 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
    }
}

Realmente gosto da solução publicada por Nirg e editada por bobobobo. Acabei de torná-lo compatível com javascript e um pouco mais legível para o meu uso:

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

Eu trabalhei nisso quando era pesquisador do Michael Stonebraker - você sabe, o professor que criou Ingres , PostgreSQL , etc.

Percebemos que a maneira mais rápida era primeiro fazer uma caixa delimitadora porque é SUPER rápido. Se estiver fora da caixa delimitadora, está fora. Caso contrário, você faz o trabalho mais difícil ...

Se você deseja um ótimo algoritmo, consulte o código-fonte PostgreSQL do projeto de código aberto para o trabalho geográfico ...

Quero ressaltar que nunca tivemos nenhuma visão sobre destros versus canhotos (também expressável como um problema "interno" vs "externo" ...

ATUALIZAR

O link do BKB forneceu um bom número de algoritmos razoáveis. Eu estava trabalhando nos problemas das ciências da terra e, portanto, precisava de uma solução que funcionasse em latitude / longitude, e ela tem o problema peculiar de destreza - é a área dentro da área menor ou maior? A resposta é que a "direção" dos vértices é importante - é canhoto ou destro e, dessa maneira, você pode indicar uma das áreas como "dentro" de qualquer polígono. Como tal, meu trabalho usou a solução três enumerada nessa página.

Além disso, meu trabalho usou funções separadas para testes "na linha".

... Como alguém perguntou: descobrimos que os testes de caixa delimitadora eram melhores quando o número de vértices ultrapassava algum número - faça um teste muito rápido antes de fazer o teste mais longo, se necessário ... Uma caixa delimitadora é criada simplesmente fazendo o x maior, menor x, maior y e menor y e juntando-os para formar quatro pontos de uma caixa ...

Outra dica para aqueles que se seguem: fizemos toda a nossa computação mais sofisticada e com "escurecimento da luz" em um espaço de grade, todos em pontos positivos em um plano e, em seguida, re-projetamos de volta em longitude / latitude "real", evitando possíveis erros de envolvendo quando cruzamos a linha 180 de longitude e ao lidar com regiões polares. Trabalhou muito bem!

A resposta de David Segond é praticamente a resposta geral padrão, e a de Richard T é a otimização mais comum, embora existam outras. Outras otimizações fortes são baseadas em soluções menos gerais. Por exemplo, se você estiver indo para verificar o mesmo polígono com muitos pontos, a triangulação do polígono pode acelerar bastante as coisas, pois existem vários algoritmos de pesquisa TIN muito rápidos. Outra é que, se o polígono e os pontos estiverem em um plano limitado em baixa resolução, digamos uma tela, você poderá pintar o polígono em um buffer de tela mapeado na memória em uma determinada cor e verificar a cor de um determinado pixel para ver se está nos polígonos.

Como muitas otimizações, elas são baseadas em casos específicos, e não gerais, e geram benefícios com base no tempo amortizado, e não no uso único.

Trabalhando neste campo, achei a Geometria de computação de Joeseph O'Rourkes no C 'ISBN 0-521-44034-3 uma grande ajuda.

A solução trivial seria dividir o polígono em triângulos e testar os triângulos como explicado aqui

Se o seu polígono for CONVEX , pode haver uma abordagem melhor. Veja o polígono como uma coleção de linhas infinitas. Cada linha dividindo o espaço em dois. para cada ponto, é fácil dizer se está de um lado ou do outro lado da linha. Se um ponto estiver no mesmo lado de todas as linhas, ele estará dentro do polígono.

Sei que isso é antigo, mas aqui está um algoritmo de conversão de raios implementado no cacau, caso alguém esteja interessado. Não tenho certeza se é a maneira mais eficiente de fazer as coisas, mas pode ajudar alguém.

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

Versão Obj-C da resposta de nirg com método de amostra para pontos de teste. A resposta de Nirg funcionou bem para mim.

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

Não há nada mais bonito do que uma definição indutiva de um problema. Por uma questão de completude, aqui você tem uma versão em prólogo que também pode esclarecer os pensamentos por trás da fundição de raios :

Com base na simulação do algoritmo de simplicidade em http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Alguns predicados auxiliares:

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]).

A equação de uma reta dada 2 pontos A e B (Linha (A, B)) é:

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

É importante que a direção de rotação da linha seja ajustada no sentido horário para limites e anti-horário para orifícios. Vamos verificar se o ponto (X, Y), ou seja, o ponto testado está no semiplano esquerdo da nossa linha (é uma questão de gosto, também pode ser o lado direito, mas também a direção dos limites Nesse caso, as linhas devem ser alteradas), para projetar o raio do ponto para a direita (ou esquerda) e reconhecer a interseção com a linha. Optamos por projetar o raio na direção horizontal (novamente, é uma questão de gosto, também pode ser feito na vertical, com restrições semelhantes), por isso temos:

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

Agora precisamos saber se o ponto está no lado esquerdo (ou direito) apenas do segmento de linha, não o plano inteiro, portanto, precisamos restringir a pesquisa apenas a esse segmento, mas isso é fácil, pois está dentro do segmento apenas um ponto na linha pode ser maior que Y no eixo vertical. Como essa é uma restrição mais forte, ela precisa ser a primeira a verificar, portanto, pegamos primeiro apenas as linhas que atendem a esse requisito e, em seguida, verificamos sua possibilidade. Pelo teorema da Curva de Jordan, qualquer raio projetado em um polígono deve se cruzar em um número par de linhas. Assim que terminamos, lançaremos o raio para a direita e sempre que ele cruzar uma linha, alternará seu estado. No entanto, em nossa implementação, vamos verificar o comprimento do pacote de soluções que atendem às restrições especificadas e decidir a adesão a ele. para cada linha no polígono, isso deve ser feito.

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).

A versão em C # da resposta de nirg está aqui: apenas compartilharei o código. Isso pode economizar tempo para alguém.

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

Versão 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;
    }

}

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

VERSÃO VBA:

Nota: Lembre-se de que, se o seu polígono for uma área dentro de um mapa, Latitude / Longitude são valores Y / X, em oposição a X / Y (Latitude = Y, Longitude = X), devido ao que eu entendo, há implicações históricas desde quando Longitude não era uma medida.

MÓDULO DE CLASSE: 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

MÓDULO:

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

Eu fiz uma implementação Python de de nirg c ++ código :

Entradas

• bounding_points: nós que compõem o polígono.

• bounding_box_positions: pontos candidatos a serem filtrados. (Na minha implementação criada a partir da caixa delimitadora.

  (As entradas são listas de tuplas no formato: [(xcord, ycord), ...])

Devoluções

• Todos os pontos que estão dentro do polígono.

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

Novamente, a idéia é tirada daqui

Surpreso, ninguém falou disso antes, mas para os pragmáticos que precisam de um banco de dados: o MongoDB tem excelente suporte para consultas geográficas, incluindo esta.

O que você está procurando é:

db.neighborhoods.findOne ({geometry: {$ geoIntersects: {$ geometry: {type: "Point", coordenadas: ["longitude", "latitude"]}}}}))

Neighborhoodsé a coleção que armazena um ou mais polígonos no formato GeoJson padrão. Se a consulta retornar nula, ela não será interceptada, caso contrário, será.

Muito bem documentado aqui: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

O desempenho de mais de 6.000 pontos classificados em uma grade irregular de 330 polígonos foi inferior a um minuto sem nenhuma otimização e incluindo o tempo para atualizar documentos com o respectivo polígono.

Aqui está um ponto no teste de polígono em C que não está usando a projeção de raios. E pode funcionar para áreas sobrepostas (auto-interseções), veja o use_holesargumento.

/* 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];
}

Nota: esse é um dos métodos menos ótimos, pois inclui muitas chamadas para atan2f, mas pode ser interessante para os desenvolvedores que estão lendo este segmento (nos meus testes é ~ 23x mais lento que o método de interseção de linhas).

Para detectar hit no Polygon, precisamos testar duas coisas:

1. Se Point estiver dentro da área de polígono. (pode ser realizado pelo algoritmo de fundição por raio)
2. Se Point estiver na borda do polígono (pode ser realizado pelo mesmo algoritmo usado para a detecção de pontos na polilinha (linha)).

Para lidar com os seguintes casos especiais no algoritmo de conversão de raios :

1. O raio se sobrepõe a um dos lados do polígono.
2. O ponto está dentro do polígono e o raio passa por um vértice do polígono.
3. O ponto está fora do polígono e o raio toca apenas um dos ângulos do polígono.

Marque Determinando se um ponto está dentro de um polígono complexo . O artigo fornece uma maneira fácil de resolvê-los, para que não haja tratamento especial necessário para os casos acima.

Você pode fazer isso verificando se a área formada conectando o ponto desejado aos vértices do polígono corresponde à área do próprio polígono.

Ou você pode verificar se a soma dos ângulos internos do seu ponto para cada par de dois vértices poligonais consecutivos até o seu ponto de verificação é de 360, mas sinto que a primeira opção é mais rápida porque não envolve divisões nem cálculos inverso de funções trigonométricas.

Não sei o que acontece se o seu polígono tiver um furo dentro dele, mas me parece que a idéia principal pode ser adaptada a essa situação

Você também pode postar a pergunta em uma comunidade de matemática. Aposto que eles têm um milhão de maneiras de fazer isso

Se você está procurando uma biblioteca de scripts java, há uma extensão javascript do google maps v3 para a classe Polygon para detectar se um ponto reside ou não nela.

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

Github de extensão do Google

Ao usar qt(Qt 4.3+), pode-se usar a função QPolygon containsPoint

A resposta depende se você possui polígonos simples ou complexos. Polígonos simples não devem ter nenhuma interseção de segmento de linha. Eles podem ter os buracos, mas as linhas não podem se cruzar. Regiões complexas podem ter interseções de linha - para que possam ter regiões sobrepostas ou regiões que se tocam apenas por um único ponto.

Para polígonos simples, o melhor algoritmo é o algoritmo de vazamento de raios (número de cruzamento). Para polígonos complexos, esse algoritmo não detecta pontos que estão dentro das regiões sobrepostas. Portanto, para polígonos complexos, é necessário usar o algoritmo de número do enrolamento.

Aqui está um excelente artigo com a implementação em C de ambos os algoritmos. Eu tentei e eles funcionam bem.

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

Versão Scala da solução por nirg (assume que a pré-verificação do retângulo delimitador é feita separadamente):

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

Aqui está a versão golang da resposta @nirg (inspirada no código C # de @@ 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
}