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

497

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.

Scott Evernden
fonte
Você esqueceu de nos contar sobre suas percepções sobre a questão da mão direita ou esquerda - que também pode ser interpretada como "dentro" versus "fora" - RT
Richard T
13
Sim, agora percebo que a pergunta deixa muitos detalhes não especificados, mas neste momento estou meio interessado em ver a variedade de respostas.
23380 Scott Evernden
4
Um polígono de 90 lados é chamado de eneacontágono e um polígono de 10.000 lados é chamado de miríade.
"Mais elegante" está fora do objetivo, pois tive problemas em encontrar um algoritmo de "trabalho de todo". Devo descobrir sozinho: stackoverflow.com/questions/14818567/...
davidkonrad

Respostas:

732

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

Polígono sem furo

E aqui está um com um buraco:

Polígono com furo

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.

Demonstrando como o raio corta um polígono

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.

Mecki
fonte
26
+1 resposta fantástica. Ao usar o hardware para fazê-lo, li em outros lugares que pode ser lento porque você precisa recuperar os dados da placa gráfica. Mas eu gosto muito do princípio de retirar a carga da CPU. Alguém tem boas referências de como isso pode ser feito no OpenGL?
Gavin
3
+1 porque é muito simples! O principal problema é que, se o seu polígono e ponto de teste estiverem alinhados em uma grade (não incomum), você precisará lidar com interseções "duplicadas", por exemplo, diretamente através de um ponto poligonal! (produzindo uma paridade de dois em vez de um). Entra nesta área estranha: stackoverflow.com/questions/2255842/… . A computação gráfica está cheia desses casos especiais: simples em teoria, cabeludos na prática.
precisa
7
@RMorrisey: Por que você acha isso? Não vejo como isso falharia em um polígono côncavo. O raio pode sair e entrar novamente no polígono várias vezes quando o polígono é côncavo, mas no final, o contador de ocorrências será ímpar se o ponto estiver dentro e mesmo se estiver fora, também para polígonos côncavos.
Mecki 31/05
6
O 'Algoritmo de número de enrolamento rápido', descrito em softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm, funciona muito rápido ...
SP
10
Seu uso de / para separar as coordenadas x e y é confuso; ele lê como x dividido por y. É muito mais claro usar x, y (ou seja, x vírgula y) No geral, uma resposta útil.
Ash
583

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.

nirg
fonte
5
Questão. Quais são exatamente as variáveis ​​que eu passo? O que eles representam?
tekknolagi 02/02/12
9
@Mick Verifica isso verty[i]e verty[j]os dois lados testy, para que nunca sejam iguais.
Peter Wood
4
Este código não é robusto e eu não recomendaria usá-lo. Aqui está um link que dá uma análise detalhada: www-ma2.upc.es/geoc/Schirra-pointPolygon.pdf
Mikola
13
Essa abordagem realmente possui limitações (casos extremos): a verificação do ponto (15,20) no polígono [(10,10), (10,20), (20,20), (20,10)] retornará falso em vez de verdadeiro. Mesmo com (10,20) ou (20,15). Para todos os outros casos, o algoritmo funciona perfeitamente bem e os falsos negativos nos casos extremos são válidos para a minha aplicação.
21814 Alexander Pacha
10
@ Alexander, isso é de fato por design: manipulando os limites esquerdo e inferior no sentido oposto aos limites superior e direito, se dois polígonos distintos compartilharem uma aresta, qualquer ponto ao longo dessa aresta estará localizado em um e apenas um polígono. ..uma propriedade útil.
wardw
69

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;
}
M Katz
fonte
5
Funciona muito bem, obrigado, eu converti para JavaScript. stackoverflow.com/questions/217578/…
Philipp Lenssen
2
Isso é> 1000x mais rápido do que usar GraphicsPath.IsVisible !! A verificação da caixa delimitadora torna a função cerca de 70% mais lenta.
James Brown,
Não apenas GraphicsPath.IsVisible () é muito mais lento, mas também não funciona bem com polígonos muito pequenos com lado na faixa 0,01f
Patrick da equipe NDepend
50

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;
}
Philipp Lenssen
fonte
32

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.

