Localização de um ponto em O (log n)

Considere o seguinte problema: você recebe um grafo planar sem vértices de grau um e zero, e um várias queries(testes). Cada query é um ponto, para o qual devemos determinar a face da subdivisão à qual pertence. Responderemos a cada query em $O(\log n)$.
Esse problema pode surgir quando você precisar localizar alguns pontos em um diagrama de Voronoi ou em algum polígono simples.

Algoritmo

Primeiramente, para cada ponto de query $p\ (x_0, y_0)$ queremos encontrar uma aresta que, se o ponto pertencer a alguma aresta, o ponto pertence a aresta que encontramos, caso contrário, essa aresta deve cruzar a linha $x = x_0$ em algum ponto único $(x_0, y)$ onde $y < y_0$ e esse $y$ é máximo dentre todas essas arestas. A imagem a seguir mostra os dois casos.

Image of Goal

Iremos resolver esse problema usando o algoritmo sweep line. Vamos iterar sobre a coordenadas x dos pontos da query e os pontos finais das arestas em ordem crescente e mantendo um conjunto de arestas $s$. Para coordenada x iremos adcionar alguns eventos antecipadamente.

Os eventos serão de quatro tipos: add, remove, vertical, get. Para cada aresta vertical (ambos pontos finais tem o mesmo x) adcionaremos um evento vertical para a coordenada x correspondente. Para todas as outras arestas adcionaremos um evento add para o mínimo de coordenadas x dos pontos finais e um evento remove para o máximo de coordenadas x dos pontos finais. Finalmente, para cada ponto da query adcionaremos um evento get para sua coordenada x.

Para cada coordenada x, ordenaremos os eventos por seus tipos: (vertical, get, remove, add). A imagem a seguir mostra todos os eventos em ordem para cada coordenada x.

Image of Events

Manteremos dois conjuntos durante o processo da sweep line. Um conjunto $ t $ para todas as arestas não verticais e um conjunto $ vert $, especialmente para as verticais. Limparemos o conjunto $vert$ no começo do processamento de cada coordenada x.

Agora iremos processar os eventos para uma coordenada x fixa.

Agora vamos escolher um comparador para o conjunto $t$. Esse comparador deve verificar se uma aresta não está acima da outra para cada coordenada x que ambas abrangem. Suponha que tenhamos duas arestas $(a, b)$ e $(c, d)$. Então o comparador é (em pseudocódigo):

$val = ((b - a)\times(c - a)) + ((b - a)\times(d - a))$
if $val \neq 0$
return $val > 0$
$val = ((d - c)\times(a - c)) + ((d - c)\times(b - c))$
return $val < 0$

Agora para cada query temos a aresta correspondente. Como encontrar o plano ? Se não conseguimos encontrar a aresta, significa que o ponto está na face externa. Se o ponto pertence à aresta que encontramos, a face não é única. Caso contrário, existem dois casos - as faces que são delimitadas por essa aresta. Como verificar qual é a resposta? Observe que a aresta não é vertical. Então a resposta é o rosto que está acima dessa aresta. Vamos encontrar tal face para cada aresta não vertical. Considere um deslocamento no sentido anti-horário de cada face. Se durante esse percurso aumentamos a coordenada x ao passar pela aresta, então essa face é a face que precisamos encontrar para essa aresta.

Notas

Na verdade, com árvores persistentes, essa abordagem pode ser usada para responder às queries.

Implementação

O código a seguir é implementado para números inteiros, mas pode ser modificado para funcionar com doubles (alterando os métodos de comparação e o tipo do ponto). Esta implementação pressupõe que o plano seja armazenado corretamente dentro de uma DCEL e a face externa seja numerada $ -1 $(não se encontra em uma face preferível).
Para cada query um par $(1, i)$ é retornado se o ponto estiver estritamente dentro do número da face $i$, e um par $(0, i)$ é retornado se o ponto estiver no número da aresta $i$.

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

    // coletando todas as coordenadas x
    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);
    }

    // mapear todas as coordenadas x para 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++;

    // criando eventos
    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});
        }
    }

    // executar o algoritmo sweep line
    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

TIMUS1848 - Fly Hunt