计算几何系列 —— 扫描线算法

大家好，欢迎来到IT知识分享网。

       本系列文章力求以简洁易懂的文字介绍计算几何中的基本概念，使读者快速入门，故不追求难度和深度，仅起到抛砖引玉的作用。     

       在我们的计算几何中，乃至基本的几何学中我们都需要解决一类问题，那就是几何图形的面积交并，周长交并问题，在计算几何中我们也常常遇到这样的问题，在这里我们介绍一种神奇的用来解决这类问题的算法，叫做扫描线（scan line algorithm）算法，是2018年公布的计算机科学技术名词，主要倾向于线段树+离散化这类思想，然后扫描求解过程有点像积分，那么我们开始讲解扫描线。

      我们先来看一个最基础的问题： 

      看到n的数据范围是：

      看到这个以后，我们是要去打ACM的话，肯定只考虑100%的数据，因此O(n^2)的算法直接被我们抛弃，这样子主要考虑O(n)或者O(nlogn)的算法，那么就是我们今天的算法:扫描线算法，时间复杂度是O(nlogn)的，接下来讲一下这个算法的运作过程。

扫描线：假设有一条扫描线从一个图形的下方扫向上方（或者左方扫到右方），那么通过分析扫描线被图形截得的线段就能获得所要的结果（相当于是每次算出每个单位长度上的微分，最后积起来这样的一个过程）。该过程可以用线段树进行加速。

       由于都是矩形，因此运用扫描线以后，面积的求法其实可以简化为 ∑截线段长度×扫过的高度，这也是扫描线算法最基础的应用。

  问题在于如何才能模拟扫描线从下向上扫过图形，并且快速计算出当前扫描线被截得的长度。

  现在假设，扫描线每次会在碰到横边的时候停下来，如图👇：

      简单来说，可对图形面积产生影响的元素，也就是这些横边左右端点的坐标。那么还有一个很严重的问题就是，如何加边减边的问题，就是在某个区间怎么判断这条边还在吗？这个很简单，我们只需要在一个矩形的上下两边加上不同的权值，按照扫描方向，假如方向是从下往上，那么下方的边权就是1，上方的边权就是-1，这样子就解决了边的存在问题了：       然后把所有的横边按照矩形y坐标升序排序。这样，对于每个矩形，扫描线总是会先碰到下边，然后再碰到上边。那么就能保证扫描线所截的长度永远非负了。这样操作以后，就可以和线段树扯上关系，我们只需要在输入的时候将每个矩形的横坐标剥离出来，放在x轴上进行离散化。（其实就是把所有点的横坐标存到X[]里，然后升序排个序，最后去重。）

       👆：在这个例子中，4个端点的纵投影总共把x轴切割成了5部分。取中间的3部分线段，建立一棵线段树，其每个端点维护一条线段（也就是一个区间）的信息：

