计算几何模板（未完待续）

Posted 2021-01-31 alphawa

目前基本都是从蓝书上摘录的。

有一部分需要线性代数的知识，但是蓝书作者并没有解释，个人觉得用数学知识推出来更有助于记忆，死记硬背板子容易忘。以后有机会的话我在这里写点注解。

二维基础操作：

 1 struct Point
 2 {
 3     double x, y;
 4     Point(double x = 0, double y = 0):x(x), y(y){}
 5 };
 6 
 7 typedef Point Vector;
 8 
 9 int dcmp(double x)//比较
10 {
11     const double eps = 1e-10;
12     if (fabs(x) < eps)    return 0;
13     else    return x < 0 ? -1 : 1;
14 }
15 Vector operator + (Vector A, Vector B) { return Vector(A.x + B.x, A.y + B.y); }
16 Vector operator - (Vector A, Vector B) { return Vector(A.x - B.x, A.y - B.y); }
17 Vector operator * (Vector A, double p) { return Vector(A.x * p, A.y * p); }
18 Vector operator / (Vector A, double p) { return Vector(A.x / p, A.y / p); }
19 bool operator < (const Point& A, const Point& B) { return A.x < B.x || (A.x == B.x && A.y < B. y); }
20 bool operator == (const Point& A, const Point& B) { return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.y) == 0; }
21 
22 const double PI = acos(-1.0);
23 double torad(double deg) { return deg/180 * PI; }//角度转弧度
24 
25 double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y; }//点积
26 
27 double Length(Vector A) { return sqrt(Dot(A, A)); }//长度
28 
29 double Angle(Vector A, Vector B) { return acos(Dot(A, B)/Length(A)/Length(B)); }//角度
30 
31 Vector Rotate(Vector A, double rad) { return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad)); }//旋转
32 
33 double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; }//叉积
34 
35 double Area2(Point A, Point B, Point C) { return Cross(B-A, C-A);}//二倍的三角形面积
36 
37 Vector Normal(Vector A) { double L = Length(A); return Vector(-A.y/L, A.x/L); }//向量的单位法线，旋转的特殊情况
38 
39 Vector Get_Line_Projection(Point P, Point A, Point B)//P在AB方向上的映射
40 {
41      Vector v = B - A;
42       return A + v*(Dot(v, P-A) / Dot(v, v));
43 }
44 
45 Vector Get_Line_Intersect(Vector A, Vector v, Vector B, Vector w)//两直线交点
46 {
47     Vector u = A - B;
48     double t = Cross(w, u) / Cross(v, w);
49     return A + v*t;
50 }
51 
52 double Distance_to_Line(Point P, Point A, Point B)//点P到直线AB的距离
53 {
54     Vector v1 = B - A, v2 = P - A;
55     return fabs(Cross(v1, v2)) / Length(v1);
56 }
57 
58 double Distance_to_Segment(Point P, Point A, Point B)//点P到线段AB的距离
59 {
60     Vector v1 = B - A, v2 = P - A, v3 = P - B;
61     if (A == B || dcmp(Dot(v1, v2)) < 0)    return Length(v2);//PA
62     else if (dcmp(Dot(v1, v3)) > 0)    return Length(v3);//PB
63     else    return fabs(Cross(v1, v2) / Length(v1));//PQ
64 }
65 
66 bool On_Segment(Point P, Point A, Point B) { return dcmp(Cross(A - P, B - P)) == 0 && dcmp(Dot(A - P, B - P)) < 0; }//点在线段上
67 
68 bool Segment_Proper_Intersection(Point a1, Point a2, Point b1, Point b2)//线段a1a2与b1b2相交（不包含端点）
69 {
70     double c1 = Cross(a2 - a1, b1 - a1), c2 = Cross(a2 - a1, b2 - a1);
71     double c3 = Cross(b2 - b1, a1 - b1), c2 = Cross(b2 - b1, a2 - b1);
72     return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0;
73 }
74 
75 double Polygon_Area(Point *p, int n)//多边形面积，逆时针取点0~n-1
76 {
77     double area = 0;
78     for (int i = 1; i < n-1; i++)
79         area += Cross(p[i] - p[0], p[i+1] - p[0]);
80     return area/2;
81 }
82 
83 //基于复数的几何运算，效率不是很高
84 #include <complex>
85 typedef complex<double> Point;
86 typedef Point Vector;
87 
88 double Dot(Vector A, Vector B) { return real(conj(A) * B); }// A（x+yi）的共轭（x-yi）乘以B，取实部
89 double Cross(Vector A, Vector B) { return imag(conj(A) * B); }//取虚部
90 Vector Rotate(Vector A, double rad) { return A * exp(Point(0, rad)); }//使用了欧拉公式：e^(0 + rad*i) = cos(rad) + sin(rad)*i，求完确实是旋转

 

