本系列文章力求以简洁易懂的文字介绍计算几何中的基本概念，使读者快速入门，故不追求难度和深度，仅起到抛砖引玉的作用。 
在我们的计算几何中，乃至基本的几何学中我们都需要解决一类问题，那就是几何图形的面积交并，周长交并问题，在计算几何中我们也常常遇到这样的问题，在这里我们介绍一种神奇的用来解决这类问题的算法，叫做扫描线（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：扫描第一条边：
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] ，所以给几道可以去试试的例题：
这个扫描线算法小有难度但是真的很有趣啊，跟旋转卡壳一样有趣不是吗[doge]，我要好好学旋转卡壳辣，今天就先到这里~
320
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
