Círculo mínimo envolvente

Considérese el siguiente problema:

Aquí, por círculo mínimo envolvente (MEC, minimum enclosing circle) entendemos un círculo de radio lo menor posible que contiene los $n$ puntos, dentro del círculo o sobre su frontera. Este problema tiene una solución aleatorizada simple que, a primera vista, parece que correría en $O(n^3)$, pero en realidad funciona en tiempo esperado $O(n)$.

Para entender mejor el razonamiento de abajo, debemos notar de inmediato que la solución del problema es única:

Algoritmo de Welzl

Por brevedad, denotemos $\operatorname{mec}(p_1,\dots,p_n)$ al MEC de ${p_1,\dots,p_n}$, y sea $P_i = {p_1,\dots,p_i}$.

El algoritmo, propuesto  inicialmente por Welzl en 1991, procede como sigue:

1. Aplicar una permutación aleatoria a la secuencia de entrada de puntos.
2. Mantener el candidato actual a ser el MEC $C$, empezando con $C = \operatorname{mec}(p_1, p_2)$.
3. Iterar sobre $i=3..n$ y comprobar si $p_i \in C$.
   1. Si $p_i \in C$ significa que $C$ es el MEC de $P_i$.
   2. En caso contrario, asignar $C = \operatorname{mec}(p_i, p_1)$ e iterar sobre $j=2..i$ y comprobar si $p_j \in C$.
      1. Si $p_j \in C$, entonces $C$ es el MEC de $P_j$ entre los círculos que pasan por $p_i$.
      2. En caso contrario, asignar $C=\operatorname{mec}(p_i, p_j)$ e iterar sobre $k=1..j$ y comprobar si $p_k \in C$.
         1. Si $p_k \in C$, entonces $C$ es el MEC de $P_k$ entre los círculos que pasan por $p_i$ y $p_j$.
         2. En caso contrario, $C=\operatorname{mec}(p_i,p_j,p_k)$ es el MEC de $P_k$ entre los círculos que pasan por $p_i$ y $p_j$.

Se puede ver que cada nivel de anidamiento aquí tiene un invariante que mantener (que $C$ es el MEC entre los círculos que además pasan por $0$, $1$ o $2$ puntos dados adicionalmente), y cada vez que se cierra el bucle interno, su invariante se vuelve equivalente al invariante de la iteración actual de su bucle padre. Esto, a su vez, asegura la corrección del algoritmo en su conjunto.

Omitiendo algunos detalles técnicos, por ahora, el algoritmo completo se puede implementar en C++ como sigue:

struct point {...};
 
// Se representa por 2 o 3 puntos sobre su circunferencia
struct mec {...};
 
bool inside(mec const& C, point p) {
    return ...;
}
 
// Elegir algún generador de aleatoriedad bueno para el shuffle
mt19937_64 gen(...);
mec enclosing_circle(vector<point> &p) {
    int n = p.size();
    ranges::shuffle(p, gen);
    auto C = mec{p[0], p[1]};
    for(int i = 0; i < n; i++) {
        if(!inside(C, p[i])) {
            C = mec{p[i], p[0]};
            for(int j = 0; j < i; j++) {
                if(!inside(C, p[j])) {
                    C = mec{p[i], p[j]};
                    for(int k = 0; k < j; k++) {
                        if(!inside(C, p[k])) {
                            C = mec{p[i], p[j], p[k]};
                        }
                    }
                }
            }
        }
    }
    return C;
}

Ahora, es de esperar que comprobar que un punto $p_i$ está dentro del MEC de $2$ o $3$ puntos se pueda hacer en $O(1)$ (lo discutiremos más adelante). Pero incluso entonces, el algoritmo de arriba parece que tomaría $O(n^3)$ en el peor caso solo por todos los bucles anidados. Entonces, ¿cómo es que afirmamos el tiempo esperado lineal? ¡Veámoslo!

Análisis de complejidad

