Skip to Content
CTRL K

Triangulación de Delaunay y diagrama de Voronoi

Consideremos un conjunto {pi}{p_i} de puntos en el plano. Un diagrama de Voronoi V({pi})V({p_i}) de {pi}{p_i} es una partición del plano en nn regiones ViV_i, donde Vi={p∈R2; ρ(p,pi)=min⁡ ρ(p,pk)}V_i = {p\in\mathbb{R}^2;\ \rho(p, p_i) = \min\ \rho(p, p_k)}. Las celdas del diagrama de Voronoi son polígonos (posiblemente infinitos). Una triangulación de Delaunay D({pi})D({p_i}) de {pi}{p_i} es una triangulación donde cada punto pip_i está fuera o en el borde de la circunferencia circunscrita de cada triángulo T∈D({pi})T \in D({p_i}).

Hay un caso degenerado desagradable cuando el diagrama de Voronoi no es conexo y la triangulación de Delaunay no existe. Este caso es cuando todos los puntos son colineales.

Propiedades

La triangulación de Delaunay maximiza el ángulo mínimo entre todas las triangulaciones posibles.

El árbol de expansión euclidiano mínimo de un conjunto de puntos es un subconjunto de aristas de su triangulación de Delaunay.

Dualidad

Supongamos que {pi}{p_i} no es colineal y entre {pi}{p_i} no hay cuatro puntos que yacen en un círculo. Entonces V({pi})V({p_i}) y D({pi})D({p_i}) son duales, así que si obtenemos uno de ellos, podemos obtener el otro en O(n)O(n). ¿Qué hacer si no es el caso? El caso colineal se puede procesar fácilmente. En caso contrario, VV y D′D’ son duales, donde D′D’ se obtiene de DD quitando todas las aristas tales que dos triángulos sobre esta arista comparten la circunferencia circunscrita.

Construir Delaunay y Voronoi

Por la dualidad, solo necesitamos un algoritmo rápido para computar solo uno de VV y DD. Describiremos cómo construir D({pi})D({p_i}) en O(nlog⁡n)O(n\log n). La triangulación se construirá vía el algoritmo de divide y vencerás debido a Guibas y Stolfi.

Estructura de datos quad-edge

Durante el algoritmo DD se guardará dentro de la estructura de datos quad-edge. Esta estructura se describe en la imagen:

Quad-Edge

En el algoritmo usaremos las siguientes funciones sobre aristas:

  1. make_edge(a, b)
    Esta función crea una arista aislada del punto a al punto b junto con su arista inversa y ambas aristas duales.
  2. splice(a, b)
    Esta es una función clave del algoritmo. Intercambia a->Onext con b->Onext y a->Onext->Rot->Onext con b->Onext->Rot->Onext.
  3. delete_edge(e)
    Esta función borra e de la triangulación. Para borrar e, podemos simplemente llamar splice(e, e->Oprev) y splice(e->Rev, e->Rev->Oprev).
  4. connect(a, b)
    Esta función crea una arista nueva e de a->Dest a b->Org de tal forma que a, b, e todos tienen la misma cara izquierda. Para hacer esto, llamamos e = make_edge(a->Dest, b->Org), splice(e, a->Lnext) y splice(e->Rev, b).

Algoritmo

El algoritmo computará la triangulación y devolverá dos quad-edges: la arista de la envolvente convexa en sentido antihorario que sale del vértice más a la izquierda y la arista de la envolvente convexa en sentido horario que sale del vértice más a la derecha.

Ordenemos todos los puntos por x, y si x1=x2x_1 = x_2 entonces por y. Resolvamos el problema para algún segmento (l,r)(l, r) (inicialmente (l,r)=(0,n−1)(l, r) = (0, n - 1)). Si r−l+1=2r - l + 1 = 2, agregaremos una arista (p[l],p[r])(p[l], p[r]) y devolveremos. Si r−l+1=3r - l + 1 = 3, primero agregaremos las aristas (p[l],p[l+1])(p[l], p[l + 1]) y (p[l+1],p[r])(p[l + 1], p[r]). También debemos conectarlas usando splice(a->Rev, b). Ahora debemos cerrar el triángulo. Nuestra siguiente acción dependerá de la orientación de p[l],p[l+1],p[r]p[l], p[l + 1], p[r]. Si son colineales, no podemos hacer un triángulo, así que simplemente devolvemos (a, b->Rev). En caso contrario, creamos una arista nueva c llamando connect(b, a). Si los puntos están orientados en sentido antihorario, devolvemos (a, b->Rev). En caso contrario devolvemos (c->Rev, c).

Ahora supongamos que r−l+1≥4r - l + 1 \ge 4. Primero, resolvamos L=(l,l+r2)L = (l, \frac{l + r}{2}) y R=(l+r2+1,r)R = (\frac{l + r}{2} + 1, r) de forma recursiva. Ahora tenemos que fusionar estas triangulaciones en una sola triangulación. Nótese que nuestros puntos están ordenados, así que al fusionar agregaremos aristas de L a R (las llamadas aristas cruzadas) y quitaremos algunas aristas de L a L y de R a R. ¿Cuál es la estructura de las aristas cruzadas? Todas estas aristas deben cruzar una recta paralela al eje y y colocada en el valor x de partición. Esto establece un orden lineal de las aristas cruzadas, así que podemos hablar de aristas cruzadas sucesivas, la arista cruzada más inferior, etc. El algoritmo agregará las aristas cruzadas en orden ascendente. Nótese que cualesquiera dos aristas cruzadas adyacentes tendrán un extremo común, y el tercer lado del triángulo que definen va de L a L o de R a R. Llamemos a la arista cruzada actual la base. El sucesor de la base o bien irá del extremo izquierdo de la base a uno de los vecinos-R del extremo derecho o viceversa. Consideremos la circunferencia circunscrita de la base y la arista cruzada anterior. Supongamos que este círculo se transforma en otros círculos que tienen la base como cuerda pero yacen más adelante en la dirección Oy. Nuestro círculo subirá por un rato, pero a menos que la base sea una tangente superior de L y R encontraremos un punto que pertenece o bien a L o a R dando lugar a un triángulo nuevo sin ningún punto en la circunferencia circunscrita. La nueva arista L-R de este triángulo es la siguiente arista cruzada agregada. Para hacer esto de forma eficiente, computamos dos aristas lcand y rcand de modo que lcand apunta al primer punto L encontrado en este proceso, y rcand apunta al primer punto R. Luego elegimos el que se encontraría primero. Inicialmente la base apunta a la tangente inferior de L y R.

Implementación

Nótese que la implementación de la función in_circle es específica de 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->rev()->rot;
    delete e->rev();
    delete e->rot;
    delete e;
}

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) {
// If there is __int128, calculate directly.
// Otherwise, calculate angles.
#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;
}

Problemas

Intersección de semiplanosTruco de la envolvente convexa