本系列文章力求以简洁易懂的文字介绍计算几何中的基本概念，使读者快速入门，故不追求难度和深度，仅起到抛砖引玉的作用。
那么书接上回（计算几何系列 ——— 闭合曲线之美（圆类-上）-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里三维计算几何的出现概率也在下降，几何已不是重点，主要考的是思维了，希望大家能在我的文章中学到真功夫，做出这些题~😄
6634
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