Para el bucle más interno (sobre $k$), claramente su tiempo esperado es $O(j)$ operaciones. ¿Y el bucle sobre $j$?

Solo dispara el siguiente bucle si $p_j$ está en la frontera del MEC de $P_j$ que además pasa por el punto $i$, y quitar $p_j$ encogería aún más el círculo. De todos los puntos de $P_j$ solo puede haber a lo sumo $2$ puntos con esa propiedad, porque si hay más de $2$ puntos de $P_j$ en la frontera, significa que después de quitar cualquiera de ellos aún habrá al menos $3$ puntos en la frontera, suficientes para definir el círculo de forma única.

En otras palabras, después de la mezcla aleatoria inicial, hay a lo sumo probabilidad $\frac{2}{j}$ de que obtengamos uno de los a lo sumo dos puntos desafortunados como $p_j$. Sumando sobre todos los $j$ de $1$ a $i$, obtenemos el tiempo esperado de

$$\sum\limits_{j=1}^i \frac{2}{j} \cdot O(j) = O(i).$$

De exactamente la misma manera ahora también podemos demostrar que el bucle más externo tiene tiempo esperado $O(n)$.

Comprobar que un punto está en el MEC de 2 o 3 puntos

Ahora veamos el detalle de implementación de point y mec. En este problema resulta particularmente útil usar std::complex  como clase para los puntos:

using ftype = int64_t;
using point = complex<ftype>;

Como recordatorio, un número complejo es un número del tipo $x+yi$, donde $i^2=-1$ y $x, y \in \mathbb R$. En C++, tal número complejo se representa por un punto bidimensional $(x, y)$. Los números complejos ya implementan operaciones lineales básicas componente a componente (suma, multiplicación por un número real), pero también su multiplicación y división tienen cierto significado geométrico.

Sin entrar en demasiado detalle, notaremos la propiedad más importante para esta tarea en particular: multiplicar dos números complejos suma sus ángulos polares (contados desde $Ox$ en sentido antihorario), y tomar el conjugado (es decir, cambiar $z=x+yi$ por $\overline{z} = x-yi$) multiplica el ángulo polar por $-1$. Esto nos permite formular algunos criterios muy simples para si un punto $z$ está dentro del MEC de $2$ o $3$ puntos específicos.

MEC de 2 puntos

Para $2$ puntos $a$ y $b$, su MEC es simplemente el círculo centrado en $\frac{a+b}{2}$ con radio $\frac{|a-b|}{2}$, en otras palabras el círculo que tiene $ab$ como diámetro. Para comprobar si $z$ está dentro de este círculo simplemente hay que comprobar que el ángulo entre $za$ y $zb$ no es agudo.

Los ángulos interiores son obtusos, los ángulos exteriores son agudos y los ángulos sobre la circunferencia son rectos

De forma equivalente, hay que comprobar que

$$I_0=(b-z)\overline{(a-z)}$$

no tiene una coordenada real positiva (correspondiente a puntos que tienen un ángulo polar entre $-90^\circ$ y $90^\circ$).

MEC de 3 puntos

Añadir $z$ al triángulo $abc$ lo convertirá en un cuadrilátero. Considérese la siguiente expresión:

$$\angle azb + \angle bca$$

En un cuadrilátero cíclico , si $c$ y $z$ están del mismo lado de $ab$, entonces los ángulos son iguales, y sumarán $0^\circ$ cuando se sumen con signo (es decir, positivos si son antihorarios y negativos si son horarios). Correspondientemente, si $c$ y $z$ están en lados opuestos, los ángulos sumarán $180^\circ$.

Los ángulos inscritos adyacentes son iguales; los ángulos opuestos se complementan a 180 grados

En términos de números complejos, podemos notar que $\angle azb$ es el ángulo polar de $(b-z)\overline{(a-z)}$ y $\angle bca$ es el ángulo polar de $(a-c)\overline{(b-c)}$. Así, podemos concluir que $\angle azb + \angle bca$ es el ángulo polar de

