Classificar pontos no sentido horário?

156

Dada uma matriz de x, y pontos, como classifico os pontos dessa matriz na ordem dos ponteiros do relógio (em torno de seu ponto central médio geral)? Meu objetivo é passar os pontos para uma função de criação de linhas para terminar com algo parecendo "sólido", o mais convexo possível, sem que as linhas se cruzem.

Pelo que vale, estou usando Lua, mas qualquer pseudocódigo seria apreciado.

Atualização: Para referência, este é o código Lua baseado na excelente resposta de Ciamej (ignore meu prefixo "app"):

function appSortPointsClockwise(points)
    local centerPoint = appGetCenterPointOfPoints(points)
    app.pointsCenterPoint = centerPoint
    table.sort(points, appGetIsLess)
    return points
end

function appGetIsLess(a, b)
    local center = app.pointsCenterPoint

    if a.x >= 0 and b.x < 0 then return true
    elseif a.x == 0 and b.x == 0 then return a.y > b.y
    end

    local det = (a.x - center.x) * (b.y - center.y) - (b.x - center.x) * (a.y - center.y)
    if det < 0 then return true
    elseif det > 0 then return false
    end

    local d1 = (a.x - center.x) * (a.x - center.x) + (a.y - center.y) * (a.y - center.y)
    local d2 = (b.x - center.x) * (b.x - center.x) + (b.y - center.y) * (b.y - center.y)
    return d1 > d2
end

