计算几何模板(入门级)
尽管才刚入门 然而最近各方面能力都在下降 还是先把这个计算几何入门级模板放上来 弃坑做点其他的
----------------------------------------------------------------------------------------------------------------------
UPD0 16.3.26
UPD1 16.4.3 新增 向量旋转 两圆求交点 余弦定理 以及一些hint
UPD2 16.4.13 新增 轴对称点
UPD3 16.7.18 新增 三角形面积 凸包上最大三角形
UPD4 16.7.27 新增 三点确定平面一般方程 点到平面距离
1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <algorithm> 5 using namespace std; 6 const double eps = 1e-8, inf = 1e8, pi = acos(-1); 7 const int N = 110; 8 struct point 9 { 10 double x, y; 11 point(){} 12 point(double _x, double _y) 13 { 14 x = _x; 15 y = _y; 16 } 17 point operator - (const point &p) const 18 { 19 return point(x - p.x, y - p.y); 20 } 21 point operator + (const point &p) const 22 { 23 return point(x + p.x, y + p.y); 24 } 25 double operator * (const point &p) const 26 { 27 return x * p.y - y * p.x; 28 } 29 double operator / (const point &p) const 30 { 31 return x * p.x + y * p.y; 32 } 33 }a[N], b[N << 1]; 34 35 struct line 36 { 37 point p1, p2; 38 double ang; 39 line(){} 40 line(point p01, point p02) 41 { 42 p1 = p01; 43 p2 = p02; 44 } 45 void getang() 46 { 47 ang = atan2(p2.y - p1.y, p2.x - p1.x); 48 } 49 }c[N], q[N]; 50 51 struct circle 52 { 53 double x, y, r; 54 }cc; 55 56 struct point3 57 { 58 double x, y, z; 59 point3(){} 60 point3(double _x, double _y, double _z) 61 { 62 x = _x; 63 y = _y; 64 z = _z; 65 } 66 point3 operator - (const point3 &p) const 67 { 68 return point3(x - p.x, y - p.y, z - p.z); 69 } 70 point3 operator + (const point3 &p) const 71 { 72 return point3(x + p.x, y + p.y, z - p.z); 73 } 74 point3 operator * (const point3 &p)const 75 { 76 return point3(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x); 77 } 78 double operator / (const point3 &p)const 79 { 80 return x * p.x + y * p.y + z * p.z; 81 } 82 }; 83 84 int n; 85 86 bool cmp(const point &aa, const point &bb) 87 { 88 return aa.x < bb.x || (aa.x == bb.x && aa.y < bb.y); 89 } 90 91 bool cmpl(const line &aa, const line &bb) 92 { 93 if(abs(aa.ang - bb.ang) > eps) 94 return aa.ang < bb.ang; 95 return (bb.p2 - aa.p1) * (bb.p1 - bb.p2) > 0; 96 } 97 98 int sgn(double x) 99 //正负判断 100 { 101 if(x > eps) 102 return 1; 103 if(x < -eps) 104 return -1; 105 return 0; 106 } 107 108 double getdist2(const point &aa) 109 { 110 return (aa.x * aa.x + aa.y * aa.y); 111 } 112 113 bool online(const point &aa, const point &bb, const point &cc) 114 //三点共线判断 115 { 116 if(!(min(bb.x, cc.x) <= aa.x && aa.x <= max(bb.x, cc.x))) 117 return 0; 118 if(!(min(bb.y, cc.y) <= aa.y && aa.y <= max(bb.y, cc.y))) 119 return 0; 120 return !sgn((bb - aa) * (cc - aa)); 121 } 122 123 bool checkline(const point &aa, const point &bb, const point &cc, const point &dd) 124 //线段是否相交 125 { 126 int c1 = sgn((bb - aa) * (cc - aa)), c2 = sgn((bb - aa) * (dd - aa)), 127 c3 = sgn((dd - cc) * (aa - cc)), c4 = sgn((dd - cc) * (bb - cc)); 128 if(c1 * c2 < 0 && c3 * c4 < 0) 129 return 1;//规范相交 130 if(c1 == 0 && c2 == 0) 131 { 132 if(min(aa.x, bb.x) > max(cc.x, dd.x) || min(cc.x, dd.x) > max(aa.x, bb.x) 133 || min(aa.y, bb.y) > max(cc.y, dd.y) || min(cc.y, dd.y) > max(aa.y, bb.y)) 134 return 0; 135 return 1;//共线有重叠 136 } 137 if(online(cc, aa, bb) || online(dd, aa, bb) 138 || online(aa, cc, dd) || online(bb, cc, dd)) 139 return 1;//端点相交 140 return 0; 141 } 142 143 point getpoint(const point &aa, const point &bb, 144 const point &cc, const point &dd) 145 //两直线交点 146 { 147 double a1, b1, c1, a2, b2, c2; 148 point re; 149 a1 = aa.y - bb.y; 150 b1 = bb.x - aa.x; 151 c1 = aa * bb; 152 a2 = cc.y - dd.y; 153 b2 = dd.x - cc.x; 154 c2 = cc * dd; 155 //以下为交点横纵坐标 156 re.x = (c1 * b2 - c2 * b1) / (a2 * b1 - a1 * b2); 157 re.y = (a2 * c1 - a1 * c2) / (a1 * b2 - a2 * b1); 158 return re; 159 } 160 161 bool protrusion(point aa[]) 162 //判断多边形凹(凸) 163 { 164 int flag = sgn((aa[1] - aa[0]) * (aa[2] - aa[0])); 165 for(int i = 1; i < n; ++i) 166 if(sgn((aa[(i + 1) % n] - aa[i]) * (aa[(i + 2) % n] - aa[i])) != flag) 167 return 1; 168 return 0; 169 } 170 171 bool inpolygon(const point &aa) 172 //点在凸多边形内(外) 173 { 174 int flag = sgn((a[0] - aa) * (a[1] - aa)); 175 for(int i = 1; i < n; ++i) 176 if(sgn((a[i] - aa) * (a[(i + 1) % n] - aa)) != flag) 177 return 0; 178 return 1; 179 } 180 181 void transline(line &aa, double dist) 182 //逆时针方向平移线段 183 { 184 double d = sqrt((aa.p1.x - aa.p2.x) * (aa.p1.x - aa.p2.x) + 185 (aa.p1.y - aa.p2.y) * (aa.p1.y - aa.p2.y)); 186 point ta; 187 ta.x = aa.p1.x + dist / d * (aa.p1.y - aa.p2.y); 188 ta.y = aa.p1.y - dist / d * (aa.p1.x - aa.p2.x); 189 aa.p2 = ta + aa.p2 - aa.p1; 190 aa.p1 = ta; 191 } 192 193 void convexhull(point aa[], point bb[]) 194 //凸包 195 { 196 int len = 0; 197 sort(aa, aa + n, cmp); 198 bb[len++] = aa[0];//bb数组要开2倍(防止出现直线) 199 bb[len++] = aa[1]; 200 for(int i = 2; i < n ;++i) 201 { 202 while(len > 1 && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0) 203 //若严格则加上等于 204 --len; 205 bb[len++] = aa[i]; 206 } 207 int t = len; 208 for(int i = n - 2; i >= 0; --i) 209 { 210 while(len > t && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0) 211 //同上 212 --len; 213 bb[len++] = aa[i]; 214 } 215 --len; 216 n = len; 217 } 218 219 bool checkout(const line &aa, const line &bb, const line &cc) 220 //检查交点是否在向量顺时针侧 221 { 222 point p = getpoint(aa.p1, aa.p2, bb.p1, bb.p2); 223 return (cc.p1 - p) * (cc.p2 - p) < 0;//如果不允许共线或算面积 则此处不取等 224 } 225 226 double halfplane(point aa[], line bb[]) 227 //半平面交 228 { 229 sort(bb, bb + n, cmpl); 230 int n2 = 1; 231 for(int i = 1; i < n; ++i) 232 { 233 if(bb[i].ang - bb[i - 1].ang > eps) 234 ++n2; 235 bb[n2 - 1]= bb[i]; 236 } 237 n = n2; 238 int front = 0, tail = 0; 239 q[tail++] = bb[0], q[tail++] = bb[1]; 240 for(int i = 2; i < n; ++i) 241 { 242 while(front + 1 < tail && checkout(q[tail - 2], q[tail - 1], bb[i])) 243 --tail; 244 while(front + 1 < tail && checkout(q[front], q[front + 1], bb[i])) 245 ++front; 246 q[tail++] = bb[i]; 247 } 248 while(front + 2 < tail && checkout(q[tail - 2], q[tail - 1], q[front])) 249 --tail; 250 while(front + 2 < tail && checkout(q[front], q[front + 1], q[tail - 1])) 251 ++front; 252 if(front + 2 >= tail) 253 return 0; 254 int j = 0; 255 for(int i = front; i < tail; ++i, ++j) 256 { 257 aa[j] = getpoint(q[i].p1, q[i].p2, q[(i != tail - 1 ? i + 1 : front)].p1, 258 q[(i != tail - 1 ? i + 1 : front)].p2); 259 } 260 double re = 0; 261 for(int i = 1; i < j - 1; ++i) 262 re += (aa[i] - aa[0]) * (aa[i + 1] - aa[0]); 263 return abs(re * 0.5); 264 } 265 266 double calipers(point aa[]) 267 //旋转卡壳 268 { 269 double re = 0; 270 aa[n] = aa[0]; 271 int now = 1; 272 for(int i = 0; i < n; ++i) 273 { 274 while((aa[i + 1] - aa[i]) * (aa[now + 1] - aa[i]) > 275 (aa[i + 1] - aa[i]) * (aa[now] - aa[i])) 276 now = (now == n - 1 ? 0 : now + 1); 277 re = max(re, getdist2(aa[now] - aa[i])); 278 } 279 return re; 280 } 281 282 line bisector(const point &aa, const point &bb) 283 //中垂线(正方形) 284 { 285 double mx, my; 286 mx = (aa.x + bb.x) / 2; 287 my = (aa.y + bb.y) / 2; 288 line cc; 289 cc.p1.x = mx - (aa.y - bb.y) / 2; 290 cc.p1.y = my + (aa.x - bb.x) / 2; 291 cc.p2 = aa + bb - cc.p1; 292 cc.getang(); 293 return cc; 294 } 295 296 point rotate(const point &p, double cost, double sint) 297 //逆时针向量旋转 298 { 299 double x = p.x, y = p.y; 300 return point(x * cost - y * sint, x * sint + y * cost); 301 } 302 303 void getpoint(circle c1, circle c2) 304 //已确保两圆有交点时求出两圆交点 305 { 306 long double dab = sqrt(getdist2(point(c1.x, c1.y) - 307 point(c2.x, c2.y))); 308 if(c1.r > c2.r) 309 swap(c1, c2); 310 long double cost = (c1.r * c1.r + dab * dab - c2.r * c2.r) / 311 (c1.r * dab * 2); 312 long double sint = sqrt(1 - cost * cost); 313 point re = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, sint); 314 re.x = c1.x + re.x * (c1.r / dab); 315 re.y = c1.y + re.y * (c1.r / dab); 316 point re2 = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, -sint); 317 re2.x = c1.x + re2.x * (c1.r / dab); 318 re2.y = c1.y + re2.y * (c1.r / dab); 319 } 320 321 double mycos(double B, double C, double A) 322 //余弦定理 给定边长 323 { 324 return (B * B + C * C - A * A) / (B * C * 2); 325 } 326 327 double mycos2(const point &aa, const point &bb, const point &cc) 328 //余弦定理 给定点坐标 329 { 330 double C2 = getdist2(aa - bb), A2 = getdist2(bb - cc), B2 = getdist2(cc - aa); 331 return (B2 + C2 - A2) / (sqrt(B2 * C2) * 2); 332 } 333 334 point mirror_point(const point &aa, const point &bb, const point &cc) 335 //轴对称点 也可用来求垂足 336 { 337 double cost, sint; 338 cost = mycos2(bb, cc, aa); 339 sint = sqrt(1 - cost * cost); 340 if((cc - bb) * (aa - bb) > 0) 341 sint = -sint; 342 point re; 343 re = rotate(aa - bb, cost, sint) + bb; 344 re = rotate(re - bb, cost, sint) + bb; 345 return re; 346 } 347 348 double area(const point &aa, const point &bb, const point &cc) 349 //求三角形面积 350 { 351 return abs((bb - aa) * (cc - aa) / 2); 352 } 353 354 void max_triangle(const point aa[]) 355 //求凸包上最大三角形 356 { 357 double ans = 0; 358 int p1, p2, p3; 359 for(int i = 0; i < n; ++i) 360 { 361 int j = (i + 1) % n; 362 int k = (j + 1) % n; 363 while(k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n])) 364 k = (k + 1) % n; 365 if(k == i) 366 continue; 367 int kk = (k + 1) % n; 368 while(j != kk && k != i) 369 { 370 if(ans < area(aa[i], aa[j], aa[k])) 371 { 372 ans = area(aa[i], aa[j], aa[k]);//三角形面积 373 p1 = i; 374 p2 = j; 375 p3 = k; 376 } 377 while(k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n])) 378 k = (k + 1) % n; 379 j = (j + 1) % n; 380 } 381 } 382 point q1, q2, q3; 383 q1 = b[p2] + b[p3] - b[p1]; 384 q2 = b[p1] + b[p3] - b[p2]; 385 q3 = b[p1] + b[p2] - b[p3]; 386 //三角形上三个点 387 } 388 389 void get_panel(const point3 &p1, const point3 &p2, const point3 &p3, 390 double &A, double &B, double &C, double &D) 391 //由三点确定一个平面方程 392 { 393 A = ((p2.y - p1.y) * (p3.z - p1.z) - (p2.z - p1.z) * (p3.y - p1.y)); 394 B = ((p2.z - p1.z) * (p3.x - p1.x) - (p2.x - p1.x) * (p3.z - p1.z)); 395 C = ((p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x)); 396 D = (-(A * p1.x + B * p1.y + C * p1.z)); 397 } 398 399 double dist_panel(const point3 &pt, double A, double B, double C, double D) 400 //点到平面距离 401 { 402 return abs(A * pt.x + B * pt.y + C * pt.z + D) / 403 sqrt(A * A + B * B + C * C); 404 } 405 406 void hints() 407 { 408 //皮克定理 ans = s - point / 2 + 1 409 //atan2一定要先作差 再取模 再取劣弧 410 //精度不够时尝试用余弦定理替代三角函数 411 //算出sin cos tan中的一种后 其他的直接公式算出 412 /* 413 四面体重心 414 x1=(s1*x1+s2*x2+s3*x3+s4*x4)/(s1+s2+s3+s4); 415 y1=(s1*y1+s2*y2+s3*y3+s4*y4)/(s1+s2+s3+s4); 416 z1=(s1*z1+s2*z2+s3*z3+s4*z4)/(s1+s2+s3+s4); 417 */ 418 } 419 420 int main() 421 { 422 return 0; 423 }
