Descomposición vertical

Panorama

La descomposición vertical es una técnica potente que se usa en varios problemas de geometría. La idea general es cortar el plano en varias franjas verticales con ciertas propiedades «buenas» y resolver el problema para estas franjas de forma independiente. Ilustraremos la idea con algunos ejemplos.

Área de la unión de triángulos

Supongamos que hay $n$ triángulos en un plano y hay que encontrar el área de su unión. El problema sería fácil si los triángulos no se intersectaran, así que eliminemos estas intersecciones dividiendo el plano en franjas verticales al trazar rectas verticales por todos los vértices y por todos los puntos de intersección de lados de triángulos distintos. Puede haber $O(n^2)$ de tales rectas, así que obtuvimos $O(n^2)$ franjas. Ahora consideremos alguna franja vertical. Cada segmento no vertical o bien la cruza de izquierda a derecha, o no la cruza en absoluto. Además, no hay dos segmentos que se intersecten estrictamente dentro de la franja. Eso significa que la parte de la unión de triángulos que yace dentro de esta franja está compuesta de trapecios disjuntos con bases sobre los lados de la franja. Esta propiedad nos permite calcular el área dentro de cada franja con el siguiente algoritmo de línea de barrido. Cada segmento que cruza la franja es superior o inferior, según si el interior del triángulo correspondiente está por encima o por debajo del segmento. Podemos visualizar cada segmento superior como un paréntesis de apertura y cada segmento inferior como un paréntesis de cierre, y descomponer la franja en trapecios descomponiendo la secuencia de paréntesis en secuencias de paréntesis correctas más pequeñas. Este algoritmo requiere tiempo $O(n^3\log n)$ y memoria $O(n^2)$ .

Optimización 1

En primer lugar reduciremos el tiempo de ejecución a $O(n^2\log n)$ . En lugar de generar trapecios para cada franja, fijemos algún lado de un triángulo (el segmento $s = (s_0, s_1)$ ) y encontremos el conjunto de franjas donde este segmento es un lado de algún trapecio. Nótese que en este caso solo hay que encontrar las franjas donde el balance de paréntesis por debajo (o por encima, en el caso de un segmento inferior) de $s$ es cero. Eso significa que, en lugar de ejecutar una línea de barrido vertical para cada franja, podemos ejecutar una línea de barrido horizontal para todas las partes de otros segmentos que afectan el balance de paréntesis respecto de $s$ . Por simplicidad mostraremos cómo hacerlo para un segmento superior; el algoritmo para segmentos inferiores es similar. Consideremos algún otro segmento no vertical $t = (t_0, t_1)$ y encontremos la intersección $[x_1, x_2]$ de las proyecciones de $s$ y $t$ sobre $Ox$ . Si esta intersección es vacía o consiste en un solo punto, $t$ se puede descartar, ya que $s$ y $t$ no intersectan el interior de la misma franja. En caso contrario, consideremos la intersección $I$ de $s$ y $t$ . Hay tres casos.

$I = \varnothing$

En este caso $t$ está o bien por encima o por debajo de $s$ en $[x_1, x_2]$ . Si $t$ está por encima, no afecta si $s$ es o no un lado de algún trapecio. Si $t$ está por debajo de $s$ , debemos sumar $1$ o $-1$ al balance de las secuencias de paréntesis para todas las franjas en $[x_1, x_2]$ , según si $t$ es superior o inferior.
$I$ consiste en un solo punto $p$

Este caso se puede reducir al anterior partiendo $[x_1, x_2]$ en $[x_1, p_x]$ y $[p_x, x_2]$ .
$I$ es algún segmento $l$

Este caso significa que las partes de $s$ y $t$ para $x\in[x_1, x_2]$ coinciden. Si $t$ es inferior, $s$ claramente no es un lado de un trapecio. En caso contrario, podría ocurrir que tanto $s$ como $t$ se puedan considerar como lado de algún trapecio. Para resolver esta ambigüedad, podemos decidir que solo el segmento con el menor índice se considere como lado (aquí suponemos que los lados de los triángulos están enumerados de alguna forma). Así, si $index(s) < index(t)$ , debemos ignorar este caso; en caso contrario debemos marcar que $s$ nunca puede ser un lado en $[x_1, x_2]$ (por ejemplo, agregando un evento correspondiente con balance $-2$ ).