圆相关（他那个两圆公切线看着有点奇怪先不贴了）：

 1 struct Circle
 2 {
 3     Point c;
 4     double r;
 5     // Circle(Point a = (0, 0), double x = 0):c(a), r(x){}
 6     Point point(double seta) { return Point(c.x + cos(seta)*r, c.y + sin(seta)*r); }//已知角度求圆上某点
 7 };
 8 
 9 struct Line
10 {
11     Point p;
12     Vector v;
13     Line(Point p, Vector v):p(p), v(v) { }
14 
15     Point point(double t) { return p + v*t; }//已知参数t求直线上一点
16     Line move(double d) { return Line(p + Normal(v)*d, v); }//平移
17 };
18 
19 int get_Tangents(Point P, Circle C, vector<Vector> &v)//过定点做圆的切线
20 {
21     Vector u = C.c - P;
22     double d = Length(u);
23 
24     if (dcmp(d - C.r) < 0)    return 0;//点在圆内部
25     else if (dcmp(d - C.r) == 0)//点在圆上
26     {
27         v.push_back(Rotate(u, PI/2));
28         return 1;
29     }
30     else//两条切线
31     {
32         double a = asin(C.r / d);
33         v.push_back(Rotate(u, a));
34         v.push_back(Rotate(u, -a));
35         return 2;
36     }
37 }
38 
39 int get_Line_Circle_Intersection(Line L, Circle C, vector<Point> &ans)//方程求解直线与圆的交点
40 {
41     double t1, t2;
42     double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
43     double e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
44     double delta = f*f - 4*e*g;//判别式
45     
46     if (dcmp(delta) < 0)    return 0;//相离
47     else if (dcmp(delta) == 0)//相切
48     {
49         t1 = t2 = -f/2/e;
50         ans.push_back(L.point(t1));
51         return 1;
52     }
53     else//相交
54     {
55         t1 = (-f + sqrt(delta)) / 2 / e;
56         t2 = (-f - sqrt(delta)) / 2 / e;
57         ans.push_back(L.point(t1)), ans.push_back(L.point(t2));
58         return 2;
59     }
60 }
61 
62 inline double angle(Vector A) { return atan2(A.y, A.x); }//某向量极角
63 
64 int get_Circle_Circle_Intersection(Circle C1, Circle C2, vector<Point> &v)//圆与圆的交点
65 {
66     double d = Length(C1.c - C2.c);
67     if (dcmp(d) == 0)
68     {
69         if (dcmp(C1.r - C2.r) == 0)    return -1;//重合
70         return 0;//同心
71     }
72     if (dcmp(C1.r + C2.r - d) < 0)    return 0;//外离
73     if (dcmp(fabs(C1.r - C2.r) - d) > 0)    return 0;//内含
74 
75     double a = angle(C2.c - C1.c);
76     double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d));//余弦公式
77     Point p1 = C1.point(a - da), p2 = C1.point(a + da);
78 
79     v.push_back(p1);
80     if (p1 == p2)    return 1;
81     v.push_back(p2);
82     return 2;
83 }