本系列文章力求以简洁易懂的文字介绍计算几何中的基本概念,使读者快速入门,故不追求难度和深度,仅起到抛砖引玉的作用。
那么书接上回(计算几何系列 ——— 闭合曲线之美(圆类-上)-CSDN博客),我们继续来讲解有关二维计算几何圆类的一些基本问题。
一.最小圆覆盖(可参考最小圆覆盖 - 洛谷)
给出 N 个点,让你画一个最小的包含所有点的圆。(1 <= N <= 1e5, xi , yi < 1e4,误差不超过10的负9次),那么在思考这个问题的时候,有几个非常常见的不一定正确的思路(本蒟蒻想的,很乐色的思路)。
1.三重循环,枚举最小圆覆盖的内接三角形。但是很明显,这种思路包包TLE(N^3),并且假如所有点都共线的情况呢?是不是不存在最大内接三角形的说法了🤔,所以是不可行的;
2.第一个思路错误了,所以更加暴力些,双重循环直接求最小圆覆盖的直径,然后判断剩下的点是否在圆内,这样的思路(⊙﹏⊙),在某些情况下又是不成立的,而且其实时间复杂度也是O(N^3),妥妥TLE。
于是乎,我们推出一种全新的玄幻的算法——随机增量法 👇
但是我们还是采用三点定圆,当然要三层循环,枚举三个点。
1.对于第一个点,我们让其为圆心。
2.对于第二点我们和第一个点连线,取其中点为圆心.连线为直径。
3.对于第三个点,我们进行三点定圆.
来了第四个点,我们判断其是否在圆中,如果在圆中,我们继续判断,如果不在圆中,我们取一二个点,和第四个点进行三点定圆. 因为第四个点不在前三个点的圆中.所以第四个点定的新圆一定比前三个定的圆大.
有关时间复杂度的问题:
显然,第三个点的时间复杂度是o(n)的,表面上时间复杂的是O(n^3) 的,实际上在第一次执行完最内层循环.已经确定了一个圆。再执行外边的两层循环,就会涉及到概率,大概有如下的关系:
所以说我们再加一个玄学优化:random_shuffle(a+1,a+n+1);把出题人善意给你的数据顺序打乱。(之所以说它是玄学,是因为我不加也过了😄)
代码是非常的简单:
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 1e5 + 10;
const double Eps = 1e-9;
int N,stk1[MAXN],stk2[MAXN],top1,top2,top,stk[MAXN];
double x,y,ans,ansi,ansj,gh;
#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));
}
Point Circum(Point Alpha,Point Beta,Point Gama)
{
double x1 = Alpha.x,y1 = Alpha.y;
double x2 = Beta.x,y2 = Beta.y;
double x3 = Gama.x,y3 = Gama.y;
double a1 = 2 * (x2 - x1);
double b1 = 2 * (y2 - y1);
double c1 = x2 * x2 + y2 * y2 - x1 * x1 - y1 * y1;
double a2 = 2 * (x3 - x2);
double b2 = 2 * (y3 - y2);
double c2 = x3 * x3 + y3 * y3 - x2 * x2 - y2 * y2;
double x = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
double y = (a1 * c2 - a2 * c1) / (a1 * b2 - a2 * b1);
return Point(x,y);
}
struct Circle{ Point Center; double R;}circles[MAXN];
int main()
{
read(N);
for(int i = 1;i <= N;i++)
{
scanf("%lf%lf",&x,&y);
point[i].x = x;
point[i].y = y;
}
circles[1].Center = point[1];
circles[1].R = 0;
for(int i = 2;i <= N;i++)
{
gh = Dis(point[i],circles[1].Center);
if(Dcmp(gh,circles[1].R) == 1)
{
circles[1].Center = point[i];
circles[1].R = 0;
for(int j = 1;j <= i - 1;j++)
{
gh = Dis(point[j],circles[1].Center);
if(Dcmp(gh,circles[1].R) == 1)
{
circles[1].Center = (point[i] + point[j]) / 2;
circles[1].R = Dis(point[i],point[j]) / 2;
for(int k = 1;k <= j - 1;k++)
{
gh = Dis(point[k],circles[1].Center);
if(Dcmp(gh,circles[1].R) == 1)
{
circles[1].Center = Circum(point[i],point[j],point[k]);
circles[1].R = Dis(point[i],circles[1].Center);
}
}
}
}
}
}
printf("%.10f\n",circles[1].R);
printf("%.10f",circles[1].Center.x);
cout<<' ';
printf("%.10f\n",circles[1].Center.y);
return 0;
}
这边给大家一个值得思考的问题:有没有一种优化的可能性,就是说1e5个点还是太多了,可不可以消去一些不必要出现的点,那么首先我们可不可以求出这个二维凸包,相信最小覆盖圆一定在这个最大二维凸包上呢?然后再次采用随机增量法来求最小覆盖圆,大家想想这个到底对不对?[doge]👆
练习:[AHOI2012] 信号塔 - 洛谷
二.求圆的面积并的两种常见思路
我们先来看一个常见的问题,就是求圆之间的面积并,假如是求两圆之间的面积并,就是非常简单的容斥原理问题:
一般思路就是我们采用余弦定理求出两个夹角,分别求出两倍角,然后把两边部分的扇形面积求出来,相加减去四边形ABCD的面积,设这部分面积为S1,两圆的面积分别为S2和S3,所以有S = S2 + S3 - S1。但是假如我们有N个圆呢?
求面积并?我们可以想到一个耳熟能详的东西,那就是扫描线算法🤔(类似于微积分):
扫描线:假设有一条扫描线从一个图形的下方扫向上方(或者左方扫到右方),那么通过分析扫描线被图形截得的线段就能获得所要的结果。该过程可以用线段树进行加速。
对于多个矩形求面积并的时候,我们一般将上下的两条边作为扫描对象:
扫描线的例题可以参考:三角形 - 洛谷 【模板】扫描线 - 洛谷
那么如何将多个圆做类似扫描线的处理呢? 一般的方法如下:
将圆与圆之间的交点,圆的左右端点,以及圆心的x值离散出来作为事件点,并按照这些值,画竖直的线切割圆。这样相邻两条竖线之间,每个圆要么不在这里,要么有一上一下两个圆弧。这些弧之间除了在竖线上就再无交点。取每条弧的中点,排序后进行扫描即可。
我们使用getUnion中的cnt记录了当前区域被多少个圆所覆盖。因此可以通过适当的修改,求出恰好被k个圆覆盖的区域的面积。
上述的这个类似于扫描线的算法,又称圆的离散化~:
👆:是对圆的各个需要离散化的坐标值进行竖线切割(貌似漏了一条线,不要在意这些细节[doge])~
👆:对于这个阴影面积的求法,就是最基本的几何问题了
深度学习圆的离散化可参考:与圆有关的离散化方法 - 百度文库
在圆的离散化中,我求出了每条弧线的中点之后即可直接计算,无需再用线段树进行维护(有些计算几何预处理的基本函数没有给出,希望能自己练一下~[doge]):
#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;
}
注释:常数中,maxn为圆的个数,MAXN为所有事件点的个数~
另外,新增了圆类Tcir和中间过程所需的区间类Tinterval ~
.norm()表示向量的模长~
时间复杂度:O(n^3 * logn)
那么对于圆的面积并,我们还有一个非常基本的思路:
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;
}
注释:.norm()表示向量的模长
OK,那么关于圆类的一些问题我们差不多都讲完了,意味着我们二维计算几何差不多将告一段落,之前说这些不过仅仅是入门的知识点,ACM-ICPC竞赛考到的会比这个难很多,深基只是为了了解基础概念,掌握基础思路,但是现场赛中的思路与算法光靠这些是远远不够的。
那么我们接下来就开始讲三维几何了吗?o.O?思考?(⊙﹏⊙)暂等蒟蒻本人好好学习一番,那么计算几何系列的文章从这里开始接轨赛场原题,并且提供更多灵活有用的思路,我会不定时的推送一些较难的计算几何真题以及题解,如果有什么问题都可以来问我~接下来推几道较为经典的计算几何:
[USACO03FALL] Beauty Contest G /【模板】旋转卡壳 - 洛谷
以上问题将会大上强度,我们先别急着学三维计算几何,我们先好好掌握一下二维的~,而且近几年的ACM里三维计算几何的出现概率也在下降,几何已不是重点,主要考的是思维了,希望大家能在我的文章中学到真功夫,做出这些题~😄