Aquí hay una representación gráfica de los tres casos.

Finalmente debemos comentar cómo procesar todas las sumas de $1$ o $-1$ en todas las franjas de $[x_1, x_2]$ . Para cada suma de $w$ en $[x_1, x_2]$ podemos crear eventos $(x_1, w),\ (x_2, -w)$ y procesar todos estos eventos con una línea de barrido.

Optimización 2

Nótese que si aplicamos la optimización anterior, ya no hace falta encontrar todas las franjas de forma explícita. Esto reduce el consumo de memoria a $O(n)$ .

Intersección de polígonos convexos

Otro uso de la descomposición vertical es computar la intersección de dos polígonos convexos en tiempo lineal. Supongamos que el plano se parte en franjas verticales por rectas verticales que pasan por cada vértice de cada polígono. Entonces, si consideramos uno de los polígonos de entrada y alguna franja, su intersección es un trapecio, un triángulo o un punto. Por tanto, podemos simplemente intersectar estas formas para cada franja vertical y fusionar estas intersecciones en un solo polígono.

Implementación

A continuación está el código que calcula el área de la unión de un conjunto de triángulos en tiempo $O(n^2\log n)$ y memoria $O(n)$ .


typedef double dbl;

const dbl eps = 1e-9;
 
inline bool eq(dbl x, dbl y){
    return fabs(x - y) < eps;
}
 
inline bool lt(dbl x, dbl y){
    return x < y - eps;
}
 
inline bool gt(dbl x, dbl y){
    return x > y + eps;
}
 
inline bool le(dbl x, dbl y){
    return x < y + eps;
}
 
inline bool ge(dbl x, dbl y){
    return x > y - eps;
}
 
struct pt{
    dbl x, y;
    inline pt operator - (const pt & p)const{
        return pt{x - p.x, y - p.y};
    }
    inline pt operator + (const pt & p)const{
        return pt{x + p.x, y + p.y};
    }
    inline pt operator * (dbl a)const{
        return pt{x * a, y * a};
    }
    inline dbl cross(const pt & p)const{
        return x * p.y - y * p.x;
    }
    inline dbl dot(const pt & p)const{
        return x * p.x + y * p.y;
    }
    inline bool operator == (const pt & p)const{
        return eq(x, p.x) && eq(y, p.y);
    }
};
 
struct Line{
    pt p[2];
    Line(){}
    Line(pt a, pt b):p{a, b}{}
    pt vec()const{
        return p[1] - p[0];
    }
    pt& operator [](size_t i){
        return p[i];
    }
};

inline bool lexComp(const pt & l, const pt & r){
	if(fabs(l.x - r.x) > eps){
		return l.x < r.x;
	}
	else return l.y < r.y;
}
 
vector<pt> interSegSeg(Line l1, Line l2){
    if(eq(l1.vec().cross(l2.vec()), 0)){
        if(!eq(l1.vec().cross(l2[0] - l1[0]), 0))
            return {};
        if(!lexComp(l1[0], l1[1]))
            swap(l1[0], l1[1]);
        if(!lexComp(l2[0], l2[1]))
            swap(l2[0], l2[1]);
        pt l = lexComp(l1[0], l2[0]) ? l2[0] : l1[0];
        pt r = lexComp(l1[1], l2[1]) ? l1[1] : l2[1];
        if(l == r)
            return {l};
        else return lexComp(l, r) ? vector<pt>{l, r} : vector<pt>();
    }
    else{
        dbl s = (l2[0] - l1[0]).cross(l2.vec()) / l1.vec().cross(l2.vec());
        pt inter = l1[0] + l1.vec() * s;
        if(ge(s, 0) && le(s, 1) && le((l2[0] - inter).dot(l2[1] - inter), 0))
            return {inter};
        else
            return {};
    }
}
inline char get_segtype(Line segment, pt other_point){
    if(eq(segment[0].x, segment[1].x))
        return 0;
    if(!lexComp(segment[0], segment[1]))
        swap(segment[0], segment[1]);
    return (segment[1] - segment[0]).cross(other_point - segment[0]) > 0 ? 1 : -1;
}
 
