Skip to Content
CTRL K

Localización de puntos en O(logn)O(log n)

Consideremos el siguiente problema: nos dan una subdivisión planar  sin ningún vértice de grado uno ni cero, y muchas consultas. Cada consulta es un punto, para el cual debemos determinar la cara de la subdivisión a la que pertenece. Responderemos cada consulta en O(log⁡n)O(\log n) de forma offline.
Este problema puede surgir cuando hay que localizar algunos puntos en un diagrama de Voronoi o en algún polígono simple.

Algoritmo

Primero, para cada punto de consulta p (x0,y0)p\ (x_0, y_0) queremos encontrar una arista tal que si el punto pertenece a alguna arista, el punto yazca en la arista que encontramos; en caso contrario esta arista debe intersectar la recta x=x0x = x_0 en algún punto único (x0,y)(x_0, y) donde y<y0y < y_0 y este yy es máximo entre todas esas aristas. La siguiente imagen muestra ambos casos.

Imagen del objetivo

Resolveremos este problema de forma offline usando el algoritmo de línea de barrido. Iteremos sobre las coordenadas x de los puntos de consulta y de los extremos de las aristas en orden creciente y mantengamos un conjunto de aristas ss. Para cada coordenada x agregaremos algunos eventos de antemano.

Los eventos serán de cuatro tipos: add, remove, vertical, get. Para cada arista vertical (ambos extremos tienen la misma coordenada x) agregaremos un evento vertical para la coordenada x correspondiente. Para cada otra arista agregaremos un evento add para el mínimo de las coordenadas x de los extremos y un evento remove para el máximo de las coordenadas x de los extremos. Finalmente, para cada punto de consulta agregaremos un evento get para su coordenada x.

Para cada coordenada x ordenaremos los eventos por sus tipos en el orden (vertical, get, remove, add). La siguiente imagen muestra todos los eventos en orden para cada coordenada xx.

Imagen de eventos

Mantendremos dos conjuntos durante el proceso de línea de barrido. Un conjunto tt para todas las aristas no verticales, y un conjunto vertvert especialmente para las verticales. Limpiaremos el conjunto vertvert al principio de procesar cada coordenada x.

Ahora procesemos los eventos para una coordenada x fija.

  • Si obtuvimos un evento vertical, simplemente insertaremos la coordenada y mínima de los extremos de la arista correspondiente en vertvert.
  • Si obtuvimos un evento remove o add, quitaremos la arista correspondiente de tt o la agregaremos a tt.
  • Finalmente, para cada evento get debemos comprobar si el punto yace en alguna arista vertical realizando una búsqueda binaria en vertvert. Si el punto no yace en ninguna arista vertical, debemos encontrar la respuesta para esta consulta en tt. Para hacer esto, otra vez hacemos una búsqueda binaria. Para manejar algunos casos degenerados (p. ej. en el caso del triángulo (0, 0)(0,~0), (0, 2)(0,~2), (1,1)(1, 1) cuando consultamos el punto (0, 0)(0,~0)), debemos responder todos los eventos get otra vez después de procesar todos los eventos para esta coordenada x y elegir la mejor de las dos respuestas.

Ahora elijamos un comparador para el conjunto tt. Este comparador debería comprobar si una arista no yace por encima de otra para cada coordenada x que ambas cubren. Supongamos que tenemos dos aristas (a,b)(a, b) y (c,d)(c, d). Entonces el comparador es (en pseudocódigo):

val=sgn((b−a)×(c−a))+sgn((b−a)×(d−a))val = sgn((b - a)\times(c - a)) + sgn((b - a)\times(d - a))
if val≠0val \neq 0
then return val>0val > 0
val=sgn((d−c)×(a−c))+sgn((d−c)×(b−c))val = sgn((d - c)\times(a - c)) + sgn((d - c)\times(b - c))
return val<0val < 0

Ahora para cada consulta tenemos la arista correspondiente. ¿Cómo encontrar la cara? Si no pudimos encontrar la arista significa que el punto está en la cara exterior. Si el punto pertenece a la arista que encontramos, la cara no es única. En caso contrario, hay dos candidatos: las caras que están acotadas por esta arista. ¿Cómo comprobar cuál es la respuesta? Nótese que la arista no es vertical. Entonces la respuesta es la cara que está por encima de esta arista. Encontremos tal cara para cada arista no vertical. Consideremos un recorrido antihorario de cada cara. Si durante este recorrido aumentamos la coordenada x al pasar por la arista, entonces esta cara es la cara que necesitamos encontrar para esta arista.

Notas

De hecho, con árboles persistentes este enfoque se puede usar para responder las consultas de forma online.

Implementación

El siguiente código está implementado para enteros, pero se puede modificar fácilmente para trabajar con doubles (cambiando los métodos de comparación y el tipo de punto). Esta implementación asume que la subdivisión está guardada correctamente en una DCEL  (lista de aristas doblemente conectada) y que la cara exterior está numerada −1-1.
Para cada consulta se devuelve un par (1,i)(1, i) si el punto yace estrictamente dentro de la cara número ii, y un par (0,i)(0, i) si el punto yace en la arista número ii.

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& x) { return le(x, 0) ? eq(x, 0) ? 0 : -1 : 1; }

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

struct Edge {
    pt l, r;
};

bool edge_cmp(Edge* edge1, Edge* edge2)
{
    const pt a = edge1->l, b = edge1->r;
    const pt c = edge2->l, d = edge2->r;
    int val = sgn(a.cross(b, c)) + sgn(a.cross(b, d));
    if (val != 0)
        return val > 0;
    val = sgn(c.cross(d, a)) + sgn(c.cross(d, b));
    return val < 0;
}

