Intersección de semiplanos

En este artículo discutiremos el problema de computar la intersección de un conjunto de semiplanos (half-planes). Tal intersección se puede representar de forma conveniente como una región/polígono convexo, donde cada punto dentro de él también está dentro de todos los semiplanos, y es este polígono el que estamos intentando encontrar o construir. Damos alguna intuición inicial para el problema, describimos un enfoque $O(N \log N)$ conocido como el algoritmo Sort-and-Incremental y damos algunas aplicaciones de ejemplo de esta técnica.

Se recomienda fuertemente que el lector esté familiarizado con primitivas y operaciones geométricas básicas (puntos, vectores, intersección de rectas). Además, el conocimiento sobre envolventes convexas o el truco de la envolvente convexa puede ayudar a entender mejor los conceptos de este artículo, pero no son un prerrequisito de ninguna forma.

Aclaraciones y definiciones iniciales

Para todo el artículo, haremos algunas hipótesis (salvo que se especifique lo contrario):

1. Definimos $N$ como la cantidad de semiplanos en el conjunto dado.
2. Representaremos rectas y semiplanos por un punto y un vector (cualquier punto que yazca en la recta dada, y el vector director de la recta). En el caso de semiplanos, asumimos que cada semiplano permite la región al lado izquierdo de su vector director. Además, definimos el ángulo de un semiplano como el ángulo polar de su vector director. Véase la imagen de abajo para un ejemplo.
3. Asumiremos que la intersección resultante es siempre o bien acotada o vacía. Si necesitamos manejar el caso no acotado, podemos simplemente agregar 4 semiplanos que definen una caja acotante lo bastante grande.
4. Asumiremos, por simplicidad, que no hay semiplanos paralelos en el conjunto dado. Hacia el final del artículo discutiremos cómo lidiar con tales casos.

El semiplano $y \geq 2x - 2$ se puede representar como el punto $P = (1, 0)$ con vector director $PQ = Q - P = (1, 2)$

Enfoque de fuerza bruta - $O(N^3)$ {data-toc-label=“Enfoque de fuerza bruta - O(N^3)”}

Una de las soluciones más directas y obvias sería computar el punto de intersección de las rectas de todos los pares de semiplanos y, para cada punto, comprobar si está dentro de todos los demás semiplanos. Como hay $O(N^2)$ puntos de intersección, y para cada uno de ellos hay que comprobar $O(N)$ semiplanos, la complejidad temporal total es $O(N^3)$. La región real de la intersección se puede reconstruir entonces usando, por ejemplo, un algoritmo de envolvente convexa sobre el conjunto de puntos de intersección que estaban incluidos en todos los semiplanos.

Es bastante fácil ver por qué esto funciona: los vértices del polígono convexo resultante son todos puntos de intersección de las rectas de los semiplanos, y cada uno de esos vértices es obviamente parte de todos los semiplanos. La principal ventaja de este método es que es fácil de entender, recordar y programar sobre la marcha si solo se necesita comprobar si la intersección está vacía o no. Sin embargo, es terriblemente lento e inadecuado para la mayoría de los problemas, así que necesitamos algo más rápido.

Enfoque incremental - $O(N^2)$ {data-toc-label=“Enfoque incremental - O(N^2)”}

Otro enfoque bastante directo es construir de forma incremental la intersección de los semiplanos, uno a la vez. Este método es básicamente equivalente a cortar un polígono convexo por una recta $N$ veces, y quitar los semiplanos redundantes en cada paso. Para hacer esto, podemos representar el polígono convexo como una lista de segmentos de recta, y para cortarlo con un semiplano simplemente encontramos los puntos de intersección de los segmentos con la recta del semiplano (solo habrá dos puntos de intersección si la recta intersecta el polígono de forma propia), y reemplazamos todos los segmentos de recta intermedios por el segmento nuevo correspondiente al semiplano. Como tal procedimiento se puede implementar en tiempo lineal, podemos simplemente empezar con una caja acotante grande y recortarla con cada uno de los semiplanos, obteniendo una complejidad temporal total de $O(N^2)$.

