Triangulação de Delaunay e diagrama de Voronoi

Considere um conjunto $\{p_i\}$ de pontos no plano.

Um diagrama de Voronoi $V(\{p_i\})$ de $\{p_i\}$ é uma partição do plano em $n$ regiões $V_i$, onde $V_i = \{p\in\mathbb{R}^2;\ \rho(p, p_i) = \min\ \rho(p, p_k)\}$. As células do diagrama de Voronoi são polígonos (possivelmente infinitos).

Uma triangulação de Delaunay $D(\{p_i\})$ de $\{p_i\}$ é uma triangulação em que todo ponto $p_i$ está fora ou no limite da circunferência de cada triângulo $T \in D(\{p_i\})$.

Há um caso desagradável e degenerado quando o diagrama de Voronoi não está conectado e a triangulação de Delaunay não existe. Nesse caso, todos os pontos são colineares.

Propriedades

A triangulação de Delaunay maximiza o ângulo mínimo entre todas as triangulações possíveis.

A árvore geradora mínima euclidiana de um conjunto de pontos é um subconjunto de arestas de sua triangulação de Delaunay.

Dualidade

Suponha que $\{p_i\}$ não seja colinear e entre $\{p_i\}$ não haja quatro pontos em um círculo. Então $V(\{p_i\})$ e $D(\{p_i\})$ são duais, portanto, se obtivermos um deles, podemos obter o outro em $O(n)$. O que fazer se não for o caso? O caso colinear pode ser processado facilmente. Caso contrário $V$ e $D'$ são duais, onde $D'$ é obtido a partir de $D$ removendo todas as arestas, de modo que dois triângulos nessa aresta compartilhem a circunferência.

Construindo Delaunay e Voronoi

Por causa da dualidade, precisamos apenas de um algoritmo rápido para calcular apenas um $V$ e $D$. Vamos descrever como criar $D(\{p_i\})$ em $O(n\log n)$. A triangulação será construída através do algoritmo de dividir e conquistar (Guibas e Stolfi).

Quad-edge - estrutura de dados

Durante o algoritmo, $D$ será armazenado dentro da estrutura de dados quad-edge data. Essa estrutura é descrita na figura:

Quad-Edge

No algoritmo, usaremos as seguintes funções nas arestas:

  1. make_edge(a, b)
    Essa função cria uma aresta isolada do ponto a ao ponto b juntamente com sua aresta reversa e as duas arestas duplas.
  2. splice(a, b)
    Esta é uma função chave do algoritmo. Ela troca(swap) a->Onext com b->Onext e a->Onext->Rot->Onext com b->Onext->Rot->Onext.
  3. delete_edge(e)
    Esta função exclui e da triangulação. Para excluir e, podemos simplesmente chamar splice(e, e->Oprev) e splice(e->Rev, e->Rev->Oprev).
  4. connect(a, b)
    Esta função cria uma nova aresta e de a->Dest para b->Org de modo que a, b, e todos tenham a mesma face esquerda. Para fazer isso, chamamos e = make_edge(a->Dest, b->Org), splice(e, a->Lnext) e splice(e->Rev, b).

Algoritmo

O algoritmo calculará a triangulação e retornará duas quad-edges: a aresta anti-horária do convex hull a partir do vértice mais à esquerda e a aresta horária do convex hull a partir do vértice mais à direita.

Vamos ordenar todos os pontos por x, e se $x_1 = x_2$ então por y. Vamos resolver o problema para algum segmento $(l, r)$ (inicialmente $(l, r) = (0, n - 1)$). Se $r - l + 1 = 2$, adicionaremos uma aresta $(p[l], p[r])$ e retornaremos. Se $r - l + 1 = 3$, primeiro adicionaremos as arestas $(p[l], p[l + 1])$ e $(p[l + 1], p[r])$. Também devemos conectá-los usando splice(a->Rev, b). Agora devemos fechar o triângulo. Nossa próxima ação dependerá da orientação de $p[l], p[l + 1], p[r]$. Se eles são colineares, não podemos fazer um triângulo, então simplesmente retornamos (a, b->Rev). Caso contrário, criamos uma nova aresta c chamando connect(b, a). Se os pontos são orientados no sentido anti-horário, retornamos (a, b->Rev). Caso contrário, retornamos (c->Rev, c).