David Segonds
fonte
1
Não, isso não é verdade. Isso funciona independentemente da convexidade do polígono.
David Segonds 20/10/08
2
@ DarenW: Apenas um acos por vértice; por outro lado, esse algoritmo deve ser o mais fácil de implementar no SIMD, pois não possui nenhuma ramificação.
Jasper Bekkers
1
@emilio, se o ponto estiver longe do triângulo, não vejo como o ângulo formado pelo ponto e dois vértices do triângulo será de 90 graus.
11119 David Segonds
2
Primeiro use a caixa delimitadora para resolver os casos "point is distante". Para trig, você pode usar tabelas pré-geradas.
JOM
3
Esta é a solução ideal, pois é O (n) e requer computação mínima. Funciona para todos os polígonos. Eu pesquisei essa solução há 30 anos no meu primeiro trabalho na IBM. Eles assinaram e ainda o usam hoje em suas tecnologias GIS.
Dominic Cerisano
24

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: insira a descrição da imagem aqui

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
Junbang Huang
fonte
21

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
Gavin
fonte
11

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
    }
}
bzz
fonte
Isso tem um potencial de dividir por zero no cálculo b. É necessário calcular apenas "b" e a próxima linha com "c" se o cálculo para "a" mostrar que existe a possibilidade de interseção. Não há possibilidade de ambos os pontos estarem acima ou ambos os pontos abaixo - o que é descrito pelo cálculo para "a".
anorskdev 23/03
11

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;
}
Dave Seidman
fonte
8

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!

Richard T
fonte
E se eu não tiver a caixa delimitadora? :)
Scott Evernden
8
Você pode criar facilmente uma caixa delimitadora - são apenas os quatro pontos que usam o maior e o menor x e o maior e o menor y. Mais cedo.
Richard T
"... evitando possíveis erros de contornar quando se cruzava a linha 180 de longitude e ao lidar com regiões polares." você pode descrever isso com mais detalhes? Eu acho que consigo descobrir como mover tudo 'para cima' para evitar o meu polígono cruzando a longitude 0, mas não estou claro sobre como lidar com polígonos que contêm um dos pólos ...
tiritea
6

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.

SmacL
fonte
4

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.

shoosh
fonte
muito rápido e pode ser aplicado a formas mais gerais. por volta de 1990, tínhamos "curvigons" onde alguns lados eram arcos circulares. Ao analisar esses lados em cunhas circulares e um par de triângulos unindo a cunha à origem (centróide poligonal), o teste de acerto foi rápido e fácil.
1111 DarenW
1
tem alguma referência nesses curvigons?
21239 shoosh
Eu adoraria um árbitro para os curvigons também.
Joel em Gö
Você perdeu uma solução importante para o caso de polígonos convexos: comparando o ponto com uma diagonal, você pode reduzir pela metade o número de vértices. E repetindo este processo, você reduz a um único triângulo no Log (N) operações, em vez de N.
Yves Daoust
4

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;
}
diatrevolo
fonte
5
Observe que se você realmente estiver fazendo isso no Cocoa, poderá usar o método [NSBezierPath containsPoint:].
Thomasw
4

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

polígono de amostra

Jon
fonte
2
Claro, no Objective-C, CGPathContainsPoint()é seu amigo.
Zev Eisenberg
@ZevEisenberg, mas isso seria fácil demais! Obrigado pela observação. Vou desenterrar esse projeto em algum momento para ver por que usei uma solução personalizada. Provavelmente eu não sabia sobreCGPathContainsPoint()
Jon
4

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).
jdavid_1385
fonte
3

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;
        }
Uğur Gümüşhan
fonte
isso funciona na maioria dos casos, mas está errado e não funciona corretamente sempre! utilizar a solução de M Katz que é correcto
Lukas Hanacek
3

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

}
YongJiang Zhang
fonte
2

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;
        }
    }
Aladar
fonte
2

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
Colin Stadig
fonte
2

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

Noresourses
fonte
2

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.

Santiago M. Quintero
fonte
1

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

ideasman42
fonte
0

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)).
VJ
fonte
0

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.

Justin
fonte
0

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

user5193682
fonte
0

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

shana
fonte
0

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

Timmy_A
fonte
0

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)
}
Michael-7
fonte
0

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
}
SamTech
fonte