Este método es un gran paso en la dirección correcta, pero se siente derrochador tener que iterar sobre $O(N)$ semiplanos en cada paso. Veremos a continuación que, haciendo algunas observaciones inteligentes, las ideas detrás de este enfoque incremental se pueden reciclar para crear un algoritmo $O(N \log N)$.

Algoritmo Sort-and-Incremental - $O(N \log N)$ {data-toc-label=“Algoritmo Sort-and-Incremental - O(N log N)”}

La primera fuente documentada adecuadamente de este algoritmo que pudimos encontrar fue la tesis de Zeyuan Zhu para el Chinese Team Selecting Contest titulada New Algorithm for Half-plane Intersection and its Practical Value , del año 2006. El enfoque que describiremos a continuación se basa en este mismo algoritmo, pero en lugar de computar dos intersecciones separadas para las mitades inferior y superior de las intersecciones, la construiremos toda de una vez en una sola pasada con un deque (cola de dos extremos).

El algoritmo mismo, como el nombre puede adelantar, aprovecha el hecho de que la región resultante de la intersección de semiplanos es convexa, y por tanto consistirá de algunos segmentos de semiplanos en orden ordenados por sus ángulos. Esto lleva a una observación crucial: si intersectamos de forma incremental los semiplanos en su orden ordenado por ángulo (como aparecerían en la forma final resultante de la intersección) y los guardamos en una cola de dos extremos, entonces solo necesitaremos quitar semiplanos del frente y del fondo del deque.

Para visualizar mejor este hecho, supongamos que estamos realizando el enfoque incremental descrito previamente sobre un conjunto de semiplanos que está ordenado por ángulo (en este caso, asumiremos que están ordenados de $-\pi$ a $\pi$), y supongamos que estamos a punto de empezar algún $k$-ésimo paso arbitrario. Esto significa que ya construimos la intersección de los primeros $k-1$ semiplanos. Ahora, como los semiplanos están ordenados por ángulo, sea cual sea el $k$-ésimo semiplano, podemos estar seguros de que formará un giro convexo con el $(K-1)$-ésimo semiplano. Por esa razón, pueden ocurrir algunas cosas:

1. Algunos (posiblemente ninguno) de los semiplanos del fondo de la intersección pueden volverse redundantes. En este caso, hay que hacer pop de estos semiplanos ahora inútiles del fondo del deque.
2. Algunos (posiblemente ninguno) de los semiplanos del frente pueden volverse redundantes. Análogo al caso 1, simplemente los hacemos pop del frente del deque.
3. La intersección puede volverse vacía (después de manejar los casos 1 y/o 2). En este caso, simplemente reportamos que la intersección está vacía y terminamos el algoritmo.

Decimos que un semiplano es “redundante” si no contribuye nada a la intersección. Tal semiplano se podría quitar y la intersección resultante no cambiaría en absoluto.

Aquí hay un pequeño ejemplo con una ilustración:

Sea $H = { A, B, C, D, E }$ el conjunto de semiplanos presentes actualmente en la intersección. Además, sea $P = { p, q, r, s }$ el conjunto de puntos de intersección de semiplanos adyacentes en H. Ahora, supongamos que queremos intersectarlo con el semiplano $F$, como se ve en la ilustración de abajo:

Nótese que el semiplano $F$ hace que $A$ y $E$ sean redundantes en la intersección. Así que quitamos tanto $A$ como $E$ del frente y del fondo de la intersección, respectivamente, y agregamos $F$ al final. Y finalmente obtenemos la intersección nueva $H = { B, C, D, F}$ con $P = { q, r, t, u }$.

Con todo esto en mente, tenemos casi todo lo que necesitamos para implementar de hecho el algoritmo, pero todavía hay que hablar de algunos casos especiales. Al principio del artículo dijimos que agregaríamos una caja acotante para cuidar de los casos en los que la intersección podría ser no acotada, así que el único caso delicado que realmente hay que manejar es el de semiplanos paralelos. Podemos tener dos subcasos: dos semiplanos pueden ser paralelos con la misma dirección o con dirección opuesta. La razón por la que este caso hay que manejarlo por separado es porque necesitaremos computar puntos de intersección de las rectas de los semiplanos para poder comprobar si un semiplano es redundante o no, y dos rectas paralelas no tienen punto de intersección, así que necesitamos una forma especial de lidiar con ellas.

