Skip to main content

二维计算几何基础

参考资料​

例题​

洛谷 P1183 多边形的面积

给出一个没有缺口的简单多边形，它的边是垂直或者水平的，要求计算多边形的面积。

Code (1)
#include <bits/stdc++.h>
using namespace std;

const int N=105;
struct Point
{
	double x,y;
}a[N];
double Cross(Point A,Point B){return A.x*B.y-A.y*B.x;}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n;
	cin>>n;
	for(int i=1;i<=n;i++)cin>>a[i].x>>a[i].y;
	double ans=0;
	for(int i=1;i<=n;i++)ans+=Cross(a[i],a[i%n+1]);
	cout<<fabs(ans)/2<<'\n';
	return 0;
}
洛谷 P2333 [SCOI2006] 一孔之见

从一个圆孔里看一个凸多边形，为了让看到的面积至少为 SS，孔的半径至少需要多大？

假设孔的圆心固定在 (0,0)(0,0)，且 (0,0)(0,0) 在多边形的内部（而不是外部或边界上）。

Solution

参考资料​

解题思路​

随着圆的半径增大，能看到的面积单调不降，所以可以二分寻找满足条件的最小半径。

将多边形的每个顶点与原点连接，构成多个三角形。每个三角形都有一个对应圆心角的扇形。

通过分类讨论每个三角形与其对应扇形的交集面积，可以计算出能看到的面积。

函数定义​

f1(A,B)​

f1(A,B)f_1(A,B)：表示三角形 △OAB\triangle OAB 的有向面积。

通过向量叉积计算 f1f_1：

f1(A,B)=OA→×OB→2f_1(A,B)=\frac{\overrightarrow{OA}\times \overrightarrow{OB}}{2} 
double f1(Point A,Point B)
{
	return Cross(A,B)/2;
}

f2(A,B,r)​

f2(A,B,r)f_2(A,B,r)：表示由 OAOA 和 OBOB 两条边围成且半径为 rr 的扇形的有向面积。

通过向量的叉积和点积计算出扇形的圆心角 θ=∠AOB\theta=\angle AOB：

θ=arctan⁡(sin⁡∠AOBcos⁡∠AOB)=atan2⁡(OA→×OB→,OA→⋅OB→)\theta=\arctan\left(\frac{\sin\angle AOB}{\cos\angle AOB}\right)=\operatorname{atan2}(\overrightarrow{OA}\times \overrightarrow{OB},\overrightarrow{OA}\cdot \overrightarrow{OB})

根据扇形面积公式计算 f2f_2：

f2(A,B,r)=θ⋅r22=atan2⁡(OA→×OB→,OA→⋅OB→)⋅r22f_2(A,B,r)=\frac{\theta\cdot r^2}{2}=\frac{\operatorname{atan2}(\overrightarrow{OA}\times \overrightarrow{OB},\overrightarrow{OA} \cdot\overrightarrow{OB})\cdot r^2}{2} 
double f2(Point A,Point B,double r)
{
	return atan2(Cross(A,B),Dot(A,B))*r*r/2;
}

calc(A,B,r)​

calc(A,B,r)calc(A,B,r)：表示每个三角形与其对应扇形的交集面积。

接下来通过分类讨论计算 calccalc，设：

S=calc(A,B,r)S=calc(A,B,r)

分类讨论​

Case 1​

如果点 A,BA,B 距离圆心的距离都小于 rr：

∣OA→∣≤r∧∣OB→∣≤r|\overrightarrow{OA}|\le r\land|\overrightarrow{OB}|\le r

此时扇形完全包含三角形，面积交集就是三角形的面积：

S=f1(A,B)S=f_1(A,B) 
if(Len(A)<=r&&Len(B)<=r)return f1(A,B);

判断相交​

接下来我们需要判断直线 ABAB 是否与圆有交点 PP。

直线 ABAB 的参数方程：

OP→=OA→+t⋅AB→\overrightarrow{OP}=\overrightarrow{OA}+t\cdot\overrightarrow{AB}

代入圆的方程：

r=∣OP→∣=∣OA→+t⋅AB→∣r=|\overrightarrow{OP}|=|\overrightarrow{OA}+t\cdot\overrightarrow{AB}|

整理为标准二次方程：

at2+bt+c=0at^2+bt+c=0

计算方程的系数：

a=∣AB∣2,b=2⋅OA→⋅AB→,c=∣OA∣2−r2a=|AB|^2,b=2\cdot\overrightarrow{OA}\cdot\overrightarrow{AB},c=|OA|^2-r^2 
double a=Len2(B-A),b=Dot(A,B-A)*2,c=Len2(A)-r*r;

计算判别式 Δ\Delta：

Δ=b2−4ac\Delta=b^2-4ac 
double d=b*b-a*c*4;

Case 2​

如果直线与圆没有交点或相切：

Δ≤0\Delta\le 0

