首页 > 代码库 > 计算几何模板
计算几何模板
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 }
大概最近用到的就这些了。。。
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。