enum EventType { DEL = 2, ADD = 3, GET = 1, VERT = 0 };

struct Event {
    EventType type;
    int pos;
    bool operator<(const Event& event) const { return type < event.type; }
};

vector<Edge*> sweepline(vector<Edge*> planar, vector<pt> queries)
{
    using pt_type = decltype(pt::x);

    // collect all x-coordinates
    auto s =
        set<pt_type, std::function<bool(const pt_type&, const pt_type&)>>(lt);
    for (pt p : queries)
        s.insert(p.x);
    for (Edge* e : planar) {
        s.insert(e->l.x);
        s.insert(e->r.x);
    }

    // map all x-coordinates to ids
    int cid = 0;
    auto id =
        map<pt_type, int, std::function<bool(const pt_type&, const pt_type&)>>(
            lt);
    for (auto x : s)
        id[x] = cid++;

    // create events
    auto t = set<Edge*, decltype(*edge_cmp)>(edge_cmp);
    auto vert_cmp = [](const pair<pt_type, int>& l,
                       const pair<pt_type, int>& r) {
        if (!eq(l.first, r.first))
            return lt(l.first, r.first);
        return l.second < r.second;
    };
    auto vert = set<pair<pt_type, int>, decltype(vert_cmp)>(vert_cmp);
    vector<vector<Event>> events(cid);
    for (int i = 0; i < (int)queries.size(); i++) {
        int x = id[queries[i].x];
        events[x].push_back(Event{GET, i});
    }
    for (int i = 0; i < (int)planar.size(); i++) {
        int lx = id[planar[i]->l.x], rx = id[planar[i]->r.x];
        if (lx > rx) {
            swap(lx, rx);
            swap(planar[i]->l, planar[i]->r);
        }
        if (lx == rx) {
            events[lx].push_back(Event{VERT, i});
        } else {
            events[lx].push_back(Event{ADD, i});
            events[rx].push_back(Event{DEL, i});
        }
    }

    // perform sweep line algorithm
    vector<Edge*> ans(queries.size(), nullptr);
    for (int x = 0; x < cid; x++) {
        sort(events[x].begin(), events[x].end());
        vert.clear();
        for (Event event : events[x]) {
            if (event.type == DEL) {
                t.erase(planar[event.pos]);
            }
            if (event.type == VERT) {
                vert.insert(make_pair(
                    min(planar[event.pos]->l.y, planar[event.pos]->r.y),
                    event.pos));
            }
            if (event.type == ADD) {
                t.insert(planar[event.pos]);
            }
            if (event.type == GET) {
                auto jt = vert.upper_bound(
                    make_pair(queries[event.pos].y, planar.size()));
                if (jt != vert.begin()) {
                    --jt;
                    int i = jt->second;
                    if (ge(max(planar[i]->l.y, planar[i]->r.y),
                           queries[event.pos].y)) {
                        ans[event.pos] = planar[i];
                        continue;
                    }
                }
                Edge* e = new Edge;
                e->l = e->r = queries[event.pos];
                auto it = t.upper_bound(e);
                if (it != t.begin())
                    ans[event.pos] = *(--it);
                delete e;
            }
        }

        for (Event event : events[x]) {
            if (event.type != GET)
                continue;
            if (ans[event.pos] != nullptr &&
                eq(ans[event.pos]->l.x, ans[event.pos]->r.x))
                continue;

            Edge* e = new Edge;
            e->l = e->r = queries[event.pos];
            auto it = t.upper_bound(e);
            delete e;
            if (it == t.begin())
                e = nullptr;
            else
                e = *(--it);
            if (ans[event.pos] == nullptr) {
                ans[event.pos] = e;
                continue;
            }
            if (e == nullptr)
                continue;
            if (e == ans[event.pos])
                continue;
            if (id[ans[event.pos]->r.x] == x) {
                if (id[e->l.x] == x) {
                    if (gt(e->l.y, ans[event.pos]->r.y))
                        ans[event.pos] = e;
                }
            } else {
                ans[event.pos] = e;
            }
        }
    }
    return ans;
}

struct DCEL {
    struct Edge {
        pt origin;
        Edge* nxt = nullptr;
        Edge* twin = nullptr;
        int face;
    };
    vector<Edge*> body;
};

vector<pair<int, int>> point_location(DCEL planar, vector<pt> queries)
{
    vector<pair<int, int>> ans(queries.size());
    vector<Edge*> planar2;
    map<intptr_t, int> pos;
    map<intptr_t, int> added_on;
    int n = planar.body.size();
    for (int i = 0; i < n; i++) {
        if (planar.body[i]->face > planar.body[i]->twin->face)
            continue;
        Edge* e = new Edge;
        e->l = planar.body[i]->origin;
        e->r = planar.body[i]->twin->origin;
        added_on[(intptr_t)e] = i;
        pos[(intptr_t)e] =
            lt(planar.body[i]->origin.x, planar.body[i]->twin->origin.x)
                ? planar.body[i]->face
                : planar.body[i]->twin->face;
        planar2.push_back(e);
    }
    auto res = sweepline(planar2, queries);
    for (int i = 0; i < (int)queries.size(); i++) {
        if (res[i] == nullptr) {
            ans[i] = make_pair(1, -1);
            continue;
        }
        pt p = queries[i];
        pt l = res[i]->l, r = res[i]->r;
        if (eq(p.cross(l, r), 0) && le(p.dot(l, r), 0)) {
            ans[i] = make_pair(0, added_on[(intptr_t)res[i]]);
            continue;
        }
        ans[i] = make_pair(1, pos[(intptr_t)res[i]]);
    }
    for (auto e : planar2)
        delete e;
    return ans;
}

Problemas

Longitud de la unión de segmentosEncontrar el par de puntos más cercano