1. 该线段被覆盖了多少次（被多少个矩形所覆盖）。
2. 该线段内被整个图形所截的长度是多少。

      显然，只要一条线段被覆盖，那么它肯定被图形所截。所以，整个问题就转化为了一个区间查询问题，即：每次将当前扫描线扫到的边对应的信息按照之前赋上的权值更新，然后再查询线段树根节点的信息，最后得到当前扫描线扫过的面积。这就可以用线段树来实现了（线段树：顾名思义，就是树上的结点为一条条线段~），接下来我们来简单看一下模拟过程：

       1.初始化，取点存线段，建立线段树：

       2：扫描第一条边：

         3：扫描第二条边：

          4：扫描第三条边：

          5：扫描第四条边：

      那么对于这样的一个基本情况就是这么简单的模拟，我们只需要建立一棵线段树来存储这些线段就OK了，但是问题来了，我们应该存储线段的那些信息？这些线段应该怎么存？（开区间 or 闭区间 or 半开半闭）

      所以为了保证线段之间的兼容性，我们一般这么考虑：考虑把线段树每个节点x对应的区间（tree[x].l,tree[x].r）不变，改变区间和横边的映射关系，具体为：节点x对应[X[tree[x].l],X[tree[x].r + 1]这条横边，如果不这样想可能会出现线段之间的端点重合问题，实在是细腻啊~，所以我们建树这样子建（保存的信息还是一样，提取的时候搞一下操作）：

        然后我们只需要在更新线段树的时候做操作~：

    OK了😄，算法思路和细节方面的问题我们都解决了，接下来就看一下这个算法的板子：

#include<bits/stdc++.h> #include<algorithm> #define ll long long using namespace std; const int MAXN = 1e6 + 10; int n,i,x1,x2,yy,y2; int X[MAXN << 1]; inline int read1(){char ch; while((ch = getchar()) < '0' || ch > '9'); int res = ch - 48; while((ch = getchar()) >= '0' && ch <= '9') res = res * 10 + ch - 48; return res;} inline long long read2(){char ch; while((ch = getchar()) < '0' || ch > '9'); long long res = ch - 48; while((ch = getchar()) >= '0' && ch <= '9') res = res * 10 + ch - 48; return res;} struct ScanLine{ long long l,r,h; int mark; bool operator<(const ScanLine &rhs) const {return h < rhs.h;} }line[MAXN << 1]; struct SegTree{int l,r,sum; long long len; }Tree[MAXN << 2]; void Build_Tree(int k,int l,int r) { Tree[k].l = l,Tree[k].r = r,Tree[k].sum = 0,Tree[k].len = 0; if(l == r) return; int mid = (l + r) >> 1; Build_Tree(k << 1,l,mid); Build_Tree(k << 1 | 1,mid + 1,r); return; } void Push_Up(int k) { int l = Tree[k].l,r = Tree[k].r; if(Tree[k].sum) Tree[k].len = X[r + 1] - X[l]; else Tree[k].len = Tree[k << 1].len + Tree[k << 1 | 1].len; } void Update_Tree(int k,long long L,long long R,int c) { int l = Tree[k].l,r = Tree[k].r; if(X[r + 1] <= L || X[l] >= R) return; if(X[l] >= L && X[r + 1] <= R) { Tree[k].sum = Tree[k].sum + c; Push_Up(k); return; } int mid = (l + r) >> 1; Update_Tree(k << 1,L,R,c); Update_Tree(k << 1 | 1,L,R,c); Push_Up(k); } int main() { n = read1(); memset(X,0,sizeof(X)); for(int i = 1;i <= n;i++) { {x1 = read2();yy = read2();x2 = read2();y2 = read2();} X[2 * i - 1] = x1,X[2 * i] = x2; line[2 * i - 1] = (ScanLine) {x1,x2,yy,1}; line[2 * i] = (ScanLine) {x1,x2,y2,-1}; } n = n << 1; sort(line + 1,line + n + 1); sort(X + 1,X + n + 1); int tot = unique(X + 1,X + n + 1) - X - 1; Build_Tree(1,1,tot - 1); long long ans = 0; for(int i = 1;i < n;i++) { Update_Tree(1,line[i].l,line[i].r,line[i].mark); ans = ans + Tree[1].len * (line[i + 1].h - line[i].h); } cout<<ans<<endl; return 0; }

好，我们趁热打铁（当然不是那个打铁~），再来看一道比较模板的扫描线题🤔：

     刚刚算过了矩形面积并之后，我们第一反应就可以想到截取线段，然后建立线段树维护，那么单单从下到上只能算出宽，不能算出高，所以很简单，我们再从左到右做一次同样的操作就行了~一样的道理，每次得到的贡献需要和上一次做差，表示这次新增或删除的贡献~✌

      由于每次扫描线扫描的过程和求面积并是相似的，所以这次我们不手动模拟了，那个之前忘记讲了，这个存放矩形其实也有一个结构体E（l,r,d,h):

      👆：感觉不是特别重要，但是这个计算几何的知识体系好歹完整一点😄

        OK,我们先来看看这一种做两次扫描线算法能不能再做简化，也就是只做一次扫描线，毕竟矩形的高与宽之间的联系还是非常大的，所以我们通过画图计算可以发现以下结论：

        1.纵边总长度=∑2×被截得的线段条数×扫过的高度

        2.横边总长度=∑∣上次截得的总长−现在截得的总长∣

          👆：大家可以自己模拟计算一下~

       那么事情看起来就简单很多了，我们只需要在原来线段树中存储的信息中再加上一个“线段的条数”，然后纵边和横边分开计算即可，代码如下😄：