Para el caso de semiplanos paralelos de orientación opuesta: Nótese que, como estamos agregando la caja acotante para lidiar con el caso no acotado, esto también lidia con el caso en el que tenemos dos semiplanos paralelos adyacentes con direcciones opuestas después de ordenar, ya que tendrá que haber al menos uno de los semiplanos de la caja acotante entre estos dos (recordemos que están ordenados por ángulo).

• Sin embargo, es posible que, después de quitar algunos semiplanos del fondo del deque, dos semiplanos paralelos de dirección opuesta terminen juntos. Este caso solo ocurre, específicamente, cuando estos dos semiplanos forman una intersección vacía, ya que este último semiplano hará que se quite todo del deque. Para evitar este problema, tenemos que comprobar manualmente si hay semiplanos paralelos, y si tienen dirección opuesta, simplemente detenemos el algoritmo al instante y devolvemos una intersección vacía.

Así el único caso que realmente hay que manejar es tener múltiples semiplanos con el mismo ángulo, y resulta que este caso es bastante fácil de manejar: solo hay que conservar el semiplano más a la izquierda y borrar el resto, ya que serán completamente redundantes de todos modos. Para resumir, el algoritmo completo se verá aproximadamente así:

1. Empezamos ordenando el conjunto de semiplanos por ángulo, lo que toma tiempo $O(N \log N)$.
2. Iteraremos sobre el conjunto de semiplanos, y para cada uno, realizaremos el procedimiento incremental, haciendo pop del frente y del fondo de la cola de dos extremos según sea necesario. Esto tomará tiempo lineal en total, ya que cada semiplano solo se puede agregar o quitar una vez.
3. Al final, el polígono convexo resultante de la intersección se puede obtener simplemente computando los puntos de intersección de semiplanos adyacentes en el deque al final del procedimiento. Esto tomará tiempo lineal también. También es posible guardar tales puntos durante el paso 2 y saltarse este paso por completo, pero creemos que es ligeramente más fácil (en términos de implementación) computarlos sobre la marcha.

En total, hemos alcanzado una complejidad temporal de $O(N \log N)$. Como el ordenamiento es claramente el cuello de botella, el algoritmo se puede hacer correr en tiempo lineal en el caso especial en el que nos dan semiplanos ya ordenados de antemano por sus ángulos (un ejemplo de tal caso sería obtener los semiplanos que definen un polígono convexo).

Implementación directa

Aquí hay una implementación de ejemplo, directa, del algoritmo, con comentarios que explican la mayoría de las partes:

Structs simples de punto/vector y semiplano:

// Redefine epsilon and infinity as necessary. Be mindful of precision errors.
const long double eps = 1e-9, inf = 1e9; 
 
// Basic point/vector struct.
struct Point { 
 
    long double x, y;
    explicit Point(long double x = 0, long double y = 0) : x(x), y(y) {}
 
    // Addition, subtraction, multiply by constant, dot product, cross product.
 
    friend Point operator + (const Point& p, const Point& q) {
        return Point(p.x + q.x, p.y + q.y); 
    }
 
    friend Point operator - (const Point& p, const Point& q) { 
        return Point(p.x - q.x, p.y - q.y); 
    }
 
    friend Point operator * (const Point& p, const long double& k) { 
        return Point(p.x * k, p.y * k); 
    } 
    
    friend long double dot(const Point& p, const Point& q) {
    	return p.x * q.x + p.y * q.y;
    }
 
    friend long double cross(const Point& p, const Point& q) { 
        return p.x * q.y - p.y * q.x; 
    }
};
 
// Basic half-plane struct.
struct Halfplane { 
 
    // 'p' is a passing point of the line and 'pq' is the direction vector of the line.
    Point p, pq; 
    long double angle;
 
    Halfplane() {}
    Halfplane(const Point& a, const Point& b) : p(a), pq(b - a) {
        angle = atan2l(pq.y, pq.x);    
    }
 
    // Check if point 'r' is outside this half-plane. 
    // Every half-plane allows the region to the LEFT of its line.
    bool out(const Point& r) { 
        return cross(pq, r - p) < -eps; 
    }
 
