Construcción de la envolvente convexa
En este artículo discutiremos el problema de construir una envolvente convexa (convex hull) a partir de un conjunto de puntos.
Considérense puntos dados en un plano, y el objetivo es generar una envolvente convexa, es decir, el menor polígono convexo que contiene todos los puntos dados.
Veremos el algoritmo de barrido de Graham (Graham’s scan) publicado en 1972 por Graham, y también el algoritmo de cadena monótona (Monotone chain) publicado en 1979 por Andrew. Ambos son , y son asintóticamente óptimos (pues está demostrado que no hay un algoritmo asintóticamente mejor), con excepción de unos pocos problemas en los que intervienen procesamiento paralelo o online.
Algoritmo de barrido de Graham
El algoritmo primero encuentra el punto más bajo . Si hay varios puntos con la misma coordenada Y, se considera el de menor coordenada X. Este paso toma tiempo .
A continuación, todos los demás puntos se ordenan por ángulo polar en sentido horario. Si el ángulo polar entre dos o más puntos es el mismo, el empate debe resolverse por distancia a , en orden creciente.
Luego iteramos por cada punto uno a uno, y nos aseguramos de que el punto actual y los dos anteriores formen un giro horario; en caso contrario se descarta el punto anterior, ya que produciría una forma no convexa. Comprobar si el giro es horario o antihorario se puede hacer comprobando la orientación.
Usamos una pila para almacenar los puntos, y una vez que alcanzamos el punto original , el algoritmo termina y devolvemos la pila que contiene todos los puntos de la envolvente convexa en sentido horario.
Si se necesitan incluir los puntos colineales al hacer un barrido de Graham, hace falta otro paso después de ordenar. Hay que tomar los puntos que tienen la mayor distancia polar desde (estos deberían estar al final del vector ordenado) y son colineales. Los puntos de esta recta deben invertirse para poder emitir todos los puntos colineales; de lo contrario el algoritmo tomaría el punto más cercano de esta recta y se detendría. Este paso no debe incluirse en la versión no colineal del algoritmo; de lo contrario no se obtendría la menor envolvente convexa.
Implementación
struct pt {
double x, y;
bool operator == (pt const& t) const {
return x == t.x && y == t.y;
}
};
int orientation(pt a, pt b, pt c) {
double v = a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y);
if (v < 0) return -1; // horario
if (v > 0) return +1; // antihorario
return 0;
}
bool cw(pt a, pt b, pt c, bool include_collinear) {
int o = orientation(a, b, c);
return o < 0 || (include_collinear && o == 0);
}
bool collinear(pt a, pt b, pt c) { return orientation(a, b, c) == 0; }
void convex_hull(vector<pt>& a, bool include_collinear = false) {
pt p0 = *min_element(a.begin(), a.end(), [](pt a, pt b) {
return make_pair(a.y, a.x) < make_pair(b.y, b.x);
});
sort(a.begin(), a.end(), [&p0](const pt& a, const pt& b) {
int o = orientation(p0, a, b);
if (o == 0)
return (p0.x-a.x)*(p0.x-a.x) + (p0.y-a.y)*(p0.y-a.y)
< (p0.x-b.x)*(p0.x-b.x) + (p0.y-b.y)*(p0.y-b.y);
return o < 0;
});
if (include_collinear) {
int i = (int)a.size()-1;
while (i >= 0 && collinear(p0, a[i], a.back())) i--;
reverse(a.begin()+i+1, a.end());
}
vector<pt> st;
for (int i = 0; i < (int)a.size(); i++) {
while (st.size() > 1 && !cw(st[st.size()-2], st.back(), a[i], include_collinear))
st.pop_back();
st.push_back(a[i]);
}
if (include_collinear == false && st.size() == 2 && st[0] == st[1])
st.pop_back();
a = st;
}Algoritmo de cadena monótona
El algoritmo primero encuentra los puntos más a la izquierda y más a la derecha A y B. Si existen varios de esos puntos, se toma como A el más bajo entre los de la izquierda (menor coordenada Y), y como B el más alto entre los de la derecha (mayor coordenada Y). Claramente, A y B deben pertenecer ambos a la envolvente convexa, pues son los más alejados y no pueden estar contenidos por ninguna recta formada por un par entre los puntos dados.
Ahora se traza una recta por AB. Esto divide el resto de los puntos en dos conjuntos, S1 y S2, donde S1 contiene todos los puntos por encima de la recta que une A y B, y S2 contiene todos los puntos por debajo de la recta que une A y B. Los puntos que yacen sobre la recta que une A y B pueden pertenecer a cualquiera de los dos conjuntos. Los puntos A y B pertenecen a ambos conjuntos. Ahora el algoritmo construye el conjunto superior S1 y el conjunto inferior S2 y luego los combina para obtener la respuesta.
Para obtener el conjunto superior, ordenamos todos los puntos por la coordenada x. Para cada punto comprobamos si o bien el punto actual es el último punto (que definimos como B), o bien la orientación entre la recta entre A y el punto actual y la recta entre el punto actual y B es horaria. En esos casos el punto actual pertenece al conjunto superior S1. Comprobar si el giro es horario o antihorario se puede hacer comprobando la orientación.
Si el punto dado pertenece al conjunto superior, comprobamos el ángulo formado por la recta que une el penúltimo punto y el último punto de la envolvente convexa superior, con la recta que une el último punto de la envolvente convexa superior y el punto actual. Si el ángulo no es horario, quitamos el punto más reciente añadido a la envolvente convexa superior, pues el punto actual podrá contener al punto anterior una vez que se añada a la envolvente convexa.
La misma lógica se aplica al conjunto inferior S2. Si o bien el punto actual es B, o bien la orientación de las rectas formadas por A y el punto actual y por el punto actual y B es antihoraria, entonces pertenece a S2.
Si el punto dado pertenece al conjunto inferior, actuamos de forma similar a un punto del conjunto superior, salvo que comprobamos una orientación antihoraria en lugar de una orientación horaria. Así, si el ángulo formado por la recta que une el penúltimo punto y el último punto de la envolvente convexa inferior, con la recta que une el último punto de la envolvente convexa inferior y el punto actual no es antihorario, quitamos el punto más reciente añadido a la envolvente convexa inferior, pues el punto actual podrá contener el punto anterior una vez añadido a la envolvente.
La envolvente convexa final se obtiene de la unión de las envolventes convexas superior e inferior, formando una envolvente en sentido horario, y la implementación es la siguiente.
Si se necesitan puntos colineales, basta con comprobarlos en las rutinas horaria/antihoraria. Sin embargo, esto permite un caso degenerado en el que todos los puntos de entrada son colineales en una sola recta, y el algoritmo emitiría puntos repetidos. Para resolverlo, comprobamos si la envolvente superior contiene todos los puntos, y si es así, simplemente devolvemos los puntos en orden inverso, que es lo que devolvería la implementación de Graham en este caso.
Implementación
struct pt {
double x, y;
};
int orientation(pt a, pt b, pt c) {
double v = a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y);
if (v < 0) return -1; // horario
if (v > 0) return +1; // antihorario
return 0;
}
bool cw(pt a, pt b, pt c, bool include_collinear) {
int o = orientation(a, b, c);
return o < 0 || (include_collinear && o == 0);
}
bool ccw(pt a, pt b, pt c, bool include_collinear) {
int o = orientation(a, b, c);
return o > 0 || (include_collinear && o == 0);
}
void convex_hull(vector<pt>& a, bool include_collinear = false) {
if (a.size() == 1)
return;
sort(a.begin(), a.end(), [](pt a, pt b) {
return make_pair(a.x, a.y) < make_pair(b.x, b.y);
});
pt p1 = a[0], p2 = a.back();
vector<pt> up, down;
up.push_back(p1);
down.push_back(p1);
for (int i = 1; i < (int)a.size(); i++) {
if (i == a.size() - 1 || cw(p1, a[i], p2, include_collinear)) {
while (up.size() >= 2 && !cw(up[up.size()-2], up[up.size()-1], a[i], include_collinear))
up.pop_back();
up.push_back(a[i]);
}
if (i == a.size() - 1 || ccw(p1, a[i], p2, include_collinear)) {
while (down.size() >= 2 && !ccw(down[down.size()-2], down[down.size()-1], a[i], include_collinear))
down.pop_back();
down.push_back(a[i]);
}
}
if (include_collinear && up.size() == a.size()) {
reverse(a.begin(), a.end());
return;
}
a.clear();
for (int i = 0; i < (int)up.size(); i++)
a.push_back(up[i]);
for (int i = down.size() - 2; i > 0; i--)
a.push_back(down[i]);
}