此时三角形完全包含扇形，面积交集就是扇形的面积：

S=f2(A,B,r)S=f_2(A,B,r) 
if(d<=0)return f2(A,B,r);

计算交点​

接下来通过求根公式，计算出 tt：

t1=−b−Δ2a,t2=−b+Δ2at_1=\frac{-b-\sqrt{\Delta}}{2a},t_2=\frac{-b+\sqrt{\Delta}}{2a} 
double t1=(-b-sqrt(d))/(a*2),t2=(-b+sqrt(d))/(a*2);

代入直线方程，计算出直线与圆的交点 PP：

OP→=OA→+t⋅AB→\overrightarrow{OP}=\overrightarrow{OA}+t\cdot\overrightarrow{AB}

我们令距离 AA 较近的交点为 CC，距离 BB 较近的交点为 DD：

OC→=OA→+t1⋅AB→,OD→=OA→+t2⋅AB→\overrightarrow{OC}=\overrightarrow{OA}+t_1\cdot\overrightarrow{AB},\overrightarrow{OD}=\overrightarrow{OA}+t_2\cdot\overrightarrow{AB} 
Point C=A+(B-A)*t1,D=A+(B-A)*t2;

Case 3​

如果点 C,DC,D 都在点 A,BA,B 的同一侧（即线段 ABAB 的延长线上）：

t1≥1∨t2≤0t_1\ge 1\lor t_2\le 0

此时三角形也完全包含扇形，面积交集就是扇形的面积：

S=f2(A,B,r)S=f_2(A,B,r) 
if(t1>=1||t2<=0)return f2(A,B,r);

Case 4​

如果点 DD 在线段 ABAB 上：

t1≤0t_1\le 0

此时 ADAD 是三角形面积，DBDB 是扇形面积：

S=f1(A,D)+f2(D,B,r)S=f_1(A,D)+f_2(D,B,r) 
if(t1<=0)return f1(A,D)+f2(D,B,r);

Case 5​

如果点 CC 在线段 ABAB 上：

t2≥1t_2\ge 1

此时 ACAC 是扇形面积，CBCB 是三角形面积：

S=f2(A,C,r)+f1(D,B)S=f_2(A,C,r)+f_1(D,B) 
if(t2>=1)return f2(A,C,r)+f1(C,B);

Case 6​

否则点 C,DC,D 都在线段 ABAB 上。

此时 ACAC 是扇形面积，CDCD 是三角形面积，DBDB 是扇形面积：

S=f2(A,C,r)+f1(C,D)+f2(D,B,r)S=f_2(A,C,r)+f_1(C,D)+f_2(D,B,r) 
return f2(A,C,r)+f1(C,D)+f2(D,B,r);

参考代码​

#include <bits/stdc++.h>
using namespace std;

const double eps=1e-8;
const int N=55;
struct Point
{
	double x,y;
	Point(){}
	Point(double _x,double _y):x(_x),y(_y){}
	Point operator+(Point B){return Point(x+B.x,y+B.y);}
	Point operator-(Point B){return Point(x-B.x,y-B.y);}
	Point operator*(double k){return Point(x*k,y*k);}
}a[N];
double Dot(Point A,Point B){return A.x*B.x+A.y*B.y;}
double Cross(Point A,Point B){return A.x*B.y-A.y*B.x;}
double Len(Point A){return sqrt(Dot(A,A));}
double Len2(Point A){return Dot(A,A);}
double f1(Point A,Point B){return Cross(A,B)/2;}
double f2(Point A,Point B,double r){return atan2(Cross(A,B),Dot(A,B))*r*r/2;}
double calc(Point A,Point B,double r)
{
	if(Len(A)<=r&&Len(B)<=r)return f1(A,B);
	double a=Len2(B-A),b=Dot(A,B-A)*2,c=Len2(A)-r*r;
	double d=b*b-a*c*4;
	if(d<=0)return f2(A,B,r);
	double t1=(-b-sqrt(d))/(a*2),t2=(-b+sqrt(d))/(a*2);
	Point C=A+(B-A)*t1,D=A+(B-A)*t2;
	if(t1>=1||t2<=0)return f2(A,B,r);
	if(t1<=0)return f1(A,D)+f2(D,B,r);
	if(t2>=1)return f2(A,C,r)+f1(C,B);
	return f2(A,C,r)+f1(C,D)+f2(D,B,r);
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n;
	double s;
	cin>>n>>s;
	for(int i=1;i<=n;i++)cin>>a[i].x>>a[i].y;
	double l=0,r=10000;
	while(r-l>eps)
	{
		double mid=(l+r)/2,sum=0;
		for(int i=1;i<=n;i++)sum+=calc(a[i],a[i%n+1],mid);
		if(fabs(sum)>=s)r=mid;
		else l=mid;
	}
	cout<<fixed<<setprecision(2)<<l<<'\n';
	return 0;
}