lrj计算几何模板

整理了一下大白书上的计算几何模板。

  1 #include <cstdio>
  2 #include <algorithm>
  3 #include <cmath>
  4 using namespace std;
  5 //lrj计算几何模板
  6 struct Point
  7 {
  8     double x, y;
  9     Point(double x=0, double y=0) :x(x),y(y) {}
 10 };
 11 typedef Point Vector;
 12 const double EPS = 1e-10;
 13
 14 //向量+向量=向量 点+向量=点
 15 Vector operator + (Vector A, Vector B)    { return Vector(A.x + B.x, A.y + B.y); }
 16
 17 //向量-向量=向量 点-点=向量
 18 Vector operator - (Vector A, Vector B)    { return Vector(A.x - B.x, A.y - B.y); }
 19
 20 //向量*数=向量
 21 Vector operator * (Vector A, double p)    { return Vector(A.x*p, A.y*p); }
 22
 23 //向量/数=向量
 24 Vector operator / (Vector A, double p)    { return Vector(A.x/p, A.y/p); }
 25
 26 bool operator < (const Point& a, const Point& b)
 27 { return a.x < b.x || (a.x == b.x && a.y < b.y); }
 28
 29 int dcmp(double x)
 30 { if(fabs(x) < EPS) return 0; else x < 0 ? -1 : 1; }
 31
 32 bool operator == (const Point& a, const Point& b)
 33 { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; }
 34
 35 /**********************基本运算**********************/
 36
 37 //点积
 38 double Dot(Vector A, Vector B)
 39 { return A.x*B.x + A.y*B.y; }
 40 //向量的模
 41 double Length(Vector A)    { return sqrt(Dot(A, A)); }
 42
 43 //向量的夹角，返回值为弧度
 44 double Angle(Vector A, Vector B)
 45 { return acos(Dot(A, B) / Length(A) / Length(B)); }
 46
 47 //叉积
 48 double Cross(Vector A, Vector B)
 49 { return A.x*B.y - A.y*B.x; }
 50
 51 //向量AB叉乘AC的有向面积
 52 double Area2(Point A, Point B, Point C)
 53 { return Cross(B-A, C-A); }
 54
 55 //向量A旋转rad弧度
 56 Vector VRotate(Vector A, double rad)
 57 {
 58     return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));
 59 }
 60
 61 //将B点绕A点旋转rad弧度
 62 Point PRotate(Point A, Point B, double rad)
 63 {
 64     return A + VRotate(B-A, rad);
 65 }
 66
 67 //求向量A向左旋转90°的单位法向量,调用前确保A不是零向量
 68 Vector Normal(Vector A)
 69 {
 70     double l = Length(A);
 71     return Vector(-A.y/l, A.x/l);
 72 }
 73
 74 /**********************点和直线**********************/
 75
 76 //求直线P + tv 和 Q + tw的交点，调用前要确保两条直线有唯一交点
 77 Point GetLineIntersection(Point P, Vector v, Point Q, Vector w)
 78 {
 79     Vector u = P - Q;
 80     double t = Cross(w, u) / Cross(v, w);
 81     return P + v*t;
 82 }//在精度要求极高的情况下，可以自定义分数类
 83
 84 //P点到直线AB的距离
 85 double DistanceToLine(Point P, Point A, Point B)
 86 {
 87     Vector v1 = B - A, v2 = P - A;
 88     return fabs(Cross(v1, v2)) / Length(v1);    //不加绝对值是有向距离
 89 }
 90
 91 //点到线段的距离
 92 double DistanceToSegment(Point P, Point A, Point B)
 93 {
 94     if(A == B)    return Length(P - A);
 95     Vector v1 = B - A, v2 = P - A, v3 = P - B;
 96     if(dcmp(Dot(v1, v2)) < 0)    return Length(v2);
 97     else if(dcmp(Dot(v1, v3)) > 0)    return Length(v3);
 98     else return fabs(Cross(v1, v2)) / Length(v1);
 99 }
100
101 //点在直线上的射影
102 Point GetLineProjection(Point P, Point A, Point B)
103 {
104     Vector v = B - A;
105     return A + v * (Dot(v, P - A) / Dot(v, v));
106 }
107
108 //线段“规范”相交判定
109 bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2)
110 {
111     double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);
112     double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
113     return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
114 }
115
116 //判断点是否在线段上
117 bool OnSegment(Point P, Point a1, Point a2)
118 {
119     Vector v1 = a1 - P, v2 = a2 - P;
120     return dcmp(Cross(v1, v2)) == 0 && dcmp(Dot(v1, v2)) < 0;
121 }
122
123 //求多边形面积
124 double PolygonArea(Point* P, int n)
125 {
126     double ans = 0.0;
127     for(int i = 1; i < n - 1; ++i)
128         ans += Cross(P[i]-P[0], P[i+1]-P[0]);
129     return ans/2;
130 }