    // Comparator for sorting. 
    bool operator < (const Halfplane& e) const { 
        return angle < e.angle;
    } 
 
    // Intersection point of the lines of two half-planes. It is assumed they're never parallel.
    friend Point inter(const Halfplane& s, const Halfplane& t) {
        long double alpha = cross((t.p - s.p), t.pq) / cross(s.pq, t.pq);
        return s.p + (s.pq * alpha);
    }
};

Algoritmo:

// Actual algorithm
vector<Point> hp_intersect(vector<Halfplane>& H) { 
 
    Point box[4] = {  // Bounding box in CCW order
        Point(inf, inf), 
        Point(-inf, inf), 
        Point(-inf, -inf), 
        Point(inf, -inf) 
    };
 
    for(int i = 0; i<4; i++) { // Add bounding box half-planes.
        Halfplane aux(box[i], box[(i+1) % 4]);
        H.push_back(aux);
    }
 
    // Sort by angle and start algorithm
    sort(H.begin(), H.end());
    deque<Halfplane> dq;
    int len = 0;
    for(int i = 0; i < int(H.size()); i++) {
 
        // Remove from the back of the deque while last half-plane is redundant
        while (len > 1 && H[i].out(inter(dq[len-1], dq[len-2]))) {
            dq.pop_back();
            --len;
        }
 
        // Remove from the front of the deque while first half-plane is redundant
        while (len > 1 && H[i].out(inter(dq[0], dq[1]))) {
            dq.pop_front();
            --len;
        }
        
        // Special case check: Parallel half-planes
        if (len > 0 && fabsl(cross(H[i].pq, dq[len-1].pq)) < eps) {
        	// Opposite parallel half-planes that ended up checked against each other.
        	if (dot(H[i].pq, dq[len-1].pq) < 0.0)
        		return vector<Point>();
        	
        	// Same direction half-plane: keep only the leftmost half-plane.
        	if (H[i].out(dq[len-1].p)) {
        		dq.pop_back();
        		--len;
        	}
        	else continue;
        }
        
        // Add new half-plane
        dq.push_back(H[i]);
        ++len;
    }
 
    // Final cleanup: Check half-planes at the front against the back and vice-versa
    while (len > 2 && dq[0].out(inter(dq[len-1], dq[len-2]))) {
        dq.pop_back();
        --len;
    }
 
    while (len > 2 && dq[len-1].out(inter(dq[0], dq[1]))) {
        dq.pop_front();
        --len;
    }
 
    // Report empty intersection if necessary
    if (len < 3) return vector<Point>();
 
    // Reconstruct the convex polygon from the remaining half-planes.
    vector<Point> ret(len);
    for(int i = 0; i+1 < len; i++) {
        ret[i] = inter(dq[i], dq[i+1]);
    }
    ret.back() = inter(dq[len-1], dq[0]);
    return ret;
}

Discusión de la implementación

Una cosa especial a notar es que, en caso de que haya múltiples semiplanos que se intersectan en el mismo punto, entonces este algoritmo podría devolver puntos adyacentes repetidos en el polígono final. Sin embargo, esto no debería tener ningún impacto en juzgar correctamente si la intersección está vacía o no, y tampoco afecta el área del polígono en absoluto. Se puede querer quitar estos duplicados dependiendo de qué tareas hay que hacer después. Se puede hacer esto muy fácilmente con std::unique. Queremos conservar los puntos repetidos durante la ejecución del algoritmo para que las intersecciones con área igual a cero se puedan computar correctamente (por ejemplo, intersecciones que consisten de un solo punto, una recta o un segmento de recta). Animo al lector a probar algunos casos pequeños hechos a mano donde la intersección resulta en un solo punto o una recta.

Una cosa más de la que se debería hablar es qué hacer si nos dan semiplanos en forma de una restricción lineal (por ejemplo, $ax + by + c \leq 0$). En tal caso, hay dos opciones. Se puede o bien implementar el algoritmo con las modificaciones correspondientes para trabajar con tal representación (esencialmente crear un struct propio de semiplano; debería ser bastante directo si se está familiarizado con el truco de la envolvente convexa), o se pueden transformar las rectas a la representación que usamos en este artículo tomando cualesquiera 2 puntos de cada recta. En general, se recomienda trabajar con la representación que se da en el problema para evitar problemas adicionales de precisión.