Agora suponha que $r - l + 1 \ge 4$. Primeiramente, vamos resolver $L = (l, \frac{l + r}{2})$ e $R = (\frac{l + r}{2} + 1, r)$ recursivamente. Agora temos que mesclar essas triangulações em uma triangulação. Observe que nossos pontos são ordenados; portanto, durante a mesclagem, adicionaremos arestas de L a R (" cross edges ") e removeremos algumas arestas de L a L e de R a R. Qual é a estrutura das cross edges? Todas essas arestas devem cruzar uma linha paralela ao eixo y e posicionadas no valor x de divisão. Isso estabelece uma ordem linear das cross edges, então podemos falar sobre cross edges sucessivas, a cross edge mais baixa, etc. O algoritmo adicionará as cross edges em ordem crescente. Observe que quaisquer duas cross edges adjacentes terão um ponto final comum, e o terceiro lado do triângulo que elas definem vai de L a L ou de R a R. Vamos chamar a atual cross edge de base. O sucessor da base irá do ponto final esquerdo da base para um dos R-vizinhos do ponto final direito ou vice-versa. Considere a circunferência da base e a cross edge anterior. Suponha que esse círculo seja transformado em outros círculos, tendo a base como uma corda, mas bem distante na direção Oy. Nosso círculo aumentará por um tempo, mas, a menos que a base seja uma tangente superior de L e R, encontraremos um ponto pertencente a L ou a R, dando origem a um novo triângulo sem pontos na circunferência. A nova aresta L-R deste triângulo é a próxima cross edge adicionada. Para fazer isso com eficiência, calculamos duas arestas lcand e rcand de modo que lcand aponte para o primeiro ponto L encontrado nesse processo e rcand aponte para o primeiro ponto R. Então escolhemos o que seria encontrado primeiro. Inicialmente os pontos base para a tangente inferior de L e R.

Implementação

A implementação da função in_circle é específica do GCC.

typedef long long ll;

bool ge(const ll& a, const ll& b) { return a >= b; }
bool le(const ll& a, const ll& b) { return a <= b; }
bool eq(const ll& a, const ll& b) { return a == b; }
bool gt(const ll& a, const ll& b) { return a > b; }
bool lt(const ll& a, const ll& b) { return a < b; }
int sgn(const ll& a) { return a >= 0 ? a ? 1 : 0 : -1; }

struct pt {
    ll x, y;
    pt() { }
    pt(ll _x, ll _y) : x(_x), y(_y) { }
    pt operator-(const pt& p) const {
        return pt(x - p.x, y - p.y);
    }
    ll cross(const pt& p) const {
        return x * p.y - y * p.x;
    }
    ll cross(const pt& a, const pt& b) const {
        return (a - *this).cross(b - *this);
    }
    ll dot(const pt& p) const {
        return x * p.x + y * p.y;
    }
    ll dot(const pt& a, const pt& b) const {
        return (a - *this).dot(b - *this);
    }
    ll sqrLength() const {
        return this->dot(*this);
    }
    bool operator==(const pt& p) const {
        return eq(x, p.x) && eq(y, p.y);
    }
};

const pt inf_pt = pt(1e18, 1e18);

struct QuadEdge {
    pt origin;
    QuadEdge* rot = nullptr;
    QuadEdge* onext = nullptr;
    bool used = false;
    QuadEdge* rev() const {
        return rot->rot;
    }
    QuadEdge* lnext() const {
        return rot->rev()->onext->rot;
    }
    QuadEdge* oprev() const {
        return rot->onext->rot;
    }
    pt dest() const {
        return rev()->origin;
    }
};

QuadEdge* make_edge(pt from, pt to) {
    QuadEdge* e1 = new QuadEdge;
    QuadEdge* e2 = new QuadEdge;
    QuadEdge* e3 = new QuadEdge;
    QuadEdge* e4 = new QuadEdge;
    e1->origin = from;
    e2->origin = to;
    e3->origin = e4->origin = inf_pt;
    e1->rot = e3;
    e2->rot = e4;
    e3->rot = e2;
    e4->rot = e1;
    e1->onext = e1;
    e2->onext = e2;
    e3->onext = e4;
    e4->onext = e3;
    return e1;
}

void splice(QuadEdge* a, QuadEdge* b) {
    swap(a->onext->rot->onext, b->onext->rot->onext);
    swap(a->onext, b->onext);
}

void delete_edge(QuadEdge* e) {
    splice(e, e->oprev());
    splice(e->rev(), e->rev()->oprev());
    delete e->rot;
    delete e->rev()->rot;
    delete e;
    delete e->rev();
}

QuadEdge* connect(QuadEdge* a, QuadEdge* b) {
    QuadEdge* e = make_edge(a->dest(), b->origin);
    splice(e, a->lnext());
    splice(e->rev(), b);
    return e;
}

bool left_of(pt p, QuadEdge* e) {
    return gt(p.cross(e->origin, e->dest()), 0);
}

bool right_of(pt p, QuadEdge* e) {
    return lt(p.cross(e->origin, e->dest()), 0);
}

template <class T>
T det3(T a1, T a2, T a3, T b1, T b2, T b3, T c1, T c2, T c3) {
    return a1 * (b2 * c3 - c2 * b3) - a2 * (b1 * c3 - c1 * b3) +
           a3 * (b1 * c2 - c1 * b2);
}

