• 计算几何模板


    1、典型的Point结构体

     1 struct point {
     2     double x, y;
     3     point(double _x = 0, double _y = 0): x(_x), y(_y) {
     4     }
     5     void input() {
     6         scanf("%lf%lf", &x, &y);
     7     }
     8     double len() const {
     9         return sqrt(x * x + y * y);
    10     }
    11     point trunc(double l) const {
    12         double r = l / len();
    13         return point(x * r, y * r);
    14     }
    15     point rotate_left() const {
    16         return point(-y, x);
    17     }
    18     point rotate_left(double ang) const {
    19         double c = cos(ang), s = sin(ang);
    20         return point(x * c - y * s, y * c + x * s);
    21     }
    22     point rotate_right() const {
    23         return point(y, -x);
    24     }
    25     point rotate_right(double ang) const {
    26         double c = cos(ang), s = sin(ang);
    27         return point(x * c + y * s, y * c - x * s);
    28     }
    29 };
    30 
    31 bool operator==(const point& p1, const point& p2) {
    32     return sgn(p1.x - p2.x) == 0 && sgn(p1.y - p2.y) == 0;
    33 }
    34 
    35 bool operator<(const point& p1, const point& p2) {
    36     return sgn(p1.x - p2.x) == 0 ? sgn(p1.y - p2.y) < 0 : p1.x < p2.x;
    37 }
    38 
    39 point operator+(const point& p1, const point& p2) {
    40     return point(p1.x + p2.x, p1.y + p2.y);
    41 }
    42 
    43 point operator-(const point& p1, const point& p2) {
    44     return point(p1.x - p2.x, p1.y - p2.y);
    45 }
    46 
    47 double operator^(const point& p1, const point& p2) {
    48     return p1.x * p2.x + p1.y * p2.y;
    49 }
    50 
    51 double operator*(const point& p1, const point& p2) {
    52     return p1.x * p2.y - p1.y * p2.x;
    53 }

    2、转角法判断点是否在简单多边形的内

     1 int get_position(const point& p,const point* pol, int n){
     2     int wn = 0;
     3     for(int i = 0; i < n; ++i){
     4         if(onsegment(p, pol[i], pol[(i+1)%n])) return -1;   //在边界上
     5         int k = sgn((pol[(i+1)%n]-pol[i]) * (p-pol[i]));
     6         int d1 = sgn(pol[i].y - p.y);
     7         int d2 = sgn(pol[(i+1)%n].y - p.y);
     8         if(k>0 && d1<=0 && d2>0)  wn++;
     9         if(k<0 && d2<=0 && d1>0)  wn--;
    10     }
    11     if(wn != 0) return 1;     //内部
    12     return 0;               //外部
    13 }

    3、凸包

     1 int dn, hd[max_n], un, hu[max_n];
     2 
     3 void get_convex_hull(point* p, int n, point* pol, int& m) {
     4     sort(p, p + n);
     5     dn = un = 2;
     6     hd[0] = hu[0] = 0;
     7     hd[1] = hu[1] = 1;
     8     for (int i = 2; i < n; ++i) {
     9         for (; dn > 1 && sgn((p[hd[dn - 1]] - p[hd[dn - 2]]) * (p[i] - p[hd[dn - 1]])) <= 0; --dn);
    10         for (; un > 1 && sgn((p[hu[un - 1]] - p[hu[un - 2]]) * (p[i] - p[hu[un - 1]])) >= 0; --un);
    11         hd[dn++] = hu[un++] = i;
    12     }
    13     m = 0;
    14     for (int i = 0; i < dn - 1; ++i)
    15         pol[m++] = p[hd[i]];
    16     for (int i = un - 1; i > 0; --i)
    17         pol[m++] = p[hu[i]];
    18 }

    4、线段交点

    1 bool get_intersection(const point& p1, const point& p2, const point& p3, const point& p4, point& c) {
    2      double d1 = (p2 - p1) * (p3 - p1), d2 = (p2 - p1) * (p4 - p1);
    3      double d3 = (p4 - p3) * (p1 - p3), d4 = (p4 - p3) * (p2 - p3);
    4      int s1 = sgn(d1), s2 = sgn(d2), s3 = sgn(d3), s4 = sgn(d4);
    5      if (s1 == 0 && s2 == 0 && s3 == 0 && s4 == 0)
    6          return false;
    7      c = (p3 * d2 - p4 * d1) / (d2 - d1);
    8      return s1 * s2 <= 0 && s3 * s4 <= 0; 
    9 }  

    5、半平面交

     1 point dp[12000];
     2 struct half_plane {
     3         point p1, p2;
     4         double ang;
     5         half_plane(){}
     6         half_plane(const point& _p1, const point& _p2): p1(_p1), p2(_p2) {
     7                  ang = atan2(p2.y - p1.y, p2.x - p1.x);
     8                  if (sgn(ang - pi) == 0)
     9                  ang = -pi;
    10         }
    11         int get_position(const point& p) const {
    12                   return sgn((p2 - p1) * (p - p1));
    13         }
    14         bool onleft(const point& p) const {
    15                   return get_position(p) > 0;
    16         }
    17 } hp[2100], dq[1200];
    18 
    19 bool operator<(const half_plane& pl1, const half_plane& pl2) {
    20         return pl1.ang < pl2.ang;
    21 }
    22 
    23 point get_intersection(const half_plane& pl1, const half_plane& pl2) {
    24         double d1 = (pl1.p2 - pl1.p1) * (pl2.p1 - pl1.p1), 
    25                d2 = (pl1.p2 - pl1.p1) * (pl2.p2 - pl1.p1);
    26         return  point((pl2.p1.x * d2  -  pl2.p2.x * d1) / (d2 - d1), 
    27                       (pl2.p1.y * d2  -  pl2.p2.y * d1) / (d2 - d1));
    28 }
    29 
    30 void get_intersection(half_plane* pl, int n, point* pol, int& m) {
    31         m = 0;
    32         sort(pl, pl + n);
    33         int first, last;
    34         dq[first = last = 0] = pl[0];
    35         for (int i = 1; i < n; ++i){
    36                while (first < last && !pl[i].onleft(dp[last-1])) --last;
    37                while (first < last && !pl[i].onleft(dp[first])) ++first;
    38                dq[++last] = pl[i];
    39                if (sgn((dq[last].p2 - dq[last].p1) * (dq[last-1].p2 - dq[last-1].p1)) == 0){
    40                      --last;
    41                      if (dq[last].onleft(pl[i].p1)) dq[last] = pl[i];
    42                }
    43                if (first < last) dp[last-1] = get_intersection(dq[last], dq[last-1]);
    44         }
    45         while (first < last && !dq[first].onleft(dp[last-1])) --last;
    46         if (last - first <= 1) return;
    47         dp[last] = get_intersection(dq[last], dq[first]);
    48         for (int i = first; i <= last; ++i) pol[m++] = dp[i];
    49 }

    6.最小圆覆盖

     1 struct circle{
     2     point c;
     3     double r;
     4     circle(){}
     5     circle(const point& _c, double _r):c(_c), r(_r){}
     6     void out(){
     7         printf("%.2lf %.2lf %.2lf
    ", c.x, c.y, r);
     8     }
     9 } C;
    10 //三点求外接圆 
    11 circle get_circumcircle(const point& p0, const point& p1, const point& p2){
    12        double a1 = p1.x - p0.x, b1 = p1.y - p0.y, c1 = (a1*a1 + b1*b1) / 2;
    13        double a2 = p2.x - p0.x, b2 = p2.y - p0.y, c2 = (a2*a2 + b2*b2) / 2;
    14        double d = a1 * b2 - a2 * b1;
    15        point c = point(p0.x + (c1 * b2 - c2 * b1) / d, p0.y + (a1 * c2 - a2 * c1) / d);
    16        return circle(c, (p0 - c).len()); 
    17 }
    18 
    19 //求最小圆覆盖,随机化算法 
    20 circle get_mincoverCircle(point* p, const int n){
    21        random_shuffle(p, p + n);  //随机化 
    22        circle C(p[0], 0);
    23        for (int i = 1; i < n; ++i)
    24            if (sgn((C.c - p[i]).len() - C.r) > 0){
    25                 C = circle(p[i], 0);
    26                 for (int j = 0; j < i; ++j)
    27                     if (sgn((C.c - p[j]).len() - C.r) > 0){ //j点不再里面,情况B 
    28                          C = circle((p[i] + p[j]) / 2, (p[i] - p[j]).len() / 2);
    29                      //    C.out(); 
    30                          for (int k = 0; k < j; ++k) 
    31                                if (sgn((C.c - p[k]).len() - C.r) > 0){ //k点也不再里面,情况C 
    32                                    if (sgn((p[k] - p[i]) * (p[k] - p[j])) == 0){ //特判三点共线无外接圆 
    33                                       //     cout << "fuck" << endl;
    34                                            if (sgn((p[k] - p[i]).len() - (p[k] - p[j]).len()) < 0)
    35                                                 C = circle((p[k] + p[j]) / 2, (p[k] - p[j]).len() / 2);       
    36                                            else 
    37                                                 C = circle((p[i] + p[k]) / 2, (p[i] - p[k]).len() / 2);
    38                                            continue;
    39                                    }
    40                                    C = get_circumcircle(p[i], p[j], p[k]);
    41                                }
    42                     }
    43            }
    44        return C;
    45 }

    7.3维凸包

     1 struct face{
     2      int a, b, c;
     3      face(int _a = 0, int _b = 0, int _c = 0):a(_a), b(_b), c(_c){}
     4 } pol[1200];
     5 const int maxn = 510, maxf = maxn * 2;
     6 int n1, n2, pos[maxn][maxn];
     7 face buf1[maxf], buf2[maxf], *p1, *p2;
     8 int get_position(const point& p, const point& p1, const point& p2, const point& p3){
     9     return sgn((p2 - p1) * (p3 - p1) ^ (p - p1));
    10 }
    11 void check(int k, int a, int b, int s){
    12     if (pos[b][a] == 0){
    13         pos[a][b] = s;
    14         return;
    15     }
    16     if (pos[b][a] != s)
    17         p2[n2++] = (s < 0 ? face(k, b, a) : face(k, a, b));
    18     pos[b][a] = 0;
    19 }
    20 
    21 void get_convex_hull(point *p, int n, face* pol, int& m){
    22     for (int i = 1; i < n; ++i)
    23         if (p[i] != p[0]){
    24             swap(p[i], p[1]);
    25             break;
    26         }
    27     for (int i = 2; i < n; ++i){
    28        if (sgn(((p[0] - p[i]) * (p[1] - p[i])).len()) != 0){
    29             swap(p[i], p[2]);
    30             break;
    31        }
    32     }
    33     for (int i = 3; i < n; ++i){
    34         if (get_position(p[i], p[0], p[1], p[2]) != 0){
    35             swap(p[i], p[3]);
    36             break;
    37         }
    38     }
    39     p1 = buf1;
    40     p2 = buf2;
    41     n1 = n2 = 0;
    42     for (int i = 0; i < n; ++i)
    43         fill(pos[i], pos[i] + n, 0);
    44     p1[n1++] = face(0, 1, 2);
    45     p1[n1++] = face(2, 1, 0);
    46     for (int i = 3; i < n; ++i){
    47         n2 = 0;
    48         for (int j = 0; j < n1; ++j){
    49              int s = get_position(p[i], p[p1[j].a], p[p1[j].b], p[p1[j].c]);
    50              if (s == 0)
    51                   s = -1;
    52              if (s <= 0) p2[n2++] = p1[j];
    53              check(i, p1[j].a, p1[j].b, s);
    54              check(i, p1[j].b, p1[j].c, s);
    55              check(i, p1[j].c, p1[j].a, s);
    56         }
    57         swap(p1, p2);
    58         swap(n1, n2);
    59     }
    60     m = n1;
    61     copy(p1, p1 + n1, pol);
    62 }

    8.三维计算几何

     1 struct point{
     2      double x, y, z;
     3      point(double _x = 0, double _y = 0, double _z = 0):x(_x), y(_y), z(_z){}
     4      void input(){
     5         scanf("%lf%lf%lf", &x, &y, &z);
     6      }
     7      double len()const{
     8          return sqrt(x * x + y * y + z * z);
     9      }
    10      point trunc(const double &l) const{
    11            double r  = l / len();
    12            return point(x * r, y * r, z * r); 
    13      }
    14      bool operator!=(const point& p1)const{
    15          return !(sgn(x - p1.x) == 0 && sgn(y - p1.y) == 0 && sgn(z - p1.z) == 0);
    16      }
    17      point operator-(const point& p1)const{
    18          return point(x - p1.x, y - p1.y, z - p1.z);
    19      }
    20      point operator+(const point& p) const{
    21          return point(x + p.x, y + p.y, z + p.z);      
    22      }
    23      double operator^(const point& p1)const{
    24          return x * p1.x + y * p1.y + z * p1.z;
    25      }
    26      point operator*(const point& p1)const{
    27          return point(y * p1.z - z * p1.y, z * p1.x - x * p1.z, x * p1.y - y * p1.x);
    28      }
    29      point operator*(const double& l)const{
    30          return point(x * l, y * l, z * l);      
    31      }
    32 } p[101];
    33 
    34 
    35 int parallel(const point& u1,const point& u2,const point& v1,const point& v2){ //直线平行
    36     return sgn(((u1 - u2) * (v1 - v2)).len()) == 0;
    37 }
    38 
    39 double line2line_dist(const point& u1, const point& u2,const point& v1,const point& v2){ //计算直线u与v距离
    40     point n =  (u1 - u2) * (v1 - v2);
    41     return fabs((u1-v1) ^ n) / n.len();
    42 }
    43 
    44 point intersection(const point& u1,const point& u2, const point& v1,const point& v2){ //直线相交,需提前判断无解情况
    45     double t=((u1.x-v1.x)*(v1.y-v2.y) - (u1.y-v1.y)*(v1.x-v2.x)) 
    46               /((u1.x-u2.x)*(v1.y-v2.y) - (u1.y-u2.y)*(v1.x-v2.x));
    47     point p = u1 + (u2 - u1) * t;
    48     return p;
    49 }

    大概最近用到的就这些了。。。

  • 相关阅读:
    Java——HashMap
    Java——Collections
    Java——String,StringBuffer,StringBuilder
    OpenModelica读取文件
    Excel 常用设置
    Ubuntu16 源码方式安装postgresql数据库
    Utunbu常见问题
    PostgreSQL 扩展开发基础教程
    《数据挖掘导论》
    高性能MySQL --- 读书笔记(2)
  • 原文地址:https://www.cnblogs.com/yzcstc/p/3742610.html
Copyright © 2020-2023  润新知