#include<bits/stdc++.h> #include <algorithm> #define lson (x << 1) #define rson (x << 1 | 1) using namespace std; #define Temp template<typename T> Temp inline void read(T &x) { x = 0; T w = 1,ch = getchar(); while(! isdigit(ch) && ! ch == '-') ch = getchar(); if(ch == '-') w = -1,ch = getchar(); while(isdigit(ch)) x = (x << 3) + (x << 1) + (ch ^ '0'),ch = getchar(); x = x * w; } const int MAXN = 1e5 + 10; int n,X[MAXN << 1]; int x1,yy,x2,y2,pre = 0; struct ScanLine{ int l, r, h, mark; if(h == rhs.h) return mark > rhs.mark; return h < rhs.h; }line[MAXN << 1]; struct SegTree{ int l,r,sum,len,c; bool lc,rc; }tree[MAXN << 2]; void Build_Tree(int x, int l, int r) { tree[x].l = l; tree[x].r = r; tree[x].lc = tree[x].rc = false; tree[x].sum = tree[x].len = 0; tree[x].c = 0; if(l == r) return; int mid = (l + r) >> 1; Build_Tree(lson, l, mid); Build_Tree(rson, mid + 1, r); } void Push_Up(int x) { int l = tree[x].l,r = tree[x].r; if(tree[x].sum) { tree[x].len = X[r + 1] - X[l]; tree[x].lc = tree[x].rc = true; tree[x].c = 1; } else { tree[x].len = tree[lson].len + tree[rson].len; tree[x].lc = tree[lson].lc, tree[x].rc = tree[rson].rc; tree[x].c = tree[lson].c + tree[rson].c; if(tree[lson].rc && tree[rson].lc) tree[x].c -= 1; } } void Update_Tree(int x, int L, int R, int c) { int l = tree[x].l, r = tree[x].r; if(X[l] >= R || X[r + 1] <= L) return; if(L <= X[l] && X[r + 1] <= R) { tree[x].sum += c; Push_Up(x); return; } Update_Tree(lson, L, R, c); Update_Tree(rson, L, R, c); Push_Up(x); } ScanLine make_line(int l, int r, int h, int mark) { ScanLine res; res.l = l, res.r = r, res.h = h, res.mark = mark; return res; } int main() { read(n); for(int i = 1; i <= n; i++) { read(x1),read(yy),read(x2),read(y2); line[i * 2 - 1] = make_line(x1, x2, yy, 1); line[i * 2] = make_line(x1, x2, y2, -1); X[i * 2 - 1] = x1, X[i * 2] = x2; } n = n << 1; sort(line + 1, line + n + 1); sort(X + 1, X + n + 1); int tot = unique(X + 1, X + n + 1) - X - 1; Build_Tree(1, 1, tot - 1); int res = 0; for(int i = 1; i < n; i++) { Update_Tree(1, line[i].l, line[i].r, line[i].mark); res = res + abs(pre - tree[1].len); pre = tree[1].len; res = res + 2 * tree[1].c * (line[i + 1].h - line[i].h); } res = res + line[n].r - line[n].l; cout<<res<<endl; return 0; }

     👆：只能说上面求周长并的代码跟面积并的代码非常相似，可以相互比较理解一下

         那么简单的扫描线算法到这里为止，我们已经会使用扫描线算法来求得多个矩形的面积并和周长并 —— 非常基础的一些计算几何问题，当然这里还涉及到了有关线段树加速的知识，如果不了解线段树的同学通过这个算法也可以加深对线段树的理解[doge]。

         OK，上点强度~✌，接下来让我们想想除了矩形，三角形，梯形，圆这些几何图形的面积并是否也可以用扫描线求出？🤔

         👆：等腰三角形的扫描线模型

          👆：圆离散后的扫描线模型

     我们仔细看一看上面的两种模型，都是离散+扫描，经典的扫描线算法啊，可行性映入眼帘了， 所以~

         ————  答案是可以的，因为之前提过扫描线算法的思想其实接近于微积分，即求出每个单位内的微分作积，而求任何几何图形的面积，在数学里我们使用的最多的就是微积分的方法，计算几何中，我们也常常使用自适应辛普森积分公式，格林公式等微积分的方法来求相关值，所以扫描线算法这种模拟微分也是可以的，并且精度方面占点优势：

             👇：辛普森积分公式

                👇：格林公式

        这边以圆的面积并为例，我们采用三种方法，一种是采用扫描线+辛普森公式，一种是圆的离散化，还有一种是基本的几何方法~，👇选择性理解一下叭，这个还是有点难度~：

           1.圆的离散化👇：