dbl union_area(vector<tuple<pt, pt, pt> > triangles){
    vector<Line> segments(3 * triangles.size());
    vector<char> segtype(segments.size());
    for(size_t i = 0; i < triangles.size(); i++){
        pt a, b, c;
        tie(a, b, c) = triangles[i];
        segments[3 * i] = lexComp(a, b) ? Line(a, b) : Line(b, a);
        segtype[3 * i] = get_segtype(segments[3 * i], c);
        segments[3 * i + 1] = lexComp(b, c) ? Line(b, c) : Line(c, b);
        segtype[3 * i + 1] = get_segtype(segments[3 * i + 1], a);
        segments[3 * i + 2] = lexComp(c, a) ? Line(c, a) : Line(a, c);
        segtype[3 * i + 2] = get_segtype(segments[3 * i + 2], b);
    }
    vector<dbl> k(segments.size()), b(segments.size());
    for(size_t i = 0; i < segments.size(); i++){
        if(segtype[i]){
            k[i] = (segments[i][1].y - segments[i][0].y) / (segments[i][1].x - segments[i][0].x);
            b[i] = segments[i][0].y - k[i] * segments[i][0].x;
        }
    }
    dbl ans = 0;
    for(size_t i = 0; i < segments.size(); i++){
        if(!segtype[i])
            continue;
        dbl l = segments[i][0].x, r = segments[i][1].x;
        vector<pair<dbl, int> > evts;
        for(size_t j = 0; j < segments.size(); j++){
            if(!segtype[j] || i == j)
                continue;
            dbl l1 = segments[j][0].x, r1 = segments[j][1].x;
            if(ge(l1, r) || ge(l, r1))
                continue;
            dbl common_l = max(l, l1), common_r = min(r, r1);
            auto pts = interSegSeg(segments[i], segments[j]);
            if(pts.empty()){
                dbl yl1 = k[j] * common_l + b[j];
                dbl yl = k[i] * common_l + b[i];
                if(lt(yl1, yl) == (segtype[i] == 1)){
                    int evt_type = -segtype[i] * segtype[j];
                    evts.emplace_back(common_l, evt_type);
                    evts.emplace_back(common_r, -evt_type);
                }
            }
            else if(pts.size() == 1u){
                dbl yl = k[i] * common_l + b[i], yl1 = k[j] * common_l + b[j];
                int evt_type = -segtype[i] * segtype[j];
                if(lt(yl1, yl) == (segtype[i] == 1)){
                    evts.emplace_back(common_l, evt_type);
                    evts.emplace_back(pts[0].x, -evt_type);
                }
                yl = k[i] * common_r + b[i], yl1 = k[j] * common_r + b[j];
                if(lt(yl1, yl) == (segtype[i] == 1)){
                    evts.emplace_back(pts[0].x, evt_type);
                    evts.emplace_back(common_r, -evt_type);
                }
            }
            else{
                if(segtype[j] != segtype[i] || j > i){
                    evts.emplace_back(common_l, -2);
                    evts.emplace_back(common_r, 2);
                }
            }
        }
        evts.emplace_back(l, 0);
        sort(evts.begin(), evts.end());
        size_t j = 0;
        int balance = 0;
        while(j < evts.size()){
            size_t ptr = j;
            while(ptr < evts.size() && eq(evts[j].first, evts[ptr].first)){
                balance += evts[ptr].second;
                ++ptr;
            }
            if(!balance && !eq(evts[j].first, r)){
                dbl next_x = ptr == evts.size() ? r : evts[ptr].first;
                ans -= segtype[i] * (k[i] * (next_x + evts[j].first) + 2 * b[i]) * (next_x - evts[j].first);
            }
            j = ptr;
        }
    }
    return ans/2;
}