Vertical decomposition¶
Overview¶
Vertical decomposition is a powerful technique used in various geometry problems. The general idea is to cut the plane into several vertical stripes with some "good" properties and solve the problem for these stripes independently. We will illustrate the idea on some examples.
Area of the union of triangles¶
Suppose that there are
Optimization 1¶
Firstly we will reduce the runtime to
-
In this case
-
This case can be reduced to the previous one by splitting
-
This case means that the parts of
Here is a graphic representation of the three cases.

Finally we should remark on processing all the additions of
Optimization 2¶
Note that if we apply the previous optimization, we no longer have to find all stripes explicitly. This reduces the memory consumption to
Intersection of convex polygons¶
Another usage of vertical decomposition is to compute the intersection of two convex polygons in linear time. Suppose the plane is split into vertical stripes by vertical lines passing through each vertex of each polygon. Then if we consider one of the input polygons and some stripe, their intersection is either a trapezoid, a triangle or a point. Therefore we can simply intersect these shapes for each vertical stripe and merge these intersections into a single polygon.
Implementation¶
Below is the code that calculates area of the union of a set of triangles in
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;
}