bool in_circle(pt a, pt b, pt c, pt d) {
// Se tiver __int128, calcular diretamente.
// Caso contrário, calcular ângulos.
#if defined(__LP64__) || defined(_WIN64)
    __int128 det = -det3<__int128>(b.x, b.y, b.sqrLength(), c.x, c.y,
                                   c.sqrLength(), d.x, d.y, d.sqrLength());
    det += det3<__int128>(a.x, a.y, a.sqrLength(), c.x, c.y, c.sqrLength(), d.x,
                          d.y, d.sqrLength());
    det -= det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), d.x,
                          d.y, d.sqrLength());
    det += det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), c.x,
                          c.y, c.sqrLength());
    return det > 0;
#else
    auto ang = [](pt l, pt mid, pt r) {
        ll x = mid.dot(l, r);
        ll y = mid.cross(l, r);
        long double res = atan2((long double)x, (long double)y);
        return res;
    };
    long double kek = ang(a, b, c) + ang(c, d, a) - ang(b, c, d) - ang(d, a, b);
    if (kek > 1e-8)
        return true;
    else
        return false;
#endif
}

pair<QuadEdge*, QuadEdge*> build_tr(int l, int r, vector<pt>& p) {
    if (r - l + 1 == 2) {
        QuadEdge* res = make_edge(p[l], p[r]);
        return make_pair(res, res->rev());
    }
    if (r - l + 1 == 3) {
        QuadEdge *a = make_edge(p[l], p[l + 1]), *b = make_edge(p[l + 1], p[r]);
        splice(a->rev(), b);
        int sg = sgn(p[l].cross(p[l + 1], p[r]));
        if (sg == 0)
            return make_pair(a, b->rev());
        QuadEdge* c = connect(b, a);
        if (sg == 1)
            return make_pair(a, b->rev());
        else
            return make_pair(c->rev(), c);
    }
    int mid = (l + r) / 2;
    QuadEdge *ldo, *ldi, *rdo, *rdi;
    tie(ldo, ldi) = build_tr(l, mid, p);
    tie(rdi, rdo) = build_tr(mid + 1, r, p);
    while (true) {
        if (left_of(rdi->origin, ldi)) {
            ldi = ldi->lnext();
            continue;
        }
        if (right_of(ldi->origin, rdi)) {
            rdi = rdi->rev()->onext;
            continue;
        }
        break;
    }
    QuadEdge* basel = connect(rdi->rev(), ldi);
    auto valid = [&basel](QuadEdge* e) { return right_of(e->dest(), basel); };
    if (ldi->origin == ldo->origin)
        ldo = basel->rev();
    if (rdi->origin == rdo->origin)
        rdo = basel;
    while (true) {
        QuadEdge* lcand = basel->rev()->onext;
        if (valid(lcand)) {
            while (in_circle(basel->dest(), basel->origin, lcand->dest(),
                             lcand->onext->dest())) {
                QuadEdge* t = lcand->onext;
                delete_edge(lcand);
                lcand = t;
            }
        }
        QuadEdge* rcand = basel->oprev();
        if (valid(rcand)) {
            while (in_circle(basel->dest(), basel->origin, rcand->dest(),
                             rcand->oprev()->dest())) {
                QuadEdge* t = rcand->oprev();
                delete_edge(rcand);
                rcand = t;
            }
        }
        if (!valid(lcand) && !valid(rcand))
            break;
        if (!valid(lcand) ||
            (valid(rcand) && in_circle(lcand->dest(), lcand->origin,
                                       rcand->origin, rcand->dest())))
            basel = connect(rcand, basel->rev());
        else
            basel = connect(basel->rev(), lcand->rev());
    }
    return make_pair(ldo, rdo);
}

vector<tuple<pt, pt, pt>> delaunay(vector<pt> p) {
    sort(p.begin(), p.end(), [](const pt& a, const pt& b) {
        return lt(a.x, b.x) || (eq(a.x, b.x) && lt(a.y, b.y));
    });
    auto res = build_tr(0, (int)p.size() - 1, p);
    QuadEdge* e = res.first;
    vector<QuadEdge*> edges = {e};
    while (lt(e->onext->dest().cross(e->dest(), e->origin), 0))
        e = e->onext;
    auto add = [&p, &e, &edges]() {
        QuadEdge* curr = e;
        do {
            curr->used = true;
            p.push_back(curr->origin);
            edges.push_back(curr->rev());
            curr = curr->lnext();
        } while (curr != e);
    };
    add();
    p.clear();
    int kek = 0;
    while (kek < (int)edges.size()) {
        if (!(e = edges[kek++])->used)
            add();
    }
    vector<tuple<pt, pt, pt>> ans;
    for (int i = 0; i < (int)p.size(); i += 3) {
        ans.push_back(make_tuple(p[i], p[i + 1], p[i + 2]));
    }
    return ans;
}