扫描线与自适应辛普森积分求多边形面积并

亦凉 2022-11-14 13:29 186阅读 0赞

### 求多边形面积并 ###

就是说在二维平面上给你一堆多边形，然后求所有面积的并集，重叠的地方只算一次。就相当于一张桌子上随意放了一堆多边形，求占桌子的总面积。  
求面积并问题一般用两种方法，扫描线和自适应辛普森积分。

--------------------

### 扫描线 ###

#### 矩形面积并 ####

就是从左到右，用一条条**竖直的线将**所有多边形**按照竖直的边**分成一个个小矩形，一般适用于各边都是直的多边形，例如矩形，三角形。具体一点可以看下图：  
![eg][]  
有四个黑色的矩形随机的放在一起，然后我们用这些红线从所有的矩形的竖直边将整个图分成若干长条，**使得每个长条内没有矩形的竖直边**。这些就是扫描线。这样操作之后，我们的面积并就好求了，只需要计算所有红线之间的小矩形面积即可。  
小矩形的面积就是长乘宽，宽就是两条红线之间的距离。长就是所有在**两条红线内的**矩形的长的**区间合并**的总长。如图：  
![eg][eg 1]  
我们挨着求完所有红线之间的面积之和，就是所有的面积并了。

#### 代码模板 ####

原题链接：[Acwing 3068. 扫描线][Acwing 3068.]

#include<iostream>
    #include<vector>
    #include<algorithm>
    #define x first //方便使用pair
    #define y second
    using namespace std;
    const int N = 1010;
    typedef pair<int,int> PII;	//用pair来存坐标，所有坐标都是整数
    typedef long long LL;		//答案要用long long存
    
    int n;
    LL ans;		//答案
    PII l[N],r[N],q[N];	//l存左下角的坐标，r存右上角的坐标，q记两条竖线之间的区间
    vector<int> xs;		//xs存所有矩形的竖直边
    
    LL range_area(int a,int b)	//计算竖直线x=a和x=b之间的面积
    { 
    	LL cnt = 0,res = 0;	//cnt记两条竖线内的矩形个数，res记区间合并后的长度
    	for(int i = 0;i < n;i++)	//循环所有矩形
    		if(l[i].x <= a && r[i].x >= b)	//判断若矩形在区间内部
    			q[cnt++] = { l[i].y,r[i].y};	//将该矩形竖直方向的区间存进q
    	if(!cnt)	//如果cnt为0，说明该两条竖线之间没有矩形
    		return 0;	//直接返回0
    	sort(q,q+cnt);	//区间合并，先对q进行排序
    	int st = q[0].x, ed = q[0].y;	//st和ed是当前区间的起点和终点
    	for(int i = 1;i < cnt;i++)	//循环所有区间
    	{ 
    		if(q[i].x <= ed)	//若某区间的起点在当前区间终点的左侧
    			ed = max(ed,q[i].y);	//更新当前区间的终点，是这两个区间终点的最大值
    		else	//如果不是，说明区间断开了
    		{ 
    			res += ed - st;	//res就加上当前区间长度
    			st = q[i].x;	//用该区间更新当前区间
    			ed = q[i].y;
    		}
    	}
    	res += ed - st;		//res一定要加上最后一个区间的长度
    	return res * (b-a);	//返回区间合并的长度乘宽度
    }
    
    int main()
    { 
    	cin >> n;
    	for(int i = 0;i < n;i++)
    	{ 
    		cin >> l[i].x >> l[i].y >> r[i].x >> r[i].y;	//输入矩形的左下角和右上角坐标
    		xs.push_back(l[i].x);		//把所有的竖直边存进xs数组
    		xs.push_back(r[i].x);
    	}
    	sort(xs.begin(),xs.end());	//排序
    	for(int i = 0;i < xs.size()-1;i++)	//循环每两条竖直边，计算之间的面积
    		if(xs[i] != xs[i+1])	//去掉重边
    			ans += range_area(xs[i],xs[i+1]);	//答案加上这之间的面积
    	cout << ans << endl;		//输出答案
    	
    	return 0;
    }

--------------------

### 自适应辛普森积分 ###

#### 圆形面积并 ####

当我们处理圆形面积并时，每个扫描线的小长条内部是许多的圆弧，他们的面积很难算，因此当小长条内的面积不好计算时，我们就用自适应辛普森积分来计算。积分的本质也是算面积，我们把他当成一个函数来积面积就行。