$$I_1 = (b-z) \overline{(a-z)} (a-c) \overline{(b-c)}$$

Si el ángulo es $0^\circ$ o $180^\circ$, significa que la parte imaginaria de $I_1$ es $0$; en caso contrario podemos deducir si $z$ está dentro o fuera del círculo envolvente de $abc$ comprobando el signo de la parte imaginaria de $I_1$. La parte imaginaria positiva corresponde a ángulos positivos, y la parte imaginaria negativa corresponde a ángulos negativos.

Pero ¿cuál de ellos significa que $z$ está dentro o fuera del círculo? Como ya notamos, tener $z$ dentro del círculo generalmente aumenta la magnitud de $\angle azb$, mientras que tenerlo fuera del círculo la disminuye. Así, tenemos los siguientes 4 casos:

1. $\angle bca > 0^\circ$, $c$ del mismo lado de $ab$ que $z$. Entonces, $\angle azb < 0^\circ$, y $\angle azb + \angle bca < 0^\circ$ para puntos dentro del círculo.
2. $\angle bca < 0^\circ$, $c$ del mismo lado de $ab$ que $z$. Entonces, $\angle azb > 0^\circ$, y $\angle azb + \angle bca > 0^\circ$ para puntos dentro del círculo.
3. $\angle bca > 0^\circ$, $c$ del lado opuesto de $ab$ respecto de $z$. Entonces, $\angle azb > 0^\circ$ y $\angle azb + \angle bca > 180^\circ$ para puntos dentro del círculo.
4. $\angle bca < 0^\circ$, $c$ del lado opuesto de $ab$ respecto de $z$. Entonces, $\angle azb < 0^\circ$ y $\angle azb + \angle bca < 180^\circ$ para puntos dentro del círculo.

En otras palabras, si $\angle bca$ es positivo, los puntos dentro del círculo tendrán $\angle azb + \angle bca < 0^\circ$; en caso contrario tendrán $\angle azb + \angle bca > 0^\circ$, asumiendo que normalizamos los ángulos entre $-180^\circ$ y $180^\circ$. Esto, a su vez, se puede comprobar por los signos de las partes imaginarias de $I_2=(a-c)\overline{(b-c)}$ y $I_1 = I_0 I_2$.

Nota: Como multiplicamos cuatro números complejos para obtener $I_1$, los coeficientes intermedios pueden ser tan grandes como $O(A^4)$, donde $A$ es la mayor magnitud de coordenada en la entrada. Del lado positivo, si la entrada es entera, ambas comprobaciones de arriba se pueden hacer enteramente en enteros.

Implementación

Ahora, para implementar realmente la comprobación, primero debemos decidir cómo representar el MEC. Como nuestros criterios trabajan con los puntos directamente, una forma natural y eficiente de hacerlo es decir que el MEC se representa directamente como un par o una terna de puntos que lo define:

using mec = variant<
    array<point, 2>,
    array<point, 3>
>;

Ahora, podemos usar std::visit para tratar ambos casos de forma eficiente de acuerdo con los criterios de arriba:

/* I < 0 si z está dentro de C,
   I > 0 si z está fuera de C,
   I = 0 si z está sobre la circunferencia de C */
ftype indicator(mec const& C, point z) {
    return visit([&](auto &&C) {
        point a = C[0], b = C[1];
        point I0 = (b - z) * conj(a - z);
        if constexpr (size(C) == 2) {
            return real(I0);
        } else {
            point c = C[2];
            point I2 = (a - c) * conj(b - c);
            point I1 = I0 * I2;
            return imag(I2) < 0 ? -imag(I1) : imag(I1);
        }
    }, C);
}
 
bool inside(mec const& C, point p) {
    return indicator(C, p) <= 0;
}
 

Ahora, por fin podemos asegurar que todo funciona enviando el problema al Library Checker: #308668 .

Problemas de práctica

• Library Checker - Minimum Enclosing Circle 
• BOI 2002 - Aliens