首页 > 代码库 > 计算几何模板

计算几何模板

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     double ang = 0;
 3     for (int i = 0; i < n; ++i) {
 4         point p1 = pol[i] - p, p2 = pol[(i + 1) % n] - p;
 5         double c = (p1 ^ p2) / (p1.len() * p2.len());
 6         to_normal(c);
 7         ang += sgn(p1 * p2) * acos(c);
 8     }
 9     ang = abs(ang);
10     return ang < 0.5 * pi ? -1 : (ang < 1.5 * pi ? 0 : 1);
11 }

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\n", 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 }

 

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