function appGetCenterPointOfPoints(points)
    local pointsSum = {x = 0, y = 0}
    for i = 1, #points do pointsSum.x = pointsSum.x + points[i].x; pointsSum.y = pointsSum.y + points[i].y end
    return {x = pointsSum.x / #points, y = pointsSum.y / #points}
end

Philipp Lenssen
fonte
1
Pense em calcular o ângulo da linha radial através desse ponto. Em seguida, classifique por ângulo.
Presidente James K. Polk
Caso você não saiba, lua possui uma função interna ipairs(tbl)que itera sobre os índices e valores de tbl de 1 a #tbl. Assim, para o cálculo soma, você pode fazer isso, o que a maioria das pessoas acham que parece mais limpo:for _, p in ipairs(points) do pointsSum.x = pointsSum.x + p.x; pointsSum.y = pointsSum.y + p.y end
Ponkadoodle
2
@ Wallacoloo Isso é altamente discutível. Além disso, na baunilha, Lua ipairsé significativamente mais lenta que o numérico para loop.
Alexander Gladysh
Eu tive que fazer algumas pequenas alterações para que funcionasse no meu caso (apenas comparando dois pontos em relação a um centro). gist.github.com/personalnadir/6624172 Todas essas comparações com 0 no código parecem assumir que os pontos estão distribuídos em torno da origem, em oposição a um ponto arbitrário. Eu também acho que a primeira condição classificará os pontos abaixo do ponto central incorretamente. Obrigado pelo código, porém, tem sido realmente útil!
personalnadir

Respostas:

192

Primeiro, calcule o ponto central. Em seguida, classifique os pontos usando o algoritmo de classificação que desejar, mas use uma rotina de comparação especial para determinar se um ponto é menor que o outro.

Você pode verificar se um ponto (a) está à esquerda ou à direita do outro (b) em relação ao centro por este cálculo simples:

det = (a.x - center.x) * (b.y - center.y) - (b.x - center.x) * (a.y - center.y)

se o resultado for zero, eles estarão na mesma linha do centro, se for positivo ou negativo, estará de um lado ou de outro, de modo que um ponto precederá o outro. Com ele, é possível construir uma relação menor que para comparar os pontos e determinar a ordem em que eles devem aparecer na matriz classificada. Mas você precisa definir onde está o início dessa ordem, quero dizer qual ângulo será o inicial (por exemplo, a metade positiva do eixo x).

O código para a função de comparação pode ser assim:

bool less(point a, point b)
{
    if (a.x - center.x >= 0 && b.x - center.x < 0)
        return true;
    if (a.x - center.x < 0 && b.x - center.x >= 0)
        return false;
    if (a.x - center.x == 0 && b.x - center.x == 0) {
        if (a.y - center.y >= 0 || b.y - center.y >= 0)
            return a.y > b.y;
        return b.y > a.y;
    }

    // compute the cross product of vectors (center -> a) x (center -> b)
    int det = (a.x - center.x) * (b.y - center.y) - (b.x - center.x) * (a.y - center.y);
    if (det < 0)
        return true;
    if (det > 0)
        return false;

    // points a and b are on the same line from the center
    // check which point is closer to the center
    int d1 = (a.x - center.x) * (a.x - center.x) + (a.y - center.y) * (a.y - center.y);
    int d2 = (b.x - center.x) * (b.x - center.x) + (b.y - center.y) * (b.y - center.y);
    return d1 > d2;
}

Isso ordenará os pontos no sentido horário a partir das 12 horas. Os pontos na mesma "hora" serão ordenados a partir dos pontos mais afastados do centro.

Se você estiver usando tipos inteiros (que não estão realmente presentes em Lua), será necessário garantir que as variáveis ​​det, d1 e d2 sejam de um tipo capaz de manter o resultado dos cálculos realizados.

Se você deseja alcançar algo com uma aparência sólida, o mais convexa possível, acho que está procurando um casco convexo . Você pode calculá-lo usando o Graham Scan . Nesse algoritmo, você também precisa classificar os pontos no sentido horário (ou anti-horário), começando em um ponto de articulação especial. Em seguida, você repete etapas simples do loop cada vez que verifica se você vira para a esquerda ou para a direita adicionando novos pontos ao casco convexo. Essa verificação é baseada em um produto cruzado, como na função de comparação acima.

Editar:

Adicionada mais uma instrução if if (a.y - center.y >= 0 || b.y - center.y >=0)para garantir que os pontos que têm x = 0 e y negativo sejam classificados começando pelos pontos mais distantes do centro. Se você não se importa com a ordem dos pontos na mesma 'hora', pode omitir esta declaração if e sempre retornar a.y > b.y.

Corrigida a primeira instrução if com adição -center.xe -center.y.

Adicionada a segunda instrução if (a.x - center.x < 0 && b.x - center.x >= 0). Era óbvio que estava faltando. As instruções if podem ser reorganizadas agora porque algumas verificações são redundantes. Por exemplo, se a primeira condição na primeira instrução se for falsa, a primeira condição da segunda se deve ser verdadeira. Decidi, contudo, deixar o código como é por uma questão de simplicidade. É bem possível que o compilador otimize o código e produza o mesmo resultado de qualquer maneira.

ciamej
fonte
25
+1: Não atan(), nenhuma raiz quadrada e até mesmo nenhuma divisão. Este é um bom exemplo de pensamento em computação gráfica. Escolha todos os casos fáceis o mais rápido possível, e mesmo nos casos difíceis, calcule o mínimo possível para saber a resposta necessária.
RBerteig
Mas isso requer comparar todos os pontos com todos os outros. Existe um método simples de inserir novos pontos?
Iterator
2
se o conjunto de pontos é conhecido a priori, são necessárias apenas comparações O (n * log n). Se você quiser adicionar pontos nesse meio tempo, precisará mantê-los em um conjunto classificado, como uma árvore de pesquisa binária equilibrada. Nesse caso, adicionar um novo ponto requer comparações O (log n) e é exatamente o mesmo para a solução que envolve coordenadas polares.
ciamej
2
Está faltando o caso: if (ax - center.x <0 && bx - center.x> = 0) retorna false;
Tom Martin
2
Olá. É bem antigo, mas: "Isso ordenará os pontos no sentido horário a partir das 12 horas". Por que 12 horas e como posso alterá-lo para 6? Alguém pode me dizer?
Ismoh
20

Uma abordagem alternativa interessante para o seu problema seria encontrar o mínimo aproximado para o Travelling Salesman Problem (TSP), ie. a rota mais curta que liga todos os seus pontos. Se os seus pontos formarem uma forma convexa, ela deve ser a solução certa; caso contrário, ainda deve ter uma boa aparência (uma forma "sólida" pode ser definida como uma que possui uma baixa relação perímetro / área, que é o que estamos otimizando aqui) .

Você pode usar qualquer implementação de um otimizador para o TSP, do qual tenho certeza de que pode encontrar muito no idioma de sua escolha.

static_rtti
fonte
Caramba. "Interessante" é um eufemismo. :)
Iterator
@Iterator: Fiquei muito feliz com a minha ideia, fiquei muito decepcionado por ter recebido votos negativos: - / Você acha que é válido?
static_rtti
1
Eu estava sugerindo o uso de uma das muitas aproximações rápidas, não o algoritmo original completo do NP, é claro.
static_rtti
6
Eu aprecio o ângulo adicional! Ter várias respostas válidas, se bem diferentes, pode ser de grande ajuda se alguém no futuro tropeçar nesse tópico procurando discutir opções.
Philipp Lenssen
1
Observe que minha abordagem é provavelmente mais lenta, mas mais correta em casos complexos: imagine o caso em que os pontos para um "8", por exemplo. As coordenadas polares não vão ajudá-lo nesse caso, e o resultado que você obterá dependerá muito do centro escolhido. A solução TSP é independente de qualquer parâmetro "heurístico".
static_rtti
19

