计算几何模板 | Kyons'Blog 




 
 
 
 



 
 



计算几何模板

计算几何模板

正如不知何方大佬所言,计算几何精妙之处,就是不用解析几何的方法去做
为了方便查找,防止自己迷路,我把函数名都写成了拼音
绝对不是因为我英语不好!!!

基本数据结构

点和向量:

  • 点和向量都可以用一个坐标$(x,y)$来表示.
  • 故向量$Vector$可以写为
  • 1
    typedef struct point vec;

完整定义如下

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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;
}
};

直线和线段:

  • 直线和线段都可以用两个点的坐标来表示
  • 故线段$Segment$可以写为
  • 1
    typedef struct Line Segment;

完整定义如下

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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;
}
};

圆:

  • 圆的表示有一个点圆心,以及其半径组成

完整定义如下

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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;
}
};

基本函数以及常量

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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);
}

AOJ相关习题

CGL_1_A:Projection

求一个点在向量$\overrightarrow{ab}$上的投影坐标
设点$c$,投影在$\overrightarrow{ab}$上为$c’$,则$c’$的坐标就是:$cos<\overrightarrow{ac},\overrightarrow{ab}>\times |\overrightarrow{ac}|+a$

代码内容
1
2
3
4
5
6
7
8
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;
}

CGL_1_B:Reflection

求一个点$c$关于向量$\overrightarrow{ab}$的对称点$c’’$
先求出$c$在$ab$上的投影,那么$c’’=2\times \overrightarrow{cc’}+c$

代码内容
1
2
3
4
5
point fanshe(Line l, point c) //求c关于直线ab的对称点c'
{
point A = touying(l, c);
return (A - c) * 2.0 + c;
}

CGL_1_C:Counter-Clockwise

就是..根据图中的判断就是了

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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;
}
}
}

CGL_2_A:Parallel/Orthogonal

先判平行,再用点积判垂直

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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;
}
}
}

CGL_2_B:Intersection

线段相交要考虑蛮多的,首先,先按照x后y从小到大排一下.
最简单的情况,$ab$穿过$cd$,那么必定有交点.
第二种,$a$在$cd$上或者$b$在$cd$上
第三种,共线时,$a$在$cd$之间或$b$在$cd$之间.
处理好以上问题,就解决了

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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;
}

CGL_2_C:Cross Point

给你两个必定相交的线段,求交点

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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);
}

CGL_2_D:Distance

  • 给定两个不相交线段,求两个线段最近距离
  • 很明显,最近距离就是两个端点到另一个线段的距离.
  • 那么两遍点到线段距离就出来了.
    • 点到线段距离有三种
    • 第一种是点在线段正上方,则距离为过点向线段作垂线
    • 第二种是点在左侧,就是左端点和该点连线
    • 第三种同第二种,不过在右侧
代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
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;
}

CGL_3_A:Area

计算多边形面积的方法蛮多的.
最暴力的当属以原点和多边形临近两点构成三角形,然后计算三角形的有向面积.
多边形内外符号不同,最后留下的就是多边形面积,然后fabs一下就完事了.这个地方建议”脑洞大开”或者拿纸画画.
不过要是会三角剖分的话,把多边形按顶点分割成一堆三角形,然后求面积也阔以.

代码内容
1
2
3
4
5
6
7
8
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;
}

CGL_3_B:Is-Convex

问所给的多边形是不是凸的.
题目给的方法是计算内角和.
嘛,感觉好麻烦的亚子.
还不如暴力求个凸包,看看所给的多边形的点数是不是和凸包点数相同来的快.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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];
}
}

CGL_3_C:Polygon-Point Containment

就,判断点和多边形的位置关系.
看网上都是角度和或者射线法.
结果就让我找到一个看起来很$nb$的象限角度法?
不用考虑角度的精度问题,还不用像射线法考虑多??

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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;
}

CGL_4_A:Convex Hull

就是…求凸包
这里有个蛋疼的地方,要求是先排y再排x

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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];
}
}

CGL_4_B:Diameter of a Convex Polygon

找到凸包距离最远的一对点.
就是旋转卡壳嘛.$O(n)$复杂度嘛.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
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;
}

CGL_4_C:Convex Cut

用一条直线切割凸包,输出得到图形的坐标.
就是逆时针找交点,按照直线的方向,$\overrightarrow{p_1p_2}$,先放入靠近$p_2$的点,然后按照叉积,向左旋转的放入点.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
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);
}

CGL_5_A:Closest Pair

在空间内找到最近的点对
分治,建议到洛谷上搜索”平面最近点对(加强版)”

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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;
}

CGL_6_A:Segment Intersections: Manhattan Geometry

扫描线算法,例题可在洛谷上搜索”扫描线”
求平面$n$条线段的交点个数.
有空单开一章来整理这个算法.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
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;
}

CGL_7_A:Intersection

求圆的切线个数

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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;
}
}

CGL_7_D:Cross Points of a Circle and a Line

求圆和直线的交点.
建议阅读”挑战程序设计竞赛2”

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
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);
}

CGL_7_E:Cross Points of Circles

求两个圆交点
建议阅读”挑战程序设计竞赛2”

代码内容
1
2
3
4
5
6
7
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));
}

CGL_7_F:Tangent to a Circle

过一个点做圆的切线
根据半径和圆心到点的距离求出夹角,旋转角度,得到交点

代码内容
1
2
3
4
5
6
7
8
9
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));
}

CGL_7_G:Common Tangent

求两个圆的公切线
根据圆和圆的位置进行判断

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
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;
}

CGL_7_H:Intersection of a Circle and a Polygon

求多边形和圆相交的面积
将多边形的边的顶点与圆心连接行成三角形.
那么面积便是三角形在圆内的有向面积.
对每个三角形在圆内进行判断来计算.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
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;
}

其他相关

辛普森法则

证明待补

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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);
}

最小圆覆盖

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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;
}

 


 
 
 





 

 


 


 

 

 

 

 
 

 

 

 //删除Valine核心代码库外链调用 //调用刚下载的本地文件以加速加载速度 //这里改为从本地加载