#include<bits/stdc++.h> using namespace std; const int maxn = 55; const int MAXN = maxn * maxn + 3 * maxn; const double Eps = 1e-10; #define Temp template<typename T> Temp inline void read(T &x) { x = 0; T w = 1,ch = getchar(); while(!isdigit(ch) && ch != '-') ch = getchar(); if(ch == '-') w = -1,ch = getchar(); while(isdigit(ch)) x = (x << 3) + (x << 1) + (ch ^ '0'),ch = getchar(); x = x * w; } int Sign(double x) { if(fabs(x) < Eps) return 0;if(x > 0) return 1;if(x < 0) return -1;} int Dcmp(double x,double y) { if(fabs(x - y) < Eps) return 0;if(x > y) return 1;if(x < y) return -1;} double sqr(double x) { return x * x;} struct Point{ double x,y; Point(double x = 0,double y = 0) : x(x) , y(y) {}; }point[MAXN]; typedef Point Vector; Vector operator + (Vector Alpha,Vector Beta) { return Vector(Alpha.x + Beta.x,Alpha.y + Beta.y);} Vector operator - (Vector Alpha,Vector Beta) { return Vector(Alpha.x - Beta.x,Alpha.y - Beta.y);} Vector operator * (Vector Alpha,double x) { return Vector(Alpha.x * x,Alpha.y * x);} Vector operator / (Vector Alpha,double x) { return Vector(Alpha.x / x,Alpha.y / x);} double Dis(Point Alpha,Point Beta) { return sqrt(sqr(Alpha.x - Beta.x) + sqr(Alpha.y - Beta.y)); } struct Tcir{ double r; Point o; }; struct Tinterval { double x,y,Area,mid; int type; Tcir owner; void area(double l,double r) { double len = sqrt(sqr(l - r) + sqr(x - y)); double d = sqrt(sqr(owner.r) - sqr(len) / 4.0); double angle = atan(len / 2.0 / d); Area = fabs(angle * sqr(owner.r) - d * len / 2.0); } }inter[maxn]; double x[MAXN],l,r; int n,N,Nn; bool cmpR(Tcir Alpha,Tcir Beta) { return Alpha.r > Beta.r; } void Get(Tcir owner,double x,double &l,double &r) { double y = fabs(owner.o.x - x); double d = sqrt(fabs(sqr(owner.r) - sqr(y))); l = owner.o.y + d; r = owner.o.y - d; } void Get_Interval(Tcir owner,double l,double r) { Get(owner,l,inter[Nn].x,inter[Nn + 1].x); Get(owner,r,inter[Nn].y,inter[Nn + 1].y); Get(owner,(l + r) / 2.0,inter[Nn].mid,inter[Nn + 1].mid); inter[Nn].owner = inter[Nn + 1].owner = owner; inter[Nn].area(l,r),inter[Nn + 1].area(l,r); inter[Nn].type = 1,inter[Nn + 1].type = -1; Nn = Nn + 2; } bool cmp(Tinterval Alpha,Tinterval Beta) { return Alpha.mid > Beta.mid + Eps; } void Add(double xx) { x[N] = xx; N = N + 1; } double dist2(Point Alpha,Point Beta) { return sqr(Dis(Alpha,Beta)); } Point Rotate(Point p,double cost,double sint) { double x = p.x,y = p.y; return Point(x * cost - y * sint,x * sint + y * cost); } pair<Point,Point>crosspoint(Point ap,double ar,Point bp,double br) { double d = (ap - bp).norm(); double cost = (ar * ar + d * d - br * br) / (2 * ar * d); double sint = sqrt(1. - cost * cost); Point v = (bp - ap) / (bp - ap).norm() * ar; return make_pair(ap + Rotate(v,cost,-sint),ap + Rotate(v,cost,sint)); } double getUnion(int n,Tcir a[]) { int p = 0; sort(a,a + n,cmpR); for(int i = 0;i < n;i++) { bool f = true; for(int j = 0;j < i;j++) if(dist2(a[i].o,a[j].o) <= sqr(a[i].r - a[j].r) + Eps) { f = false; break; } if(f == true) { a[p] = a[i]; p = p + 1; } } n = p; N = 0; for(int i = 0;i < n;i++) { Add(a[i].o.x - a[i].r); Add(a[i].o.x + a[i].r); Add(a[i].o.x); for(int j = i + 1;j < n;j++) if(dist2(a[i].o,a[j].o) <= sqr(a[i].r + a[j].r) + Eps) { pair<Point,Point> cross=crosspoint(a[i].o,a[i].r,a[j].o,a[j].r); Add(cross.first.x); Add(cross.second.x); } } sort(x,x + N); p = 0; for(int i = 0;i < N;i++) if(! i || fabs(x[i] - x[i - 1]) > Eps) { x[p] = x[i]; p = p + 1; } N = p; double ans = 0; for(int i = 0;i + 1 < N;i++) { l = x[i],r = x[i + 1]; Nn = 0; for(int j = 0;j < n;j++) if(fabs(a[j].o.x - l) < a[j].r + Eps && fabs(a[j].o.x - r) < a[j].r + Eps) Get_Interval(a[j],l,r); if(Nn) { sort(inter,inter + Nn,cmp); int cnt = 0; for(int i = 0;i < Nn;i++) { if(cnt > 0) { ans = ans + (fabs(inter[i - 1].x - inter[i].x) + fabs(inter[i - 1].y - inter[i].y)) * (r - l) / 2.0; ans = ans + inter[i - 1].type * inter[i - 1].Area; ans = ans - inter[i].type * inter[i].Area; } cnt = cnt + inter[i].type; } } } return ans; } int main() { read(n); Tcir Circle[MAXN]; for(int i = 0;i < n;i++) { double xxx,yyy,rrr; scanf("%lf%lf%lf",&xxx,&yyy,&rrr); Circle[i].o.x = xxx,Circle[i].o.y = yyy; Circle[i].r = rrr; } printf("%.10f\n",getUnion(n,Circle)); return 0; }

         2：👇基本的几何方法：

    1.如图所示，将圆的面积剖分成若干个多边形的面积与若干个弓形的面积。多边形的边就是圆的交点构成的不被其他圆覆盖的弦。计算有向面积的话，可以看到中间“洞”的面积恰好被顺时针的多边形包围，因此会被减去。请注意，必须去重复的圆，否则答案会有重复计算的面积。

        代码：