O que você está pedindo é um sistema conhecido como coordenadas polares . A conversão de coordenadas cartesianas em polares é feita facilmente em qualquer idioma. As fórmulas podem ser encontradas nesta seção .

Depois de converter para coordenadas polares, basta classificar pelo ângulo, theta.

Iterador
fonte
4
Isso funcionará, mas também terá o defeito de fazer mais cálculos do que o necessário para responder à pergunta sobre pedidos. Na prática, você realmente não se importa com os ângulos reais ou as distâncias radiais, apenas com a ordem relativa deles. A solução de ciamej é melhor porque evita divisões, raízes quadradas e trigonométricas.
RBerteig
1
Não sei qual é o seu critério para "melhor". Por exemplo, comparar todos os pontos entre si é uma espécie de desperdício de computação. Trig não é algo que assusta os adultos, é?
Iterator
3
Não é que trigonometria seja assustadora. A questão é que o trig é caro para calcular e não era necessário para determinar a ordem relativa dos ângulos. Da mesma forma, você não precisa usar as raízes quadradas para colocar os raios em ordem. Uma conversão completa de coordenadas cartesianas em polares fará uma tangente ao arco e uma raiz quadrada. Portanto, sua resposta está correta, mas no contexto de computação gráfica ou geometria computacional é provável que não seja a melhor maneira de fazê-lo.
RBerteig
Entendi. No entanto, o OP não publicou como comp-geo, que foi uma tag de outra pessoa. Ainda assim, parece que a outra solução é polinomial no número de pontos, ou estou enganado? Nesse caso, isso queima mais ciclos do que trigonométricas.
Iterator
Na verdade, eu não havia notado a etiqueta comp-geo, apenas presumi que as únicas aplicações racionais para a pergunta deviam ser uma ou outra. Afinal, a questão do desempenho se torna discutível se houver apenas alguns pontos, e / ou a operação será feita raramente o suficiente. Nesse ponto, saber como fazê-lo se torna importante e é por isso que concordo que sua resposta está correta. Explica como calcular a noção de "ordem no sentido horário" em termos que podem ser explicados a praticamente qualquer pessoa.
RBerteig
3

Outra versão (retorne true se a vier antes de b no sentido anti-horário):

    bool lessCcw(const Vector2D &center, const Vector2D &a, const Vector2D &b) const
    {
        // Computes the quadrant for a and b (0-3):
        //     ^
        //   1 | 0
        //  ---+-->
        //   2 | 3

        const int dax = ((a.x() - center.x()) > 0) ? 1 : 0;
        const int day = ((a.y() - center.y()) > 0) ? 1 : 0;
        const int qa = (1 - dax) + (1 - day) + ((dax & (1 - day)) << 1);

        /* The previous computes the following:

           const int qa =
           (  (a.x() > center.x())
            ? ((a.y() > center.y())
                ? 0 : 3)
            : ((a.y() > center.y())
                ? 1 : 2)); */

        const int dbx = ((b.x() - center.x()) > 0) ? 1 : 0;
        const int dby = ((b.y() - center.y()) > 0) ? 1 : 0;
        const int qb = (1 - dbx) + (1 - dby) + ((dbx & (1 - dby)) << 1);

        if (qa == qb) {
            return (b.x() - center.x()) * (a.y() - center.y()) < (b.y() - center.y()) * (a.x() - center.x());
        } else {
            return qa < qb;
       } 
    }