Problemas, tareas y aplicaciones

Muchos problemas que se pueden resolver con intersección de semiplanos también se pueden resolver sin ella, pero con enfoques (usualmente) más complicados o poco comunes. En general, la intersección de semiplanos puede aparecer al lidiar con problemas relacionados con polígonos (principalmente convexos), visibilidad en el plano y programación lineal bidimensional. Aquí hay algunas tareas de ejemplo que se pueden resolver con esta técnica:

Intersección de polígonos convexos

Una de las aplicaciones clásicas de la intersección de semiplanos: Dados $N$ polígonos, computar la región que está incluida dentro de todos los polígonos.

Como la intersección de un conjunto de semiplanos es un polígono convexo, también podemos representar un polígono convexo como un conjunto de semiplanos (cada arista del polígono es un segmento de un semiplano). Generar estos semiplanos para cada polígono y computar la intersección de todo el conjunto. La complejidad temporal total es $O(S \log S)$, donde S es el número total de lados de todos los polígonos. El problema también se puede resolver teóricamente en $O(S \log N)$ fusionando los $N$ conjuntos de semiplanos usando un heap y luego ejecutando el algoritmo sin el paso de ordenamiento, pero tal solución tiene un factor constante mucho peor que el ordenamiento directo y solo da ganancias menores de velocidad para $N$ muy pequeño.

Visibilidad en el plano

Los problemas que requieren algo en la línea de “determinar si algunos segmentos de recta son visibles desde algún(os) punto(s) en el plano” usualmente se pueden formular como problemas de intersección de semiplanos. Tómese, por ejemplo, la siguiente tarea: Dado algún polígono simple (no necesariamente convexo), determinar si hay algún punto dentro del polígono tal que todo el borde del polígono se pueda observar desde ese punto. Esto también se conoce como encontrar el kernel de un polígono  y se puede resolver por simple intersección de semiplanos, tomando cada arista del polígono como un semiplano y luego computando su intersección.

Aquí hay un problema relacionado, más interesante, que presentó Artem Vasilyev en una de sus clases de la Brazilian ICPC Summer School : Dado un conjunto $p$ de puntos $p_1, p_2\ \dots \ p_n$ en el plano, determinar si hay algún punto $q$ en el que se pueda estar de pie tal que se puedan ver todos los puntos de $p$ de izquierda a derecha en orden creciente de su índice.

Tal problema se puede resolver notando que poder ver algún punto $p_i$ a la izquierda de $p_j$ es lo mismo que poder ver el lado derecho del segmento de recta de $p_i$ a $p_j$ (o de forma equivalente, poder ver el lado izquierdo del segmento de $p_j$ a $p_i$). Con eso en mente, podemos simplemente crear un semiplano para cada segmento de recta $p_i p_{i+1}$ (o $p_{i+1} p_i$ según la orientación que se elija) y comprobar si la intersección de todo el conjunto está vacía o no.

Intersección de semiplanos con búsqueda binaria

Otra aplicación común es utilizar la intersección de semiplanos como herramienta para validar el predicado de un procedimiento de búsqueda binaria. Aquí hay un ejemplo de tal problema, también presentado por Artem Vasilyev en la misma clase que se mencionó antes: Dado un polígono convexo $P$, encontrar la mayor circunferencia que se puede inscribir dentro de él.

En lugar de buscar algún tipo de solución en forma cerrada, fórmulas molestas o soluciones algorítmicas oscuras, intentemos en su lugar hacer búsqueda binaria sobre la respuesta. Nótese que, para algún $r$ fijo, un círculo de radio $r$ se puede inscribir dentro de $P$ solo si existe algún punto dentro de $P$ que tiene distancia mayor o igual que $r$ a todos los puntos del borde de $P$. Esta condición se puede validar “encogiendo” el polígono hacia adentro en una distancia $r$ y comprobando que el polígono permanece no degenerado (o es un punto/segmento en sí mismo). Tal procedimiento se puede simular tomando los semiplanos de los lados del polígono en orden antihorario, trasladando cada uno de ellos una distancia $r$ en la dirección de la región que permiten (es decir, ortogonal al vector director del semiplano), y comprobando si la intersección no está vacía.