double Cross(Point Alpha,Point Beta) { return Alpha.x * Beta.y - Alpha.y * Beta.x; } struct Circle{ Point p; double r; bool operator < (Circle o) { if(Sign(r - o.r) != 0) return Sign(r - o.r) == -1; if(Sign(p.x - o.p.x) != 0) { return Sign(p.x - o.p.x) == -1; } return Sign(p.y - o.p.y) == -1; } bool operator == (Circle o) { return Sign(r - o.r) == 0 && Sign(p.x - o.p.x) == 0 && Sign(p.y - o.p.y) == 0; } }; inline pair<Point,Point>crosspoint(Circle Alpha,Circle Beta) { return crosspoint(Alpha.p,Alpha.r,Beta.p,Beta.r); } Circle c[1000],tc[1000]; int n,m; struct Node{ Point p; double a; int d; Node(Point p,double a,int d) : p(p) , a(a) , d(d) {} bool operator < (Node o) { return a < o.a; } }; double arg(Point p) { return arg(complex<double>(p.x,p.y)); } double solve() { sort(tc,tc + m); m = unique(tc,tc + m) - tc; for(int i = m - 1;i >= 0;i--) { bool ok = true; for(int j = i + 1;j < m;j++) { double d = (tc[i].p - tc[j].p).norm(); if(Sign(d - abs(tc[i].r - tc[j].r)) <= 0) { ok = false; break; } } if(ok == true) { c[n] = tc[i]; n = n + 1; } } double ans = 0; for(int i = 0;i < n;i++) { vector<Node> event; Point boundary = c[i].p + Point(-c[i].r,0); event.push_back(Node(boundary,-pi,0)); event.push_back(Node(boundary,pi,0)); for(int j = 0;j < n;j++) { if(i == j) continue; double d = (c[i].p - c[j].p).norm(); if(Sign(d - (c[i].r + c[j].r)) < 0) { pair<Point,Point> ret = crosspoint(c[i],c[j]); double x = arg(ret.first - c[i].p); double y = arg(ret.second - c[i].p); if(Sign(x - y) > 0) { event.push_back(Node(ret.first,x,1)); event.push_back(Node(boundary,pi,-1)); event.push_back(Node(boundary,-pi,1)); event.push_back(Node(ret.second,y,-1)); }else { event.push_back(Node(ret.first,x,1)); event.push_back(Node(ret.second,y,-1)); } } } sort(event.begin(),event.end()); int sum = event[0].d; for(int j = 1;j <(int)event.size();j++) { if(sum == 0) { ans = ans + cross(event[j - 1].p,event[j].p) / 2; double x = event[j - 1].a; double y = event[j].a; double area = c[i].r * c[i].r * (y - x) / 2; Point v1 = event[j - 1].p - c[i].p; Point v2 = event[j].p - c[i].p; area = area - Cross(v1,v2) / 2; ans = ans + area; } } } return ans; }

          3：辛普森积分👇：

      那么小有难度的圆的面积并也被我们轻松解决了，ok，这篇博文实在是有点长了，再讲三角形和梯形就没必要了毕竟大家思路都懂了对叭~[doge] ，所以给几道可以去试试的例题：

           三角形 – 洛谷

           [HNOI2012] 三角形覆盖问题 – 洛谷 

           POJ 1177

           POJ 1151 

    这个扫描线算法小有难度但是真的很有趣啊，跟旋转卡壳一样有趣不是吗[doge]，我要好好学旋转卡壳辣，今天就先到这里~