Isso é mais rápido, porque o compilador (testado no Visual C ++ 2015) não gera salto para calcular dax, day, dbx, dby. Aqui o conjunto de saída do compilador:

; 28   :    const int dax = ((a.x() - center.x()) > 0) ? 1 : 0;

    vmovss  xmm2, DWORD PTR [ecx]
    vmovss  xmm0, DWORD PTR [edx]

; 29   :    const int day = ((a.y() - center.y()) > 0) ? 1 : 0;

    vmovss  xmm1, DWORD PTR [ecx+4]
    vsubss  xmm4, xmm0, xmm2
    vmovss  xmm0, DWORD PTR [edx+4]
    push    ebx
    xor ebx, ebx
    vxorps  xmm3, xmm3, xmm3
    vcomiss xmm4, xmm3
    vsubss  xmm5, xmm0, xmm1
    seta    bl
    xor ecx, ecx
    vcomiss xmm5, xmm3
    push    esi
    seta    cl

; 30   :    const int qa = (1 - dax) + (1 - day) + ((dax & (1 - day)) << 1);

    mov esi, 2
    push    edi
    mov edi, esi

; 31   : 
; 32   :    /* The previous computes the following:
; 33   : 
; 34   :    const int qa =
; 35   :        (   (a.x() > center.x())
; 36   :         ? ((a.y() > center.y()) ? 0 : 3)
; 37   :         : ((a.y() > center.y()) ? 1 : 2));
; 38   :    */
; 39   : 
; 40   :    const int dbx = ((b.x() - center.x()) > 0) ? 1 : 0;

    xor edx, edx
    lea eax, DWORD PTR [ecx+ecx]
    sub edi, eax
    lea eax, DWORD PTR [ebx+ebx]
    and edi, eax
    mov eax, DWORD PTR _b$[esp+8]
    sub edi, ecx
    sub edi, ebx
    add edi, esi
    vmovss  xmm0, DWORD PTR [eax]
    vsubss  xmm2, xmm0, xmm2

; 41   :    const int dby = ((b.y() - center.y()) > 0) ? 1 : 0;

    vmovss  xmm0, DWORD PTR [eax+4]
    vcomiss xmm2, xmm3
    vsubss  xmm0, xmm0, xmm1
    seta    dl
    xor ecx, ecx
    vcomiss xmm0, xmm3
    seta    cl

; 42   :    const int qb = (1 - dbx) + (1 - dby) + ((dbx & (1 - dby)) << 1);

    lea eax, DWORD PTR [ecx+ecx]
    sub esi, eax
    lea eax, DWORD PTR [edx+edx]
    and esi, eax
    sub esi, ecx
    sub esi, edx
    add esi, 2

; 43   : 
; 44   :    if (qa == qb) {

    cmp edi, esi
    jne SHORT $LN37@lessCcw

; 45   :        return (b.x() - center.x()) * (a.y() - center.y()) < (b.y() - center.y()) * (a.x() - center.x());

    vmulss  xmm1, xmm2, xmm5
    vmulss  xmm0, xmm0, xmm4
    xor eax, eax
    pop edi
    vcomiss xmm0, xmm1
    pop esi
    seta    al
    pop ebx

; 46   :    } else {
; 47   :        return qa < qb;
; 48   :    }
; 49   : }

    ret 0
$LN37@lessCcw:
    pop edi
    pop esi
    setl    al
    pop ebx
    ret 0
?lessCcw@@YA_NABVVector2D@@00@Z ENDP            ; lessCcw

Aproveitar.

AGPX
fonte
1
As duas instruções de retorno no comutador são matematicamente equivalentes. Existe uma razão para ter a opção?
unagi
0
  • vetor3 a = novo vetor3 (1, 0, 0) .............. wrt X_axis
  • vector3 b = qualquer ponto - Centro;
- y = |a * b|   ,   x =  a . b

- Atan2(y , x)...............................gives angle between -PI  to  + PI  in radians
- (Input % 360  +  360) % 360................to convert it from  0 to 2PI in radians
- sort by adding_points to list_of_polygon_verts by angle  we got 0  to 360

Finalmente, você recebe verts classificados no Anticlockwize

list.Reverse () .................. Orderwise_order

Pavan
fonte