若学过微元法就会很熟悉，如果要求一个随意函数的定积分，我们将这个函数分成很多个小矩形，然后我们可以取这个矩形的左上角或右上角的值来当作这个矩形的高，然后乘上矩形的宽就是面积，所有矩形的面积和就是积分。当这个矩形被分的无限多的时候，他们的面积和就越趋近于函数的积分，这就是微分，如下图：  
![eg][eg 2]  
而辛普森积分也是用微分的方法来积分，只是他不是选的是矩形的一个点，而是将每个矩形的顶端看作一个**二次函数**，然后计算这个的面积。当把矩形的顶看作二次函数时更加精确。辛普森积分也是自适应的，他会**自适应你给的精度**，因为他是一个**递归**的过程。当我们要求\[l，r\]的积分时，他先计算出整段的中点mid，然后计算\[l，r\]的积分记作**S**，然后计算\[l，mid\]的积分记作**left**，然后计算\[mid，r\]的积分记作**right**，然后判断若 S-(left+right) 的误差在eps精度内，就返回这个S。他会判断当前求的积分的精度够不够，不够的话他会再递归的去求left和right两个区间，不够的话再往下递归，最终会被分成很多很细的区间，精度就会满足要求了。  
那么有一个问题，然后求顶部为二次函数的矩形的面积呢？这里直接给出结论：若矩形的左右两条边是x=l和x=r，矩形的面积就等于 ( r − l ) × f ( l ) + 4 × f ( l + r 2 ) + f ( r ) 6 (r-l) \\times \\frac \{f(l) +4\\times f(\\frac \{l+r\} 2) +f(r)\} 6 (r−l)×6f(l)\+4×f(2l\+r)\+f(r)。因为三点确定一条二次函数，就用f(l)和f(r )和f((l+r)/2)三点，通过不定积分变换就能得到上面的公式了。一般直接背过公式用就可以了。  
其实辛普森积分并不难，可以看代码理解

#### 代码模板 ####

原题链接：[Acwing 3069. 圆的面积并][Acwing 3069.]

**思路：** 用自适应辛普森积分去积f(x)即可，主要是求f(x)，f(a)的意思就是x=a的这条直线和图中所有圆相交的长度，也就是这条直线在圆内的总长度，然后我们从左右两个边界去积f(x)就好了。这里求f(x)用计算几何的方法去求。用勾股定理求所有圆和直线相交的区间，然后将这些区间合并就是总长度。

#include<iostream>
    #include<cmath>
    #include<algorithm>
    #define x first 
    #define y second
    using namespace std;
    const int N = 1010;
    const double eps = 1e-8;
    typedef pair<double,double> PDD;	//pair存坐标
    
    struct Circle{ 		//存圆的结构体
    	PDD o;		//圆心o
    	double r;	//半径r
    }c[N];
    int n;
    PDD q[N];	//存所有区间
    
    int judge(double a,double b)	//判断函数
    { 
    	if(fabs(a-b) < eps)	//ab相等返回0
    		return 0;
    	if(a < b)		//a小于b返回-1
    		return -1;
    	return 1;		//a大于b返回1
    }
    
    double f(double x)	//f函数
    { 
    	int cnt = 0;	//存所有区间
    	for(int i = 0;i < n;i++)	//循环所有圆
    	{ 
    		double X = fabs(c[i].o.x - x), R = c[i].r;	//X是圆心到直线距离，R是半径
    		if(judge(X,R) < 0)	//如果X小于半径R，说明直线和圆相交
    		{ 
    			double Y = sqrt(R*R - X*X);	//用勾股定理算出Y
    			q[cnt++] = { c[i].o.y - Y, c[i].o.y + Y};	//圆心的y加减Y就是区间
    		}
    	}
    	if(!cnt)	//如果没有相交的区间
    		return 0;	//就返回0
    	sort(q,q+cnt);		//区间排序
    	double st = q[0].x, ed = q[0].y, res = 0;	//区间合并
    	for(int i = 1;i < cnt;i++)
    	{ 
    		if(q[i].x <= ed)
    			ed = max(ed,q[i].y);
    		else
    		{ 
    			res += ed - st;
    			st = q[i].x;
    			ed = q[i].y;
    		}
    	}
    	return res + ed - st;	//不要忘记加最后一个区间
    }
    
    double simpson(double l,double r)	//辛普森积分
    { 	
    	double mid = (l+r) / 2;
    	return (r-l) * (f(l) + 4*f(mid) + f(r)) / 6;	//公式
    }
    
    double asr(double l,double r,double s)	//递归自适应
    { 
    	double mid = (l+r) / 2;
    	double left = simpson(l,mid), right = simpson(mid,r);
    	if(fabs(s - left - right) < eps)
    		return left + right;
    	return asr(l,mid,left) + asr(mid,r,right);
    }
    
    int main()
    { 
    	cin >> n;
    	for(int i = 0;i < n;i++)
    		cin >> c[i].o.x >> c[i].o.y >> c[i].r;
    	double l = -2000, r = 2000;	//题目边界，积分区间[-2000,2000]
    	printf("%.3lf\n",asr(l,r,simpson(l,r)));
    
    	return 0;
    }

[eg]: /images/20221022/29bbaca59c9f42fc85254cb560434450.png
[eg 1]: /images/20221022/2ba4568ed92f4381bd0836a7f8dc5351.png
[Acwing 3068.]: https://www.acwing.com/problem/content/3071/
[eg 2]: /images/20221022/f4dd6ddd572746bfb8888d66a13817c2.png
[Acwing 3069.]: https://www.acwing.com/problem/content/3072/