zoukankan      html  css  js  c++  java
  • 计算几何模板(未完待续)

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

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

    二维基础操作:

     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 }
  • 相关阅读:
    ASP.NET 3.5 的 ListView 控件与 CSS Friendly
    从 Adobe SHARE 说到 Silverlight 的 XPS 支持
    编写 iPhone Friendly 的 Web 应用程序 (Part 5 交互入门)
    初为项目经理
    管理的最高境界不是完美
    url传递中文的解决方案总结
    我想跟什么样的人合作
    异步Socket通信总结
    让机器自动支持中文文件名
    Socket基本编程
  • 原文地址:https://www.cnblogs.com/AlphaWA/p/10201050.html
Copyright © 2011-2022 走看看