望远镜问题
解题思路:
本题是求一个简单多边形和圆的交集的面积,这里可以用三角剖分来求。
三角剖分是一个思想,可以用来求一些计算几何的问题。这里讲解如何用三角剖分的思想来求圆和简单多边形的面积并
设圆的圆心为 $p$,半径为 $r$
我们为多边形的每条边定义一个逆时针的方向,则对于每条边 $\vec{ab}$ 都能和圆心构成一个三角形 $pab$,我们定义这个三角形的有向面积为 $\vec{pa} \times \vec{pb}$,则此时整个多边形的面积就是所有边和圆心构成的三角形的有向面积的和。这是一个求简单多边形面积的常用方法。
这样我们能很容易得到多边形和圆的面积,而他们的面积并就应该是多边形的面积加上圆的面积再减去多边形和圆的交集的面积。
因此我们接下来的问题就是求出多边形和圆的交集的面积。我们还是看多边形每条边和圆心构成的三角形,此时我们不计算每个三角形的面积,我们计算每个三角形和圆的交集的有向面积,如果我们能求出每个三角形和圆的交集的有向面积,则所有三角形和圆的交集的有向面积的和恰好就是多边形和圆的交集的面积。可以自行画图理解,重点在于有向面积。
然后问题就变成了对于每个三角形,求一下它和圆的交集的面积,这样我们将问题从求一个多边形和圆的交集变成了求一个三角形和圆的交集,求多边形和圆的交集情况非常的多,完全考虑不完,但是三角形就比较固定了,能考虑的情况也比较少,这样我们就可以进一步分情况来讨论,这也是三角剖分的核心思想。
然后考虑三角形 $pab$ 和圆的交集的所有情况。
第一种情况,$a,b$ 两点都在圆内,此时三角形和圆的交集就是三角形的面积,直接用叉积计算即可。
第二种情况,$a,b$ 都在圆外且边 $ab$ 不经过圆,设 $pb$ 和 $pa$ 的夹角为 $\theta$,此时的交集就是一个半径为 $r$ 角度为 $\theta$ 的扇形,也可以直接计算。
第三种情况,$a$ 在圆外,$b$ 在圆内,设 $ab$ 和圆的交点是 $q$,则此时我们可以将交集分成一个三角形和一个扇形,三角形用叉积,扇形求一下夹角套扇形公式计算。
第四种情况,$a$ 在圆内,$b$ 在圆外,这一个情况和第三种情况类似,不再多说。
第五种情况,$a,b$ 都在圆外且边 $ab$ 经过圆,此时的交集可以分成两个扇形和一个三角形,三个部分分开计算即可。
到此我们就将所有三角形和圆的交集的情况都考虑清楚了,接下来只需要根据不同的情况来计算交集的面积即可,最终将每个三角形和圆的交集的有向面积加在一起就是多边形和圆的交集的面积了。
在求三角形和圆的交集的时候,我们需要快速的知道边和圆的两个交点 $a,b$ 的位置,可以联立方程然后解方程求出来,这里用向量的方式来求,我们可以先求一下 $p$ 到 $ab$ 的垂足 $e$,$p$ 加上 $ab$ 顺时针旋转 $90^。$ 后的向量就是 $p$ 到 $ab$ 的垂线,求一下垂线和 $ab$ 的交点就是 $e$。然后求一下 $p$ 到 $ab$ 的距离,就可以用勾股定理求出 $ea$ 的距离,我们知道 $ea$ 的距离和方向,就能通过 $e$ 的坐标得到 $a$ 的坐标。$b$ 就是 $e$ 往反方向走同样距离得到的坐标。
另外,区分第二种情况和第五种情况我们可以通过判断 $p$ 在 $ab$ 的垂足 $e$ 是否在 $ab$ 上,如果在 $ab$ 上,说明是第五种情况,否则是第二种情况。
以上就是三角剖分求圆和多边形的面积并的具体思路和细节,其中涉及很多计算几何的常用操作,可以参考 计算几何的基础知识
注意本题给出的多边形的点可能是顺时针方向的也可能是逆时针方向的,如果是逆时针方向的那么我们最终求出的有向面积应该是正的,如果是顺时针方向的那么我们最终求出的有向面积应该是负的,因此最终答案需要取一个绝对值。
C++ 代码
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
const int N = 55;
const double eps = 1e-8;
const double PI = acos(-1);
double R;
int n;
PDD q[N]; //存储多边形所有顶点
PDD p; //圆心
int sign(double x) //符号函数
{
if(fabs(x) < eps) return 0;
if(x < 0) return -1;
return 1;
}
int dcmp(double x, double y) //比较函数
{
if(fabs(x - y) < eps) return 0;
if(x < y) return -1;
return 1;
}
PDD operator+(PDD a, PDD b) //重载向量加法
{
return {a.x + b.x, a.y + b.y};
}
PDD operator-(PDD a, PDD b) //重载向量减法
{
return {a.x - b.x, a.y - b.y};
}
PDD operator*(PDD a, double t) //重载数乘
{
return {a.x * t, a.y * t};
}
PDD operator/(PDD a, double t) //重载数除
{
return {a.x / t, a.y / t};
}
double operator*(PDD a, PDD b) //重载叉乘
{
return a.x * b.y - a.y * b.x;
}
double operator&(PDD a, PDD b) //重载点乘
{
return a.x * b.x + a.y * b.y;
}
double area(PDD a, PDD b, PDD c) //求 ab 和 ac 构成的平行四边形的有向面积
{
return (b - a) * (c - a);
}
double get_len(PDD a) //求向量 a 的长度
{
return sqrt(a & a);
}
double get_dist(PDD a, PDD b) //求 a 和 b 的距离
{
return get_len(b - a);
}
double project(PDD a, PDD b, PDD c) //求 ac 在 ab 上的投影的长度
{
return ((b - a) & (c - a)) / get_len(b - a);
}
PDD rotate(PDD a, double b) //将 a 顺时针旋转 b 弧度
{
return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}
PDD norm(PDD a) //求向量 a 的单位向量
{
return a / get_len(a);
}
bool on_segment(PDD p, PDD a, PDD b) //判断 p 是否在 ab 上
{
return !sign((a - p) * (b - p)) && sign((a - p) & (b - p)) <= 0; //pa, pb 平行并且方向相反,说明 p 在 ab 上
}
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w) //求两个点向式直线的交点
{
PDD u = p - q;
double t = (w * u) / (v * w);
return p + v * t;
}
//求圆和直线 ab 的两个交点 pa, pb,并返回 mind
//如果是情况 5,mind 表示 p 到 ab 的垂线的长度,否则 mind 表示 p 到 a 和 b 的最短距离
double get_circle_line_intersection(PDD a, PDD b, PDD &pa, PDD &pb)
{
PDD e = get_line_intersection(a, b - a, p, rotate(b - a, PI / 2)); //求 p 到 ab 的垂足
double mind = get_dist(p, e); //计算 p 到 e 的距离
if(!on_segment(e, a, b)) mind = min(get_dist(p, a), get_dist(p, b)); //如果不是情况 5,修改 mind 的定义
if(dcmp(R, mind) <= 0) return mind; //如果是情况 2,则不需要用到 pa, pb,直接返回
double len = sqrt(R * R - get_dist(p, e) * get_dist(p, e)); //勾股定理求 e 到 ab 和圆的交点的距离
pa = e + norm(a - b) * len; //求出 pa 的坐标
pb = e + norm(b - a) * len; //求出 pb 的坐标
return mind;
}
double get_sector(PDD a, PDD b) //求 a 和 b 构成的半径为 R 的有向扇形面积
{
double angle = acos((a & b) / get_len(a) / get_len(b)); //计算 a 和 b 的夹角弧度
if(sign(a * b) < 0) angle = -angle; //如果 a 和 b 构成的有向面积是负的,需要将弧度取反
return R * R * angle / 2;
}
double get_circle_trangle_area(PDD a, PDD b) //求圆和三角形 pab 的交集的面积
{
double da = get_dist(p, a), db = get_dist(p, b); //求 a 和 b 到圆心的距离
if(dcmp(R, da) >= 0 && dcmp(R, db) >= 0) return a * b / 2; //情况 1
if(!sign(a * b)) return 0; //如果 p, a, b 三点共线,面积为 0
PDD pa, pb; //记录 ab 和圆的交点,pa 距离 a 更近,pb 距离 b 更近
double mind = get_circle_line_intersection(a, b, pa, pb); //求 pa, pb,并求出 p 到 ab 的最短距离 mind
if(dcmp(R, mind) <= 0) return get_sector(a, b); //情况 2
if(dcmp(R, db) >= 0) return pa * b / 2 + get_sector(a, pa); //情况 3
if(dcmp(R, da) >= 0) return a * pb / 2 + get_sector(pb, b); //情况 4
return get_sector(a, pa) + pa * pb / 2 + get_sector(pb, b); //情况 5
}
double work() //求圆和多边形的交集的面积
{
double res = 0; //记录交集的面积
for(int i = 0; i < n; i++) res += get_circle_trangle_area(q[i], q[(i + 1) % n]);
return fabs(res); //返回面积的绝对值
}
int main()
{
while(scanf("%lf%d", &R, &n) != -1)
{
for(int i = 0; i < n; i++) scanf("%lf%lf", &q[i].x, &q[i].y);
printf("%.2lf\n", work());
}
return 0;
}