正如不知何方大佬所言,计算几何精妙之处,就是不用解析几何的方法去做
为了方便查找,防止自己迷路,我把函数名都写成了拼音
绝对不是因为我英语不好!!!
点和向量都可以用一个坐标\((x,y)\)来表示.
故向量\(Vector\)可以写为
typedef struct point vec;
完整定义如下
{% fold 代码内容 %}
typedef struct point vec; //向量vec
struct point { //点的基本数据结构
double x, y;
point(double _x=0, double _y=0):
x(_x),y(_y)
{
}
point operator*(const point& i_T) const
{
return point(x * i_T.x, y * i_T.y);
}
point operator*(double u) const
{
return point(x * u, y * u);
}
bool operator==(const point& i_T) const
{
return x == i_T.x && y == i_T.y;
}
point operator/(double u) const
{
return point(x / u, y / u);
}
point operator+(const point& i_T)
{
return point(x + i_T.x, y + i_T.y);
}
point operator-(const point& i_T)
{
return point(x - i_T.x, y - i_T.y);
}
friend bool operator<(point a, point b)
{
return a.y == b.y ? a.x < b.x : a.y < b.y;
}
friend ostream& operator<<(ostream& out, point& a)
{
//cout << a.x << ‘ ‘ << a.y;
printf("%.8f %.8f", a.x, a.y);
return out;
}
friend istream& operator>>(istream& in, point& a)
{
in >> a.x >> a.y;
return in;
}
};
{%endfold%}
直线和线段都可以用两个点的坐标来表示
故线段\(Segment\)可以写为
typedef struct Line Segment;
完整定义如下
{% fold 代码内容 %}
typedef struct Line Segment; //线段Segment
struct Line { //直线
vec a, b;
Line(point _a = point(), point _b = point())
: a(_a)
, b(_b)
{
}
friend istream& operator>>(istream& in, Line& a)
{
cin >> a.a >> a.b;
return in;
}
friend ostream& operator<<(ostream& out, Line& a)
{
out << a.a << ‘ ‘ << a.b;
return out;
}
};
{%endfold%}
完整定义如下
{%fold 代码内容 %}
struct cirles {
point o;
double r;
cirles(point _o = point(), double _r = 0.0)
: r(_r)
, o(_o)
{
}
point Point(double t) //圆上任意一点
{
return point(o.x + r * cos(t), o.y + r * sin(t));
}
friend istream& operator>>(istream& in, cirles& a)
{
in >> a.o >> a.r;
return in;
}
friend ostream& operator<<(ostream& out, cirles& a)
{
out << a.o << ‘ ‘ << a.r;
return out;
}
};
{%endfold%}
{%fold 代码内容 %}
const double pi = acos(-1.0);
const double eps = 1e-12;
int zhengfu(double d) //判断正负,即sign/dcmp
{
if (fabs(d) < eps)
return 0;
if (d > 0)
return 1;
return -1;
}
int bijiao(double x, double y) //判断x和y的大小关系
{
if (fabs(x - y) < eps)
return 0;
if (x > y)
return 1;
return -1;
}
int dayu_dengyu(double x, double y) //x>=y否?
{
if (fabs(x - y) < eps || x > y)
return 1;
return 0;
}
double chaji(vec A, vec B) //求向量叉积,即cross
{
return A.x * B.y - A.y * B.x; // 正为A->B左旋
}
double xuanzhuan(point a, point b, point c) //求三点叉积,即side
{
return chaji(b - a, c - a);
}
bool cmp(vec a, vec b) //极角排序,运用叉积的极角排序,相比于atan2慢,但是精度高
{
vec c(0, 0);
if (chaji(a - c, b - c) == 0)
return a.x < b.x;
return chaji(a - c, b - c) > 0;
}
double dianji(vec A, vec B) //点积
{
return A.x * B.x + A.y * B.y;
}
double changdu(vec a) //向量长度
{
return sqrt(a.x * a.x + a.y * a.y);
}
{%endfold%}
求一个点在向量\(\overrightarrow{ab}\)上的投影坐标
设点\(c\),投影在\(\overrightarrow{ab}\)上为\(c‘\),则\(c‘\)的坐标就是:\(cos<\overrightarrow{ac},\overrightarrow{ab}>\times |\overrightarrow{ac}|+a\)
{%fold 代码内容 %}
point touying(Line l, point c) //c投影在直线ab上的位置
{
vec A = l.b - l.a;
vec B = c - l.a;
double La = changdu(A);
double Lc = dianji(A, B) / (La * La);
return A * Lc + l.a;
}
{%endfold%}
求一个点\(c\)关于向量\(\overrightarrow{ab}\)的对称点\(c‘‘\)
先求出\(c\)在\(ab\)上的投影,那么\(c‘‘=2\times \overrightarrow{cc‘}+c\)
{%fold 代码内容 %}
point fanshe(Line l, point c) //求c关于直线ab的对称点c‘
{
point A = touying(l, c);
return (A - c) * 2.0 + c;
}
{%endfold%}
就是..根据图中的判断就是了
{%fold 代码内容 %}
void Counter_Clockwise(point p,point p1,point p2)
{
if (zhengfu(chaji(p2 - p1, p - p1)) == 1)
cout << "COUNTER_CLOCKWISE" << endl;
else if (zhengfu(chaji(p2 - p1, p - p1)) == -1)
cout << "CLOCKWISE" << endl;
else {
if (zhengfu(dianji(p2 - p1, p - p1)) == -1)
cout << "ONLINE_BACK" << endl;
else {
double j = changdu(p2 - p1);
double k = changdu(p - p1);
if (bijiao(j, k) >= 0)
cout << "ON_SEGMENT" << endl;
else
cout << "ONLINE_FRONT" << endl;
}
}
}
{%endfold%}
先判平行,再用点积判垂直
{%fold 代码内容 %}
bool pingxing(vec a,vec b)
{
return bijiao(a.x*b.y,a.y*b.x)==0;
}
void Parallel/Orthogonal(Line l1,Line l2)
{
if (pingxing(l1.b-l1.a,l2.b-l2.a))
cout << "2" << endl;
else {
if (zhengfu(dianji(l1.b-l1.a,l2.b-l2.a))==0)
cout << "1" << endl;
else {
cout << "0" << endl;
}
}
}
{%endfold%}
线段相交要考虑蛮多的,首先,先按照x后y从小到大排一下.
最简单的情况,\(ab\)穿过\(cd\),那么必定有交点.
第二种,\(a\)在\(cd\)上或者\(b\)在\(cd\)上
第三种,共线时,\(a\)在\(cd\)之间或\(b\)在\(cd\)之间.
处理好以上问题,就解决了
{%fold 代码内容 %}
bool bijiao3(vec a, vec b, vec c)
{
if (a.x <= c.x && c.x <= b.x && ((a.y <= c.y && c.y <= b.y || b.y <= c.y && c.y <= a.y)))
return 1;
return 0;
}
bool xianduan_xiangjiao(Segment l1,Segment l2) //两线段是否有交点
{
if (l1.a.x == l1.b.x) {
if (l1.a.y > l1.b.y)
swap(l1.a, l1.b);
} else if (l1.a.x > l1.b.x)
swap(l1.a, l1.b);
if (l2.a.x == l2.b.x) {
if (l2.a.y > l2.b.y)
swap(l2.a, l2.b);
} else if (l2.a.x > l2.b.x)
swap(l2.a, l2.b);
double c1 = chaji(l1.b - l1.a, l2.a - l1.a), d1 = chaji(l1.b - l1.a, l2.b - l1.a);
double c2 = chaji(l2.b - l2.a, l1.a - l2.a), d2 = chaji(l2.b - l2.a, l1.b - l2.a);
if (zhengfu(c1 * d1) < 0 && zhengfu(c2 * d2) < 0) //ab横穿cd
return 1;
else if (zhengfu(c1 * d1) != 0 && zhengfu(c2 * d2) == 0) { //ab不穿过cd
if (zhengfu(c2) == 0) {
if (bijiao3(l2.a, l2.b, l1.a))
return 1;
}
if (zhengfu(d2) == 0)
if (bijiao3(l2.a, l2.b, l1.b))
return 1;
} else if (zhengfu(c1 * d1) == 0 && zhengfu(c2 * d2) != 0) { //cd不穿过ab
if (c1 == 0)
if (bijiao3(l1.a, l1.b, l2.a))
return 1;
if (d1 == 0)
if (bijiao3(l1.a, l1.b, l2.b))
return 1;
} else if (zhengfu(c1 * d1) == 0 && zhengfu(c2 * d2) == 0) { //平行
if (l1.a == l2.a || l1.a == l2.b || l1.b == l2.a || l1.b == l2.b)
return 1;
if (bijiao3(l1.a, l1.b, l2.a) == 1)
return 1;
if (bijiao3(l2.a, l2.b, l1.a) == 1)
return 1;
}
return 0;
}
{%endfold%}
给你两个必定相交的线段,求交点
{%fold 代码内容 %}
point xianduan_jiaodian(Segment l1,Segment l2)//两线段交点
{
double tmpLeft, tmpRight, x = inf, y = inf;
if (xianduan_xiangjiao(l1,l2)) {
tmpLeft = (l2.b.x - l2.a.x) * (l1.a.y - l1.b.y) - (l1.b.x - l1.a.x) * (l2.a.y - l2.b.y);
tmpRight = (l1.a.y - l2.a.y) * (l1.b.x - l1.a.x) * (l2.b.x - l2.a.x) + l2.a.x * (l2.b.y - l2.a.y) * (l1.b.x - l1.a.x) - l1.a.x * (l1.b.y - l1.a.y) * (l2.b.x - l2.a.x);
x = tmpRight / tmpLeft;
tmpLeft = (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) - (l1.b.y - l1.a.y) * (l2.a.x - l2.b.x);
tmpRight = l1.b.y * (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) + (l2.b.x - l1.b.x) * (l2.b.y - l2.a.y) * (l1.a.y - l1.b.y) - l2.b.y * (l2.a.x - l2.b.x) * (l1.b.y - l1.a.y);
y = tmpRight / tmpLeft;
}
return point(x, y);
}
{%endfold%}
{%fold 代码内容 %}
double dian_dao_xianduan(Segment l, point c) //点到线段的距离
{
double L = changdu(l.b - l.a);
double r = dianji(l.b - l.a, c - l.a) / (L * L);
if (zhengfu(r) == -1) {
return changdu(c - l.a);
} else if (dayu_dengyu(r, 1)) {
return changdu(c - l.b);
} else {
double L = r * changdu(l.b - l.a), l2 = changdu(c - l.a);
return sqrt(l2 * l2 - L * L);
}
}
double xianduanjuli(Segment l1,Segment l2) //两线段距离
{
if (xianduan_xiangjiao(l1,l2))
return 0.0;
double minn = inf;
double l = dian_dao_xianduan(l1, l2.a);
minn = dayu_dengyu(minn, l) ? l : minn;
l = dian_dao_xianduan(l1, l2.b);
minn = dayu_dengyu(minn, l) ? l : minn;
l = dian_dao_xianduan(l2, l1.a);
minn = dayu_dengyu(minn, l) ? l : minn;
l = dian_dao_xianduan(l2, l1.b);
minn = dayu_dengyu(minn, l) ? l : minn;
return minn;
}
{%endfold%}
计算多边形面积的方法蛮多的.
最暴力的当属以原点和多边形临近两点构成三角形,然后计算三角形的有向面积.
多边形内外符号不同,最后留下的就是多边形面积,然后fabs一下就完事了.这个地方建议"脑洞大开"或者拿纸画画.
不过要是会三角剖分的话,把多边形按顶点分割成一堆三角形,然后求面积也阔以.
{%fold 代码内容 %}
double duobianxingmianji(int n) //多边形面积
{
double ans = 0;
for (int i = 1; i <= n; i++)
ans += chaji(p[i], p[(i + 1) % n]);
ans = fabs(ans) * 0.5;
return ans;
}
{%endfold%}
问所给的多边形是不是凸的.
题目给的方法是计算内角和.
嘛,感觉好麻烦的亚子.
还不如暴力求个凸包,看看所给的多边形的点数是不是和凸包点数相同来的快.
{%fold 代码内容 %}
void Andrew(int& tail) //求凸包
{
sort(p + 1, p + 1 + n);
tail = 1;
q[1] = p[1];
for (int i = 2; i <= n; i++) {
while (tail > 1 && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
tail--;
q[++tail] = p[i];
}
int basic = tail;
for (int i = n - 1; i >= 1; i--) {
while (tail > basic && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
tail--;
q[++tail] = p[i];
}
}
{%endfold%}
就,判断点和多边形的位置关系.
看网上都是角度和或者射线法.
结果就让我找到一个看起来很\(nb\)的象限角度法?
不用考虑角度的精度问题,还不用像射线法考虑多??
{%fold 代码内容 %}
bool zaibianshang(point& t,int n) //点在多边形边上否
{
for (int i = 1; i <=n; i++)
if (dian_zai_xianshang(Line(p[i], p[(i + 1)%n]), t))
return 1;
return 0;
}
bool duobianxingnei(point& t,int n) //点在多边形内
{
int t1, t2, sum, i;
double f;
p[0] = p[n];
for (i = 0; i <= n; i++)
p[i] = p[i] - t; // 坐标平移
t1 = p[0].x >= 0 ? (p[0].y >= 0 ? 0 : 3) : (p[0].y >= 0 ? 1 : 2); // 计算象限
for (sum = 0, i = 1; i <= n; i++) {
if (fabs(p[i].x)<eps && fabs(p[i].y)<eps)
break;
// 被测点为多边形顶点
f = chaji(p[i - 1], p[i - 1]);
// 计算叉积
if (fabs(f)<eps && p[i - 1].x * p[i].x <= 0 && p[i - 1].y * p[i].y <= 0)
break; // 点在边上
t2 = p[i].x >= 0 ? (p[i].y >= 0 ? 0 : 3) : (p[i].y >= 0 ? 1 : 2); // 计算象限
if (t2 == (t1 + 1) % 4)
sum += 1;
// 情况1
else if (t2 == (t1 + 3) % 4)
sum -= 1;
// 情况2
else if (t2 == (t1 + 2) % 4)
// 情况3
if (f > 0)
sum += 2;
else
sum -= 2;
t1 = t2;
}
bool tf = 0;
if (i <= n || sum)
tf = 1;
for (i = 0; i <= n; i++)
p[i] = p[i] + t; // 恢复坐标
return tf;
}
{%endfold%}
就是...求凸包
这里有个蛋疼的地方,要求是先排y再排x
{%fold 代码内容 %}
void Andrew(int& tail)
{
sort(p + 1, p + 1 + n);
tail = 1;
q[1] = p[1];
for (int i = 2; i <= n; i++) {
while (tail > 1 && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
tail--;
q[++tail] = p[i];
}
int basic = tail;
for (int i = n - 1; i >= 1; i--) {
while (tail > basic && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
tail--;
q[++tail] = p[i];
}
}
{%endfold%}
找到凸包距离最远的一对点.
就是旋转卡壳嘛.\(O(n)\)复杂度嘛.
{%fold 代码内容 %}
double tubao_zhijing(int tail) //求出凸包直径
{
double re = 0;
if (tail == 2) //仅有两个点
return changdu(q[2] - q[1]);
q[0] = q[tail]; //把最后的点放到最前面
for (int i = 0, j = 2; i < tail; i++) //枚举边
{
while (xuanzhuan(q[i], q[i + 1], q[j]) < xuanzhuan(q[i], q[i + 1], q[j + 1]))
j = (j + 1) % tail;
re = max(re, max(changdu(q[j] - q[i]), changdu(q[j] - q[i + 1])));
}
return re;
}
{%endfold%}
用一条直线切割凸包,输出得到图形的坐标.
就是逆时针找交点,按照直线的方向,\(\overrightarrow{p_1p_2}\),先放入靠近\(p_2\)的点,然后按照叉积,向左旋转的放入点.
{%fold 代码内容 %}
int zhixian_xianduan_xiangjiao(Line l1, Segment l2) //直线与线段是否有交点
{
double c1 = chaji(l1.b - l1.a, l2.a - l1.a), d1 = chaji(l1.b - l1.a, l2.b - l1.a);
if (zhengfu(c1) == 0 && zhengfu(d1) == 0) { //重合
return -1;
} else if (zhengfu(c1 * d1) <= 0) //有交点
return 1;
return 0; //平行
}
void qiegetubao(Line l)
{
int tail = 0;
q[0] = q[n];
int tf = -1, x;
for (int i = 0; i < n; i++) {
if (zhengfu(xuanzhuan(l.a,l.b,q[i])) == 1)
p[++tail] = q[i];
int f = zhixian_xianduan_xiangjiao(l,Segment( q[i], q[i + 1]));
if (f == 1) {
tf = 1;
p[++tail] = zhixian_xianduan_jiaodian(l, Segment(q[i], q[i + 1]));
} else if (f == -1) {
tf = 0;
x = i;
break;
}
}
if (tf == 0) {
//cout<<q[x+1]<<‘ ‘<<q[x]<<‘ ‘<<endl;
if (zhengfu(dianji(l.b - l.a, q[x + 1] - q[x])) == -1)
printf("%.8f\n", 0.0);
else {
for (int i = 1; i <= n; i++)
p[i] = q[i];
printf("%.8f\n", duobianxingmianji(n));
}
} else {
/*for (int i = 1; i <= tail; i++)
cout << p[i] << endl;*/
printf("%.8f\n", duobianxingmianji(tail));
}
q[0] = point(0, 0);
}
{%endfold%}
在空间内找到最近的点对
分治,建议到洛谷上搜索"平面最近点对(加强版)"
{%fold 代码内容 %}
double solve(int l, int r)
{
if (l == r)
return inf;
int mid = (l + r) >> 1;
double ans = solve(l, mid);
ans = min(ans, solve(mid + 1, r));
int tot = 0;
for (int i = l; i <= r; i++)
if (fabs(p[mid].x - p[i].x) <= ans)
temp[tot++] = p[i];
sort(temp,temp+tot, cmp);
for (int i = 0; i < tot; i++)
for (int j = i + 1; j <tot; j++) {
if (temp[j].y - temp[i].y > ans)
break;
ans = min(ans, changdu(temp[j] - temp[i]));
}
return ans;
}
{%endfold%}
扫描线算法,例题可在洛谷上搜索"扫描线"
求平面\(n\)条线段的交点个数.
有空单开一章来整理这个算法.
{%fold 代码内容 %}
struct SegTree {
int l, r;
ll len;
// sum: 被完全覆盖的次数；
// len: 区间内被截的长度。
} tree[MAXN << 2];
void build_tree(int x, int l, int r)
{
tree[x].l = l, tree[x].r = r;
tree[x].len = 0;
if (l == r)
return;
int mid = (l + r) >> 1;
build_tree(lson(x), l, mid);
build_tree(rson(x), mid + 1, r);
return;
}
void pushup(int x)
{
tree[x].len = tree[x << 1].len + tree[x << 1 | 1].len;
// 合并儿子信息
}
void edit_tree(int x, ll id, int c)
{
int l = tree[x].l, r = tree[x].r;
if (l == id && id == r) {
tree[x].len += c;
return;
}
int mid = (l + r) >> 1;
if (id <= mid)
edit_tree(lson(x), id, c);
else
edit_tree(rson(x), id, c);
pushup(x);
}
int query(int x, int L, int R)
{
int l = tree[x].l, r = tree[x].r;
if (L <= l && r <= R)
return tree[x].len;
int mid = (l + r) >> 1;
int res = 0;
if (L <= mid)
res += query(lson(x), L, R);
if (R > mid)
res += query(rson(x), L, R);
return res;
}
vector<Line> l;
vector<int> X;
int main()
{
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
point a, b;
cin >> a >> b;
if (a.x != b.x) {
if (a.x > b.x)
swap(a, b);
X.push_back(a.x), X.push_back(b.x);
l.push_back(Line(a, b, 2));
} else {
if (a.y > b.y)
swap(a, b);
X.push_back(a.x);
l.push_back(Line(a, a, 1));
l.push_back(Line(b, b, 3));
}
}
sort(l.begin(), l.end());
sort(X.begin(), X.end());
n = unique(X.begin(), X.end()) - X.begin(); //去重
X.erase(X.begin() + n, X.end());
if(n==1){
cout<<0<<endl;
return 0;
}
build_tree(1, 1, n); //根据X[]的坐标建立线段树
int ans = 0;
for (int i = 0; i < l.size(); i++) {
if (l[i].mark!=2) {
int id = lower_bound(X.begin(), X.end(), l[i].a.x) - X.begin();
int c = l[i].mark == 3 ? -1 : 1;
edit_tree(1, id+1, c);
} else {
int id1 = lower_bound(X.begin(), X.end(), l[i].a.x) - X.begin();
int id2 = lower_bound(X.begin(), X.end(), l[i].b.x) - X.begin();
ans += query(1, id1+1, id2+1);
}
}
printf("%d\n", ans);
return 0;
}
{%endfold%}
求圆的切线个数
{%fold 代码内容 %}
int yuan_yuan_xiangjiao(cirles a, cirles b) //询问圆和圆的切线个数
{
double l = changdu(b.o - a.o);
if (bijiao(l, a.r + b.r) == 1)
return 4;
else if (bijiao(l, a.r + b.r) == 0)
return 3;
else {
if (bijiao(l + min(a.r, b.r), max(a.r, b.r)) == 1)
return 2;
else if (bijiao(l + min(a.r, b.r), max(a.r, b.r)) == 0)
return 1;
else
return 0;
}
}
{%endfold%}
求圆和直线的交点.
建议阅读"挑战程序设计竞赛2"
{%fold 代码内容 %}
pair<point, point> zhixian_yuan_jiaodian(Line l, cirles a) //求直线与圆交点
{
if (dayu_dengyu(a.r, zhixian_yuanxin_juli(l, a))) {
vec a_i = touying(l, a.o);
double L = changdu(l.b - l.a);
double lt = changdu(a.o - a_i);
double lr = sqrt(a.r * a.r - lt * lt);
vec p = (l.b - l.a) / L * lr + a_i;
return make_pair(a_i - (l.b - l.a) / L * lr, a_i + (l.b - l.a) / L * lr);
} else
return make_pair(0, 0);
}
{%endfold%}
求两个圆交点
建议阅读"挑战程序设计竞赛2"
{%fold 代码内容 %}
pair<vec, vec> yuan_yuan_jiaodian(cirles a, cirles b) //求圆和圆的交点
{
double l = changdu(a.o - b.o);
double x = acos((a.r * a.r + l * l - b.r * b.r) / (2.0 * a.r * l));
double t = atan2((b.o - a.o).y, (b.o - a.o).x);
return make_pair(a.o + vec(cos(t - x) * a.r, sin(t - x) * a.r), a.o + vec(cos(x + t) * a.r, sin(t + x) * a.r));
}
{%endfold%}
过一个点做圆的切线
根据半径和圆心到点的距离求出夹角,旋转角度,得到交点
{%fold 代码内容 %}
pair<point, point> guodian_yuan_qiedian(cirles a, point p) //过一点做圆的切线求切点
{
double l = changdu(a.o - p);
double t = asin(a.r / l);
double lb = l * cos(t);
vec x = a.o - p;
x = x / l * lb;
return make_pair(p + xiangliang_xuanzhuan(x, t), p + xiangliang_xuanzhuan(x, -t));
}
{%endfold%}
求两个圆的公切线
根据圆和圆的位置进行判断
{%fold 代码内容 %}
int yuanyuan_gongqiexian(cirles a, cirles b, point* u, point* v) //求圆和圆公切线以及切线个数
{
int cnt = 0;
if (a.r < b.r) {
swap(a, b);
swap(u, v);
}
double l = changdu(a.o - b.o);
double rdiff = a.r - b.r;
double rsum = a.r + b.r;
if (zhengfu(l - rdiff) < 0)
return 0;
double base = atan2((b.o - a.o).y, (b.o - a.o).x);
if (zhengfu(l) == 0)
return -1;
if (zhengfu(l - rdiff) == 0) {
u[cnt] = v[cnt] = a.Point(base);
cnt++;
return 1;
}
double ang = acos((a.r - b.r) / l);
u[cnt] = a.Point(base + ang);
v[cnt] = b.Point(base + ang);
cnt++;
u[cnt] = a.Point(base - ang);
v[cnt] = b.Point(base - ang);
cnt++;
if (zhengfu(l - rsum) == 0) {
u[cnt] = v[cnt] = a.Point(base);
cnt++;
} else if (zhengfu(l - rsum) > 0) {
double ang = acos((a.r + b.r) / l);
u[cnt] = a.Point(base + ang);
v[cnt] = b.Point(pi + base + ang);
cnt++;
u[cnt] = a.Point(base - ang);
v[cnt] = b.Point(pi + base - ang);
cnt++;
}
return cnt;
}
{%endfold%}
求多边形和圆相交的面积
将多边形的边的顶点与圆心连接行成三角形.
那么面积便是三角形在圆内的有向面积.
对每个三角形在圆内进行判断来计算.
{%fold 代码内容 %}
point zhixian_zhixian_jiaodian(Line l1, Line l2) //两直线交点
{
double t = ((l1.a.x - l2.a.x) * (l2.a.y - l2.b.y) - (l1.a.y - l2.a.y) * (l2.a.x - l2.b.x)) / ((l1.a.x - l1.b.x) * (l2.a.y - l2.b.y) - (l1.a.y - l1.b.y) * (l2.a.x - l2.b.x));
return l1.a + (l1.b - l1.a) * t;
}
point xianduan_duandian_dian(point p, Segment l) //线段距离点p最近的端点
{
point t = p;
t.x += l.a.y - l.b.y;
t.y += l.b.x - l.a.x;
if (chaji(l.a - p, t - p) * chaji(l.b - p, t - p) > eps)
return changdu(p - l.a) < changdu(p - l.b) ? l.a : l.b;
return zhixian_zhixian_jiaodian(Line(p, t), l);
}
double distp(Line l) //长度的平方
{
return (l.a.x - l.b.x) * (l.a.x - l.b.x) + (l.a.y - l.b.y) * (l.a.y - l.b.y);
}
double yuanxin_dian_sanjiao(Line l, cirles c) //求圆心与两点所成三角形有向面积
{
double sign = 1.0;
l.a = l.a - c.o;
l.b = l.b - c.o;
c.o = point(0.0, 0.0);
if (fabs(chaji(l.a - c.o, l.b - c.o)) < eps)
return 0.0;
if (distp(Line(l.a, c.o)) > distp(Line(l.b, c.o))) {
swap(l.a, l.b);
sign = -1.0;
}
if (distp(Line(l.a, c.o)) < c.r * c.r + eps) { //a在圆内
if (distp(Line(l.b, c.o)) < c.r * c.r + eps) //b也在圆内,返回叉积/2
return chaji(l.a - c.o, l.b - c.o) / 2.0 * sign;
point p1, p2;
pair<point, point> q = zhixian_yuan_jiaodian(l, c); //oa和ob与圆的交点
p1 = q.first;
p2 = q.second;
if (changdu(p1 - l.b) > changdu(p2 - l.b))
swap(p1, p2);
double ret1 = fabs(chaji(l.a - c.o, p1 - c.o));
double ret2 = acos((p1.x * l.b.x + p1.y * l.b.y) / changdu(p1) / changdu(l.b)) * c.r * c.r;
double ret = (ret1 + ret2) / 2.0;
if (chaji(l.a - c.o, l.b - c.o) < eps && sign > 0.0 || chaji(l.a - c.o, l.b - c.o) > eps && sign < 0.0)
ret = -ret;
return ret;
}
point ins = xianduan_duandian_dian(c.o, l);
if (distp(Line(c.o, ins)) > c.r * c.r - eps) {
double ret = acos((l.a.x * l.b.x + l.a.y * l.b.y) / changdu(l.a) / changdu(l.b)) * c.r * c.r / 2.0;
if (chaji(l.a - c.o, l.b - c.o) < eps && sign > 0.0 || chaji(l.a - c.o, l.b - c.o) > eps && sign < 0.0)
ret = -ret;
return ret;
}
point p1, p2;
pair<point, point> q = zhixian_yuan_jiaodian(l, c); //oa和ob与圆的交点
p1 = q.first;
p2 = q.second;
double cm = c.r / (changdu(c.o - l.a) - c.r);
point m = point((c.o.x + cm * l.a.x) / (1 + cm), (c.o.y + cm * l.a.y) / (1 + cm));
double cn = c.r / (changdu(c.o - l.b) - c.r);
point n = point((c.o.x + cn * l.b.x) / (1 + cn), (c.o.y + cn * l.b.y) / (1 + cn));
double ret1 = acos((m.x * n.x + m.y * n.y) / changdu(m) / changdu(n)) * c.r * c.r;
double ret2 = acos((p1.x * p2.x + p1.y * p2.y) / changdu(p1) / changdu(p2)) * c.r * c.r - fabs(chaji(p1 - c.o, p2 - c.o));
double ret = (ret1 - ret2) / 2.0;
if (chaji(l.a - c.o, l.b - c.o) < eps && sign > 0.0 || chaji(l.a - c.o, l.b - c.o) > eps && sign < 0.0)
ret = -ret;
return ret;
}
double duobianxing_yuan_xiangjiao(cirles c, point p[], int n) //多边形与圆相交面积
{
double sum = 0;
for (int i = 0; i < n; i++)
sum += yuanxin_dian_sanjiao(Line(p[i], p[i + 1]), c);
return sum;
}
{%endfold%}
证明待补
{% fold 代码内容 %}
double fun(double x)//原积分函数
{
return x;
}
double simpson(double l,double r) //辛普森法则
{
double mid=(l+r)/2.0;
return (fun(l)+fun(r)+4*fun(mid))*(r-l)/6.0;
}
double solve(double l,double r,double ans) //调整精度求答案
{
double mid=(l+r)/2.0;
double ls=simpson(l,mid),rs=simpson(mid,r);
if(fabs(ls+rs-ans)<=15.0*eps)
return ls+rs+(ls+rs-ans)/15.0;
return solve(l,mid,ls)+solve(mid,r,rs);
}
{%endfold%}
{% fold 代码内容 %}
point zhixian_zhixian_jiaodian(Line l1, Line l2) //两直线交点
{
double t = ((l1.a.x - l2.a.x) * (l2.a.y - l2.b.y) - (l1.a.y - l2.a.y) * (l2.a.x - l2.b.x)) / ((l1.a.x - l1.b.x) * (l2.a.y - l2.b.y) - (l1.a.y - l1.b.y) * (l2.a.x - l2.b.x));
return l1.a + (l1.b - l1.a) * t;
}
point sanjiaoxing_waixin(point a, point b, point c) //三角形外心
{
Line u, v;
u.a.x = (a.x + b.x) / 2;
u.a.y = (a.y + b.y) / 2;
u.b.x = u.a.x - a.y + b.y;
u.b.y = u.a.y + a.x - b.x;
v.a.x = (a.x + c.x) / 2;
v.a.y = (a.y + c.y) / 2;
v.b.x = v.a.x - a.y + c.y;
v.b.y = v.a.y + a.x - c.x;
return zhixian_zhixian_jiaodian(u, v);
}
cirles zuixiaoyuan_fugai(point p[], int n) //最小圆覆盖
{
random_shuffle(p + 1, p + 1 + n);
cirles c(p[1], 0.0);
for (int i = 2; i <= n; i++)
if (zhengfu(changdu(c.o - p[i]) - c.r) > 0) {
c = cirles(p[i], 0.0);
for (int j = 1; j < i; j++)
if (zhengfu(changdu(c.o - p[j]) - c.r) > 0) {
c.o = (p[i] + p[j]) / 2;
c.r = changdu(c.o - p[i]);
for (int k = 1; k < j; k++)
if (zhengfu(changdu(c.o - p[k]) - c.r) > 0) {
c.o = sanjiaoxing_waixin(p[i], p[j], p[k]);
c.r = changdu(c.o - p[i]);
}
}
}
return c;
}
{%endfold%}
原文：https://www.cnblogs.com/CrossAutomaton/p/14379616.html