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.
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.
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.
- Se tivermos um evento vertical, vamos simplesmente inserir a coordenada y mínima dos finais das arestas correspondentes para $vert$.
- Se tivermos um evento de remove ou add, removeremos a aresta correspondente de $t$ ou adcionaremos em $t$.
- Finalmente, para cada evento get devemos verificar se o ponto está em alguma aresta vertical executando uma pesquisa binária em $vert$. Se o ponto não estiver em nenhuma aresta vertical, devemos encontrar a resposta para essa query em $t$. Para fazer isso, novamente, executamos uma busca binária. Para lidar com alguns casos do tipo (e.g. um triângulo $(0,~0)$, $(0,~2)$, $(1, 1)$ quando recebemos o ponto na query $(0,~0)$), precisamos responder todos os eventos get novamente depois de processarmos todos os eventos dessa coordenada x e escolher a melhor das duas respostas.
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;
}