Claramente, si podemos inscribir un círculo de radio $r$, también podemos inscribir cualquier otro círculo de radio menor que $r$. Así que podemos realizar una búsqueda binaria sobre el radio $r$ y validar cada paso usando intersección de semiplanos. También, nótese que los semiplanos de un polígono convexo ya están ordenados por ángulo, así que el paso de ordenamiento se puede saltar en el algoritmo. Así obtenemos una complejidad temporal total de $O(NK)$, donde $N$ es el número de vértices del polígono y $K$ es el número de iteraciones de la búsqueda binaria (el valor real dependerá del rango de respuestas posibles y de la precisión deseada).

Programación lineal bidimensional

Una aplicación más de la intersección de semiplanos es la programación lineal en dos variables. Todas las restricciones lineales para dos variables se pueden expresar en la forma $Ax + By + C \leq 0$ (el comparador de desigualdad puede variar). Claramente, estos son solo semiplanos, así que comprobar si existe una solución factible para un conjunto de restricciones lineales se puede hacer con intersección de semiplanos. Además, para un conjunto dado de restricciones lineales, es posible computar la región de soluciones factibles (es decir, la intersección de los semiplanos) y luego responder múltiples consultas de maximizar/minimizar alguna función lineal $f(x, y)$ sujeta a las restricciones en $O(\log N)$ por consulta usando búsqueda binaria (muy similar al truco de la envolvente convexa).

Vale la pena mencionar que también existe un algoritmo aleatorizado bastante simple que puede comprobar si un conjunto de restricciones lineales tiene una solución factible o no, y maximizar/minimizar alguna función lineal sujeta a las restricciones dadas. Este algoritmo aleatorizado también lo explicó bien Artem Vasilyev en la clase mencionada antes. Aquí hay algunos recursos adicionales sobre él, por si el lector está interesado: CG - Lecture 4, parts 4 and 5  y el blog de Petr Mitrichev (que incluye la solución al problema más difícil de la lista de problemas de práctica de abajo) .

Problemas de práctica

Problemas clásicos, aplicación directa

• Codechef - Animesh decides to settle down 
• POJ - How I mathematician Wonder What You Are! 
• POJ - Rotating Scoreboard 
• POJ - Video Surveillance 
• POJ - Art Gallery 
• POJ - Uyuw’s Concert 

Problemas más difíciles

• POJ - Most Distant Point from the Sea - Medium 
• Baekjoon - Jeju’s Island - Same as above but seemingly stronger test cases 
• POJ - Feng Shui - Medium 
• POJ - Triathlon - Medium/hard 
• DMOJ - Arrow - Medium/hard 
• POJ - Jungle Outpost - Hard 
• Codeforces - Jungle Outpost (alternative link, problem J) - Hard 
• Yandex - Asymmetry Value (need virtual contest to see, problem F) - Very Hard 

Problemas adicionales

• 40th Petrozavodsk Programming Camp, Winter 2021 - Day 1: Jagiellonian U Contest, Grand Prix of Krakow - Problem B: (Almost) Fair Cake-Cutting. Al momento de escribir el artículo, este problema era privado y solo accesible para los participantes del Programming Camp.

Referencias, bibliografía y otras fuentes

Fuentes principales

• New Algorithm for Half-plane Intersection and its Practical Value.  Artículo original del algoritmo.
• Artem Vasilyev’s Brazilian ICPC Summer School 2020 lecture.  Clase excelente sobre intersección de semiplanos. También cubre otros temas de geometría.

Buenos blogs (chino)

• Fundamentals of Computational Geometry - Intersection of Half-planes. 
• Detailed introduction to the half-plane intersection algorithm. 
• Summary of Half-plane intersection problems. 
• Sorting incremental method of half-plane intersection. 

Algoritmo aleatorizado

• Linear Programming and Half-Plane intersection - Parts 4 and 5. 
• Petr Mitrichev’s Blog: A half-plane week.